C @(#)rfotresid.for	19.1 (ES0-DMD) 02/25/03 13:30:15
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
