C @(#)rfotresid.for 17.1.1.1 (ES0-DMD) 01/25/02 17:18:17 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C PROGRAM RESIDU C+++ C.IDENTIFICATION: RFOTRESID C.PURPOSE: Compute reconstructed image and the difference with the original C.AUTHOR: R. Buonanno, G. Buscema, C. Corsi, I. Ferraro, G. Iannicola C Osservatorio Astronomico di Roma C.VERSION: 881024 Creation R.B. C 890808 ST interfaces, MIDAS table file system implemented R.H.W. C---- IMPLICIT NONE C INTEGER NX INTEGER NY PARAMETER (NX=2048) PARAMETER (NY=2048) C INTEGER MADRID(1) INTEGER IAV, ISTAT INTEGER LX, LY INTEGER KUN, KNUL INTEGER NPL, NL INTEGER EC,ED, EL INTEGER IROW INTEGER TIDREG INTEGER IMFI,IMFA INTEGER IPI INTEGER*8 IPA INTEGER LD, LC INTEGER IX0, IY0 INTEGER I, J, K INTEGER NAXIS INTEGER NPIX(3) INTEGER TINULL INTEGER ICOL(12) INTEGER NCO,NRO,NSC,KW,NSA C DOUBLE PRECISION TDNULL,TDTRUE,TDFALS DOUBLE PRECISION BEGIN(3),STEP(3) C REAL VIV(NX,NY) REAL RIV(NX) REAL AL REAL DX2, DY2 REAL BETA, SIG, SI2 REAL TAB(12) REAL TRNULL REAL TBLSEL C CHARACTER IDENT*72,CUNIT*80 CHARACTER FRAMI*60,FRAMA*60 CHARACTER REGFIL*60 CHARACTER STRING*80 C LOGICAL SFLAG LOGICAL NUL(12) C INCLUDE 'MID_INCLUDE:TABLES.INC' COMMON /VMR/MADRID INCLUDE 'MID_INCLUDE:TABLED.INC' C DATA ICOL/2,3,4,5,6,7,8,9,10,11,12,13/ C C *** Hello, is MIDAS out there? CALL STSPRO('RESIDU') C C *** get the input frame CALL STKRDC('IN_A',1,1,60,IAV,FRAMI,KUN,KNUL,ISTAT) C C *** get the frame and read all relevant descriptors CALL STIGET(FRAMI,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, 2 3,NAXIS,NPIX,BEGIN,STEP,IDENT,CUNIT,IPI,IMFI,ISTAT) NPL = NPIX(1) NL = NPIX(2) C C create constructed frame CALL STKRDC('OUT_A',1,1,60,IAV,FRAMA,KUN,KNUL,ISTAT) IDENT = 'Reconstructed frame' CALL STIPUT(FRAMA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIX,BEGIN,STEP,IDENT, + CUNIT,IPA,IMFA,ISTAT) C C *** Input table CALL STKRDC('IN_B',1,1,60,IAV,REGFIL,KUN,KNUL,ISTAT) CALL STECNT('GET',EC,ED,EL) CALL STECNT('PUT',1,0,0) CALL TBTOPN(REGFIL,0,TIDREG,ISTAT) C C *** get the registration table opened CALL TBMNUL(TINULL,TRNULL,TDNULL) CALL TBMCON(TBLSEL,TDTRUE,TDFALS) CALL TBIGET(TIDREG,NCO,NRO,NSC,KW,NSA,ISTAT) IF (ISTAT.NE.0) THEN STRING = '*** FATAL: Problems with opening '// 2 'the registration file' CALL STTPUT(STRING,ISTAT) CALL STSEPI ENDIF CALL STECNT('PUT',EC,ED,EL) C AL = 4*ALOG(2.) C DO 100 IROW = 1,NRO CALL TBSGET(TIDREG,IROW,SFLAG,ISTAT) IF (SFLAG) THEN CALL TBRRDR(TIDREG,IROW,12,ICOL,TAB,NUL,ISTAT) BETA = TAB(10) IF (TAB(8).LT.0.) THEN SIG = TAB(9) SI2 = SIG**2 IF (BETA.GT.0) THEN LD = SI2*SQRT(ABS((.5/TAB(3))**(-1./BETA)-1))+.5 ELSE LD = SI2*SQRT(ABS(ALOG(2*TAB(3))))/AL+.5 END IF C C Convert from world to pixel coordinates C IX0 = ((TAB(1)-BEGIN(1))/STEP(1)+1.)-LD IY0 = ((TAB(2)-BEGIN(2))/STEP(2)+1.)-LD LX = 2*LD+1 LY = 2*LD+1 IF (IX0.LT.1) THEN LX = LX+IX0-1 IX0 = 1 END IF IF (IX0+LX-1.GT.NPL) THEN LX = NPL-IX0+1 ENDIF IF (IY0.LT.1) THEN LY = LY+IY0-1 IY0 = 1 END IF IF (IY0+LY-1.GT.NL) THEN LY=NL-IY0+1 ENDIF DO 110 J = IY0,LY+IY0 LC = 2 DY2 = (J-((TAB(2)-BEGIN(2))/STEP(2)+1.))**2 DO 120 K = IX0,LX+IX0 DX2 = (K-((TAB(1)-BEGIN(1))/STEP(1)+1.))**2 IF (BETA.GT.0) THEN VIV(K,J) = VIV(K,J)+ 2 (TAB(3)*(1+(DX2+DY2)/SI2)**(-BETA)) ELSE VIV(K,J) = VIV(K,J)+ 2 (TAB(3)*EXP(-AL*(DX2+DY2)/SI2)) END IF 120 CONTINUE LC = 1 110 CONTINUE END IF ENDIF 100 CONTINUE C DO 200 I = 1,NL DO 210 J = 1,NPL RIV(J) = VIV(J,I) 210 CONTINUE CALL WRILIN(NPL,NL,I,1,NPL,MADRID(IPA),RIV) 200 CONTINUE C CALL STDCOP(IMFI,IMFA,3,' ',ISTAT) CALL STSEPI END