C @(#)rebdec.for	19.1 (ESO-DMD) 02/25/03 13:31:23
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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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 REBDEC
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  program    REBDEC			version 1.00	   880930
C  Basic algorithm: L. Lucy		ESO - Garching
C  MIDAS implementation: D. Baade       ESO - Garching
C  K. Banse					1.10	   890517
C  1.20	    900220
C
C.KEYWORDS
C  (nonlinear) rebinning, deconvolution, point spread function, undersampling
C
C.PURPOSE
C  Rebin an image to smaller stepsize to remove effects of coarse PSF sampling 
C  ("pixelation"). Optionally deconvolve with PSD. 
C
C.ALGORITHM
C  Get the names of input + output frames from IN_A + OUT_A (both 2-D)
C  Do the rebinning (which is nonlinear in the way the flux is redistributed 
C  whereas the spatial coordinate transformation is strictly linear) while 
C  simultaneously deconvolving with user-supplied point spread function. 
C	
C.INPUT/OUTPUT
C  the following keys are used:
C
C  IN_A/C/1/60			input frame
C  IN_B/C/1/60			point spread function 
C  OUT_A/C/1/60			output frame
C  INPUTI/I/1/2			zoom factors in x and y
C  INPUTI/I/3/1			number of iterations 
C
C.VERSIONS
C  1.10		convert to Portable Midas
C  1.20		fix bug in STIPUT (missing NAXIS..)
C 
C 010201		last modif
C	
C --------------------------------------------------------------------
C
      IMPLICIT NONE
C      
      INTEGER    STAT,ZOOM(2),NITER,NAXIN,NPXIN(2)
      INTEGER    INNO,PSFNO,OUTNO,DUMNO,D0NO,D1NO,D2NO,D3NO,D4NO
      INTEGER    NAXPSF,NPXPSF(2),IAV,OFF(2)
      INTEGER    NPXPSR(4),NAXPSR,I,SUBLO(2),GS
      INTEGER    NPXOUT(2),RESPIX(2),SIZE,NAXOUT,NPIXL(2)
      INTEGER    NPIXZ(2),MSIZ,INPUTI(4)
      INTEGER    UNI(1),NULO,MADRID(1)
      INTEGER*8    INPTR,PSFPTR,OUTPTR,DUMPTR,D0PTR,D1PTR,D2PTR,
     +             D3PTR,D4PTR
C
      CHARACTER  INFRM*60,OUTFRM*60,PSF*60
      CHARACTER  CUNIN*64,CUNPSF*64,IDENT*72
      CHARACTER  IDENTP*72,CUNOUT*64,IDENTO*72
C      
      DOUBLE PRECISION   STAIN(2),STPIN(2),STAPSF(2)
      DOUBLE PRECISION   STPPSF(2),STPOUT(2),STPPSR(4)
      DOUBLE PRECISION   STAPSR(4),STAOUT(2),TEST,EPS
C 
      REAL       CUTS(2),CUTVALS(4)
C 
      INCLUDE    'MID_INCLUDE:ST_DEF.INC'
C 
      COMMON   /VMR/  MADRID
C 
      DATA    CUTVALS /4*0./,  SUBLO /1,1/
C      
      INCLUDE    'MID_INCLUDE:ST_DAT.INC'
C 
C  set up MIDAS environment + enable automatic error abort
      CALL STSPRO('REBDEC')
C
C  get name of input image  + map frame 
      CALL STKRDC('IN_A',1,1,60,IAV,INFRM,UNI,NULO,STAT)
      CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXIN,NPXIN,STAIN,STPIN,
     +            IDENT,CUNIN,INPTR,INNO,STAT)
      IF (NAXIN.NE.2) CALL STETER
     +   (1,'ERROR: Currently only 2-D frames supported!')
C      
C  get and check zoom factors in x + y ... and the number of iterations 
      CALL STKRDI('INPUTI',1,4,IAV,INPUTI,UNI,NULO,STAT)
      ZOOM(1) = INPUTI(1)
      ZOOM(2) = INPUTI(2)
      NITER = INPUTI(3)
      GS = INPUTI(4)
C
      IF ((((ZOOM(1)*ZOOM(2))/2)*2).EQ.(ZOOM(1)*ZOOM(2))) 
     +   CALL STETER
     +     (2,'ERROR: Both zoom factors need to be odd!')
      IF ((ZOOM(1).LT.1).OR.(ZOOM(2).LT.1))
     +   CALL STETER
     +     (3,'ERROR: Zoom factors less than unity are illegal!')
C
      IF (NITER.LT.1) CALL STETER 
     +     (4,'ERROR: Number of iterations must at least be 1!')
C      
C  get name of point spread function frame + map it
      CALL STKRDC('IN_B',1,1,60,IAV,PSF,UNI,NULO,STAT)
      CALL STIGET(PSF,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXPSF,NPXPSF,STAPSF,
     +            STPPSF,IDENTP,CUNPSF,PSFPTR,PSFNO,STAT)
C
C  Further error traps: make sure dimensions are the same ... 
      IF (NAXPSF.NE.NAXIN)
     +   CALL STETER
     +   (5,'Dimensions of inframe and PSF must be the same!')
C
C  ... steps have the right relative proportions, 
        IF ((NINT(STPIN(1)/STPPSF(1)).NE.ZOOM(1)).OR.
     +      (NINT(STPIN(2)/STPPSF(2)).NE.ZOOM(2)))
     +     CALL STETER
     +     (6,'Steps of inframe and PSF must differ by zoom fact!')
C
C  ... the zoom factor is not even when the number of pixels in the input PSF 
C  is odd (check both X and Y) 
      DO 100 I=1,2
         IF ((NPXPSF(I)/2)*2.EQ.NPXPSF(I)) 
     +      CALL STETER
     +        (7,'NPIX(PSF) must be odd in both X and Y!') 
100   CONTINUE
C
C  Get name of output frame, ...
         CALL STKRDC('OUT_A',1,1,60,IAV,OUTFRM,UNI,NULO,STAT)
         IF (OUTFRM.EQ.INFRM)             ! in- and output frame must be different
     +      CALL STETER
     +      (8,'Names of inframe + outframe must be different!')
C
C  ... depending on its existence or not, 
C
C  ... open it for reading and writing or create it for writing to it 
C
      IF (GS.EQ.0) THEN             ! we start without an initial guess
C
C  Calculate further defining parameters of output frame ...
         STPOUT(1)=STPIN(1)/ZOOM(1)
         STPOUT(2)=STPIN(2)/ZOOM(2)
         NPXOUT(1)=NPXIN(1)*ZOOM(1)
         NPXOUT(2)=NPXIN(2)*ZOOM(2)
C
C  The zooming
         STAOUT(1)=STAIN(1)-STPOUT(1)*(ZOOM(1)/2)
         STAOUT(2)=STAIN(2)-STPOUT(2)*(ZOOM(2)/2)
C      
C  ... and map it:
         IDENTO(1:) = IDENT(1:)
         CUNOUT(1:) = CUNIN(1:)
         NAXOUT = 2
         CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NAXOUT,NPXOUT,STAOUT,
     +               STPOUT,IDENT,CUNIN,OUTPTR,OUTNO,STAT)
C
      ELSE                        ! there is a guess
C
C  If there is a guess, we first map the frame ...
         CALL STIGET(OUTFRM,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
     +               2,NAXOUT,NPXOUT,STAOUT,
     +               STPOUT,IDENTO,CUNOUT,OUTPTR,OUTNO,STAT)
C
C  ... and then check if it is technically compatible with the 
C  input frame and the zoom factors:
C 
         IF (NAXOUT.NE.2) CALL STETER
     +      (9,'Output frame, too, must be 2-dimensional!')
         IF ((NPXIN(1)*NPXIN(2)*ZOOM(1)*ZOOM(2)).NE.
     +       (NPXOUT(1)*NPXOUT(2))) 
     +      CALL STETER
     +      (10,'NPIX(OUT) incompatible with NPIX(IN) + zoom factors!')
         TEST = (STAOUT(1)+STPOUT(1)*(ZOOM(1)/2))*
     +          (STAOUT(2)+STPOUT(2)*(ZOOM(2)/2))
         EPS = STPIN(1)*STPIN(2)*0.01
         IF (ABS(TEST-(STAIN(1)*STAIN(2))).GT.EPS)
     +      CALL STETER
     +      (11,'START values of in- and outframe differ!')
         EPS = STPIN(1)*STPIN(2)*0.0001
         IF (ABS((STPOUT(1)*ZOOM(1)*STPOUT(2)*ZOOM(2))-
     +       (STPIN(1)*STPIN(2))).GT.EPS)
     +      CALL STETER
     +      (12,'steps of in- and outframe must differ by factor ZOOM!')
      ENDIF                        ! guess or no guess?
C
C  We also need six working buffers for ... 
C
C  a) the output frame with extrapolation beyond its edges
      DO 400 I=1,2
C
C  The radius of the PSF determines the width of the margin to be added 
C  around inframe, for outframe additionally the zoom factor must be considered:
         OFF(I) = NPXPSF(I)/(2*ZOOM(I)) + 2
         NPIXL(I) = NPXIN(I)+2*OFF(I)
         NPIXZ(I) = NPIXL(I)*ZOOM(I)
400   CONTINUE
C 
      MSIZ = NPIXZ(1) * NPIXZ(2)
      CALL STFCRE('DUMMY',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            MSIZ,DUMNO,STAT)
      CALL STFMAP(DUMNO,F_X_MODE,1,MSIZ,IAV,DUMPTR,STAT)
C
C  If there is a guess, we transfer it right away to the enlarged frame:
      IF (GS.EQ.1) 
     +   CALL EXTND(MADRID(OUTPTR),MADRID(DUMPTR),
     +              NPXOUT(1),NPXOUT(2),NPIXZ(1),NPIXZ(2),OFF,
     +              ZOOM(1),ZOOM(2))
C
C  In order not to waste virtual memory space, we temporarily close OUTFRM
      CALL STFCLO(OUTNO,STAT)
C
C  b) the input frame after extrapolation beyond its edges
      MSIZ = NPIXL(1) * NPIXL(2)
      CALL STFCRE('DUMMY0',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            MSIZ,D0NO,STAT)
      CALL STFMAP(D0NO,F_X_MODE,1,MSIZ,IAV,D0PTR,STAT)
C
C  We do the enlargement here and then close inframe (no longer needed)
C 
      CALL EXTND (MADRID(INPTR),MADRID(D0PTR),NPXIN(1),NPXIN(2),
     +            NPIXL(1),NPIXL(2),OFF,1,1)
      CALL STFCLO(INNO,STAT)
C 
C  c) ratio of observed to calculated image
      MSIZ = NPIXL(1) * NPIXL(2)
      CALL STFCRE('DUMMY1',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            MSIZ,D1NO,STAT)
      CALL STFMAP(D1NO,F_X_MODE,1,MSIZ,IAV,D1PTR,STAT)
C
C  d) correction factors to be applied to previous intermediate result
      MSIZ = NPIXZ(1) * NPIXZ(2)
      CALL STFCRE('DUMMY2',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            MSIZ,D2NO,STAT)
      CALL STFMAP(D2NO,F_X_MODE,1,MSIZ,IAV,D2PTR,STAT)
C
C  e) re-calculated frame on coarse grid
      MSIZ = NPIXL(1) * NPIXL(2)
      CALL STFCRE('DUMMY3',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            MSIZ,D3NO,STAT)
      CALL STFMAP(D3NO,F_X_MODE,1,MSIZ,IAV,D3PTR,STAT)
C
C  f) PSF, but re-sampled on coarse grid for various centerings. 
C
C  First, the dimension of this new PSF need to be properly determined:
      NAXPSR = 4
      DO 1000 I=1,2
C
C  The number of pixels of the coarse PSF (see subroutine PSFRESAMP) must be 
C  odd. If ZOOM divides NPIX(PSF), this is automatically fulfilled because 
C  both numers are odd. 
         NPXPSR(2*I-1) = NPXPSF(I) / ZOOM(I)
C
C  Otherwise, we just make it odd:
         IF (((NPXPSF(I)/ZOOM(I))*ZOOM(I)).NE.NPXPSF(I))
     +      NPXPSR(2*I-1) = (NPXPSR(2*I-1)/2)*2 + 1
C
C  Finally, an additional margin of 1 is added at either limit in order 
C  to avoid truncations of input-PSF:
         NPXPSR(2*I-1) = NPXPSR(2*I-1)+2
         NPXPSR(2*I) = ZOOM(I)
         STPPSR(2*I-1) = STPPSF(I)
         STPPSR(2*I) = 1.D0/ZOOM(I)
         STAPSR(2*I-1) = STAPSF(I)
         STAPSR(2*I) = 1.D0/ZOOM(I)
1000  CONTINUE
C
      MSIZ = NPXPSR(1) * NPXPSR(2) * NPXPSR(3) * NPXPSR(4)
      CALL STFCRE('DUMMY4',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            MSIZ,D4NO,STAT)
      CALL STFMAP(D4NO,F_X_MODE,1,MSIZ,IAV,D4PTR,STAT)
C
C  Do the resampling of the PSF here and then close the input PSF 
C  because it's no longer used
      CALL PSFRES(NPXPSF(1),NPXPSF(2),NPXPSR(1),NPXPSR(3),
     +            ZOOM(1),ZOOM(2),MADRID(PSFPTR),MADRID(D4PTR))
      CALL STFCLO(PSFNO,STAT)
C      
C  Now comes the real thing:
      CALL REBDC(NPIXL(1),NPIXL(2),NPIXZ(1),NPIXZ(2),NITER,
     +           NPXPSR(1),NPXPSR(3),ZOOM(1),ZOOM(2),GS,
     +           MADRID(D0PTR),MADRID(D1PTR),MADRID(D2PTR),
     +           MADRID(D3PTR),MADRID(D4PTR),MADRID(DUMPTR))
C
C  Close dispensible intermediate files and re-open outframe . Cut off 
C  extrapolated regions of intermediate result as to match size of OUTFRM, 
C  insert into outframe:
C 
      CALL STFCLO(D0NO,STAT)
      CALL STFCLO(D1NO,STAT)
      CALL STFCLO(D2NO,STAT)
      CALL STFCLO(D3NO,STAT)
      CALL STFCLO(D4NO,STAT)
C 
      CALL STIGET(OUTFRM,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
     +            2,NAXOUT,NPXOUT,STAOUT,
     +            STPOUT,IDENTO,CUNOUT,OUTPTR,OUTNO,STAT)
      CALL TRUNC(MADRID(DUMPTR),MADRID(OUTPTR),NPIXZ(1),NPIXZ(2),
     +           NPXOUT(1),NPXOUT(2),OFF,ZOOM(1),ZOOM(2))
C
C  Last, we determine physical minimum and maximum of OUTFRM
      CALL STVALS('MIN',MADRID(OUTPTR),2,NPXOUT,SUBLO,
     +            NPXOUT,CUTVALS,CUTS,RESPIX,SIZE,STAT)
C
C  Copy descriptors and update LHCUTS + HISTORY
      CALL DSCUPT(INNO,OUTNO,' ',STAT)
      CALL STDWRR(OUTNO,'LHCUTS',CUTS,3,2,UNI,STAT)
C
C  Bye, Bye
      CALL STSEPI
      END

      SUBROUTINE REBDC(ICX,ICY,I1X,I1Y,NITER,OPSFX,OPSFY,ZX,ZY,GS,
     1                 PHT,RT,CF,PHI,PSF,PSI)
C
C      Removal of pixel grid from CCD data plus seeing by appropriate 
C      rebinning (including 2-D deconvolution) to smaller step size. 
C      The deconvolution algorithm follows the description given in 
C      The Astronomical Journal, Vol. 79, p. 745-754 (1974) but the 
C      rebinning to a smaller step is an extra feature. 
C
C      L.B. Lucy              ESO - Garching        1988 August
C      Modifications applied for MIDAS version (incl. generalisation 
C      to non-square images and PSF's):
C       
C      D. Baade               ESO - Garching        880930
C
      IMPLICIT NONE
C
      INTEGER   ICX,ICY,GS,ZX,ZY,NITER,STAT,I1X,I1Y,IT,I,J
      INTEGER   OPSFX,OPSFY
C
      REAL      PHT(ICX,ICY),PHI(ICX,ICY),CD,CD1,RT(ICX,ICY)
      REAL      PSF(OPSFX,ZX,OPSFY,ZY),CF(I1X,I1Y),PSI(I1X,I1Y)
C
      CHARACTER  TIME*25,CBUF*80
C
C     PHT        inframe, but extrapolated beyond edges, dims: ICX x ICY 
C     PHI        calculated frame in coarse grid
C     RT         basically, the ratio observed/calculated in coarse grid
C     CF         correction factor to be applied to current fine-grid result
C     PSI        fine-grid output frame, STAINg from some trivial guess
C     IC         dimensions of input frame: ICX x ICY
C     Z          zoom factors coarse -> fine grid: ZX, ZY
C     I1         dimensions of output frame: ZX*ICX x ZY*ICY
C
C
C ...      If an initial guess has not been provided, we start with a 
C ...      smooth trial-image
C
      CD = ZX*ZY
      IF (GS.EQ.0) THEN 
         CD1 = CD/(I1X*I1Y)
         DO 200 I=1,I1X                              
            DO 100 J=1,I1Y
               PSI(I,J) = CD1                  ! initial smooth trial - image
100         CONTINUE
200      CONTINUE
      ENDIF                                    ! start frame prepared 
C
      DO 1000 IT=1,NITER                       ! iteration loop
C
C  Form PSF-mapped coarse-grid image of guess PSI (result in frame PHI):
         CALL INTEGD(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY,
     +               PSF,PSI,PHI) 
C
C  Compute ratio of observed to calculated (guessed) image 
C
         DO 400 I=1,ICX
            DO 300 J=1,ICY
               RT(I,J) = 1.0                     ! initialize
               IF (PHI(I,J).NE.0.0) 
     +            RT(I,J) = CD*PHT(I,J)/PHI(I,J) ! obs/calc
300         CONTINUE
400      CONTINUE
C
C  Determine correction factor to current intermediate result
         CALL INTEGF(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY,
     +               PSF,RT,CF)  
C
C  Apply the correction, thereby completing this iteration
         DO 600 I=1,I1X
            DO 500 J=1,I1Y
               PSI(I,J) = PSI(I,J)*CF(I,J)       ! apply correc factor
500         CONTINUE
600      CONTINUE
C
C  The code is slow. Progress messages seem a good idea:
         CALL GENTIM(TIME)
         WRITE (CBUF,10000) TIME,IT
         CALL STTPUT(CBUF,STAT)
C
1000  CONTINUE
C
      RETURN
C 
10000 FORMAT (A,'  ',I2,'. iteration completed')
      END

      SUBROUTINE EXTND (RAW,PHT,IRX,IRY,ICX,ICY,OFF,ZX,ZY)
C
C     To avoid ugly effects at the frame limits, repeat the outermost 
C     columns (rows) RX (RY) times, i.e. by the halfwidth of the 
C     PSF resampled on the coarse grid.
C
      IMPLICIT NONE
C
      INTEGER     IRX,IRY,ICX,ICY,I,J,IQ,JQ,OFF(2),ZX,ZY,OFX,OFY
C
      REAL        RAW(IRX,IRY),PHT(ICX,ICY)
C
      OFX = OFF(1)*ZX
      OFY = OFF(2)*ZY
C
C  Insert raw frame at center of output frame 
C 
      DO 200 I=1,IRX
         IQ = I+OFX
         DO 100 J=1,IRY
            JQ = J+OFY
            PHT(IQ,JQ) = RAW(I,J)
100      CONTINUE
200   CONTINUE
C
C ...      Do the extrapolation, first in X:
C
      DO 500 I=1,OFX
         DO 300 J=1,IRY
            PHT(I,J+OFY) = RAW(1,J)
            PHT(IRX+OFX+I,J+OFY) = RAW(IRX,J)
300      CONTINUE
C
C  Here we also deal with the four corners
C
         DO 400 J=1,OFY
            PHT(I,J) = RAW(1,1)
            PHT(IRX+OFX+I,J) = RAW(IRX,1)
            PHT(I,IRY+OFY+J) = RAW(1,IRY)
            PHT(IRX+OFX+I,IRY+OFY+J) = RAW(IRX,IRY)
400      CONTINUE
500   CONTINUE
C
C  The continuation in Y:
C
      DO 700 I=1,IRX
         DO 600 J=1,OFY
            PHT(I+OFX,J) = RAW(I,1)
            PHT(I+OFX,IRY+OFY+J) = RAW(I,IRY)
600      CONTINUE
700   CONTINUE
C
      RETURN
      END

      SUBROUTINE PSFRES (IPSFX,IPSFY,OPSFX,OPSFY,ZX,ZY,PSFIN,PSF)
C
C       Author: L. Lucy                   ESO - Garching     1988 August
C       Modifications:        D. Baade    ESO - Garching     880930
C
C      The subroutine resamples the user-supplied high-resolution PSF 
C      on the coarser grid of the image to be deconvolved. The array 
C      holding the resulting PSF is four-dimensional because the resampling 
C      is done ZX x ZY times where each time the input PSF is 
C      shifted by 1/ZX and/or 1/ZY as to cover all cases of the 
C      relative positioning of the two PSF grids. 
C
      IMPLICIT NONE
C
      INTEGER   II,JJ,IP,JP,K,L,ZX,ZY,INDX,INDY,IPSFX,IPSFY
      INTEGER   OPSFX,OPSFY,DX,DY,OFFX,OFFY,CX,CY
C
      REAL      PSFIN(IPSFX,IPSFY),PSFBUF,CD,PSF(OPSFX,ZX,OPSFY,ZY)
C
      CD = 1./(ZX*ZY)
      DX = (OPSFX*ZX-IPSFX)/2       ! the excess in width of the resampled PSF
      DY = (OPSFY*ZY-IPSFY)/2      ! is always an even number of small pixels. 
C
C  Shift by one small pixel in X ...
C
      DO 1000 II=1,ZX
         OFFX = II-(ZX+1)/2+DX
C
C  and/or Y
         DO 800 JJ=1,ZY
            OFFY = JJ-(ZY+1)/2+DY
C
C  For all coarse-grid pixels ... 
            DO 600 IP=1,OPSFX
               CX = (IP-1)*ZX
               DO 500 JP=1,OPSFY
                  CY = (JP-1)*ZY
                  PSFBUF = 0.                        ! initialize 
C
C  co-add all fine-grid pixels covered by them at the current offsets
                  DO 400 K=1,ZX
                     INDX = CX+K-OFFX
C
C  Skip fine-grid pixels for which the input-PSF is undefined:
                     IF ( (INDX.LT.1) .OR. (INDX.GT.IPSFX) ) GOTO 400
                     DO 300 L=1,ZY
                        INDY = CY+L-OFFY
                        IF ( (INDY.GE.1) .AND. (INDY.LE.IPSFY) )
     +                     PSFBUF = PSFBUF + PSFIN(INDX,INDY)
300                  CONTINUE
400               CONTINUE 
C 
                  PSF(IP,II,JP,JJ) = PSFBUF*CD
500            CONTINUE
600         CONTINUE
C 
800      CONTINUE
1000  CONTINUE
C
      RETURN
      END

      SUBROUTINE INTEGD(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY,
     +                  PSF,PSI,PHI)
C
C       Author: L. Lucy                   ESO - Garching     1988 August
C       Modifications: D. Baade           ESO - Garching     880930
C
C      This subroutine maps the current fine-grid guess, PSI, 
C      on to the coarse grid (result in array PHI). The mapping 
C      uses the stack of PSF's that have been constructed by 
C      re-sampling the input PSF onto the coarse grid, i.e. basically 
C      a convolution with the psf is done.
C
      IMPLICIT NONE
C
      INTEGER   ICX,ICY,I1X,I1Y,OPSFX,OPSFY,IMIN,IMAX,JMIN,JMAX
      INTEGER   ID,ID0,ID1,ID2,II,IP,JD,JD0,JD1,JD2,JP,JJ,ZX,ZY
      INTEGER   I,J,CX(10000),CY(10000),RX,RY,RX1,RY1,RX1I,RY1J
      INTEGER   AUXX(2000),AUXY(2000)
C
      REAL      PSF(OPSFX,ZX,OPSFY,ZY),PSI(I1X,I1Y),PHIBUF,PHI(ICX,ICY) 
C
      RX = OPSFX/2                        ! 'radius' of re-sampled PSF, 
      RY = OPSFY/2                        !  full diameter is 2*R+1
      RX1 = RX+1
      RY1 = RY+1
C
C  For all fine pixels, determine X and Y coordinates on the coarse grid.
C  Also, introduce auxiliary array to further save computing time.
C
      DO 100 I=1,I1X
         CX(I) = (I-1)/ZX + 1
         AUXX(I) = (CX(I)-1)*ZX
100   CONTINUE
      DO 200 J=1,I1Y
         CY(J) = (J-1)/ZY + 1
         AUXY(J) = (CY(J)-1)*ZY
200   CONTINUE
C
C  Proceed pixel by pixel on the coarse grid
C
      IMIN = RX1
      IMAX = ICX-RX
      JMIN = RY1
      JMAX = ICY-RY
      ID0 = 1-RX1*ZX
      JD0 = 1-RY1*ZY
C
      DO 2000 I=IMIN,IMAX                         ! pixel of coarse grid
         ID1=ID0+I*ZX
         ID2=(I+RX)*ZX
         RX1I=RX1-I
         DO 1800 J=JMIN,JMAX                   ! pixel of coarse grid
            JD1=JD0+J*ZY
            JD2=(J+RY)*ZY
            RY1J=RY1-J
C
C  Then sum over all relevant pixels in fine grid. the PSF serves as 
C  the weighting function. For a small pixel (II,JJ) within a given large 
C  pixel (I,J) the value PSF(IP,II,JP,JJ) is taken. IP and JP point 
C  to the coarse pixels contributing (via the PSF) to the currently 
C  treated coarse pixel. IP and JP are no absolute coordinates but 
C  rather are the distances of the contributing large pixels from the 
C  large pixel whose flux is just being computed. II and JJ are also 
C  relative indices and only range from 1 to ZX and ZY, respectively.
C
            PHIBUF = 0.                        ! initialize 
C
            DO 1600 ID=ID1,ID2                 ! relevant range in fine grid
               II=ID-AUXX(ID)                  ! position within large pixel
               IP=CX(ID)+RX1I                  ! position of large pixel in PSF
               DO 1500 JD=JD1,JD2              ! relevant range in fine grid 
                  JJ = JD-AUXY(JD)             ! position within large pixel
                  JP = CY(JD)+RY1J             ! position of large pixel in PSF
                  PHIBUF = PHIBUF+PSF(IP,II,JP,JJ)*PSI(ID,JD)  ! off we go
1500           CONTINUE
1600        CONTINUE
            PHI(I,J) = PHIBUF                  ! this coarse bin finished 
C
1800     CONTINUE
2000  CONTINUE
C
C  Now extrapolate solution achieved so far to edges of frame by 
C  constant extrapolation
C
      DO 3000 I=1,RX
         DO 2200 J=JMIN,JMAX
            PHI(I,J)=PHI(IMIN,J)
            PHI(IMAX+I,J)=PHI(IMAX,J)
2200     CONTINUE
         DO 2400 J=1,RY
            PHI(I,J)=PHI(IMIN,JMIN)
            PHI(IMAX+I,J)=PHI(IMAX,JMIN)
            PHI(I,JMAX+J)=PHI(IMIN,JMAX)
            PHI(IMAX+I,JMAX+J)=PHI(IMAX,JMAX)
2400     CONTINUE
3000  CONTINUE
      DO 4000 I =IMIN,IMAX
         DO 3300 J=1,RY
            PHI(I,J)=PHI(I,JMIN)
            PHI(I,JMAX+J)=PHI(I,JMAX)
3300     CONTINUE
4000  CONTINUE
C
      RETURN
      END

      SUBROUTINE INTEGF(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY,PSF,RT,CF)  
C
C       Author: L. Lucy                   ESO - Garching     1988 August
C       Modifications: D. Baade            ESO - Garching     880930
C
C  Here the correction factors (pixelwise on the fine grid, array CF) 
C  are computed as the PSF-mapped ratio (pixelwise on the coarse grid, 
C  array RT) of the observed image to the result of the last iteration 
C  after its re-sampling on the coarse grid. 
C
      IMPLICIT NONE 
C
      INTEGER   I,J,I1X,I1Y,ICX,ICY,II,ID1,ID2,IMIN,IMAX,JMIN,JMAX
      INTEGER   JJ,JD1,JD2,ID,IP,JD,JP,OPSFX,OPSFY,ZX,ZY,RX,RY
      INTEGER   CX(10000),CY(10000),AUXX1(10000),AUXY1(10000)
      INTEGER   AUXX2(10000),AUXY2(10000)
C
      REAL      PSF(OPSFX,ZX,OPSFY,ZY),CF(I1X,I1Y),RT(ICX,ICY)
      REAL      CFBUF
C
      RX=OPSFX/2                        ! 'radius' of coarse PSF, 
      RY=OPSFY/2                        ! full diameter is 2*R+1
C
      IMIN=RX*ZX+1
      IMAX=I1X-RX*ZX
      JMIN=RY*ZY+1
      JMAX=I1Y-RY*ZY
C
C  Determine for all fine pixels coordinates on coarse grid (taking 
C  into account that the input frame has been enlarged). 
C  Also, introduce auxiliary arrays to further reduce computing time. 
C
      DO 100 I=IMIN,IMAX
         CX(I)=(I-1)/ZX+1                  ! X-coordinates on coarse grid 
         AUXX1(I)=RX-CX(I)+1
         AUXX2(I)=(CX(I)-1)*ZX
100   CONTINUE
      DO 200 J=JMIN,JMAX
         CY(J)=(J-1)/ZY+1                  ! Y-coordinates on coarse grid
         AUXY1(J)=RY-CY(J)+1
         AUXY2(J)=(CY(J)-1)*ZY
200   CONTINUE
C
      DO 1000 I=IMIN,IMAX                        ! pixel on fine grid
         II=I-AUXX2(I)                  ! pixel-# within coarse pixel
         ID1=CX(I)-RX
         ID2=CX(I)+RX
C
         DO 900 J=JMIN,JMAX                  ! pixel on fine grid
            JJ=J-AUXY2(J)                  ! pixel-# within coarse pixel
            JD1=CY(J)-RY
            JD2=CY(J)+RY
C
            CFBUF = 0.                     ! initialize 
C
            DO 800 ID=ID1,ID2                  ! relevant coarse-grid pixels 
               IP=AUXX1(I)+ID                  ! X-position within coarse PSF
               DO 700 JD=JD1,JD2                  ! relevant coarse-grid pixels 
                  JP=AUXY1(J)+JD            ! Y-position within coarse PSF
                  CFBUF=CFBUF+PSF(IP,II,JP,JJ)*RT(ID,JD)      ! the correction
700            CONTINUE
800         CONTINUE
            CF(I,J) = CFBUF
C
C  At this point the corrections for fine-grid pixel (I,J) to the 
C  ratios at positions (ID,JD) on the coarse grid have been fully 
C  evaluated for all pixels (ID,JD) concerned by pixel (I,J). 
C
900      CONTINUE
1000  CONTINUE
C
C  Fill up rest of frame by constant extrapolation of the edges
C
      DO 2000 I=1,RX*ZX
         DO 1100 J=JMIN,JMAX
            CF(I,J)=CF(IMIN,J)
            CF(IMAX+I,J)=CF(IMAX,J)
1100     CONTINUE
         DO 1200 J=1,RY*ZY
            CF(I,J)=CF(IMIN,JMIN)
            CF(IMAX+I,J)=CF(IMAX,JMIN)
            CF(I,JMAX+J)=CF(IMIN,JMAX)
            CF(IMAX+I,JMAX+J)=CF(IMAX,JMAX)
1200     CONTINUE
2000  CONTINUE
      DO 3000 I =IMIN,IMAX
         DO 2500 J=1,RY*ZY
            CF(I,J)=CF(I,JMIN)
            CF(I,JMAX+J)=CF(I,JMAX)
2500     CONTINUE
3000  CONTINUE
C
      RETURN
      END

      SUBROUTINE TRUNC(PSI,OUT,I1X,I1Y,IOX,IOY,OFF,ZX,ZY) 
C
C  This subroutine truncates the result frame such as to correspond 
C  to the true dimensions (in world coordinates) of the input frame.
C
      IMPLICIT NONE
C
      INTEGER   I1X,I1Y,IOX,IOY,I,J,IQ,JQ,OFF(2),ZX,ZY,OFX,OFY
C
      REAL      PSI(I1X,I1Y),OUT(IOX,IOY)
C
      OFX = OFF(1)*ZX
      OFY = OFF(2)*ZY
C
      DO 100 I=1,IOX
         IQ=I+OFX
         DO 50 J=1,IOY
            JQ = J+OFY
            OUT(I,J) = PSI(IQ,JQ)
50       CONTINUE
100   CONTINUE
C
      RETURN
      END
