C @(#)rctint.for	19.1 (ES0-DMD) 02/25/03 13:30:47
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
C +++
C.PURPOSE:     Interpolation
C.ALGORYTHM:   Each input pixel is divided into 9 subpixels
C              Their indices are defined as follows:
C              1,3   2,3   3,3       FIRST INDEX:  X
C              1,2   2,2   3,2
C              1,1   2,1   3,1     SECOND INDEX: Y
C
C              For the interpolation, the following weighting scheme is used:
C             .1667       .1667      .333      .1667       .1667
C                      --------------------------------
C             .1667    |  .5         .667      .5     |    .1667
C             .333     |  .667      (1.)       .667   |    .333
C             .1667    |  .5         .667      .5     |    .1667
C                      --------------------------------
C             .1667       .1667      .333      .1667       .1667
C
C                (the box marks the edges of one input pixel)
C
C.NOTE:        This "deconvolution" makes an implicit assumption about
C              the point spread function, namely that its fwhm is smaller 
C              than or roughly equal to one input pixel. future improvements
C              could include reading in a more realistic or even exact 
C              point spread function and doing a proper deconvolution.
C ---
      SUBROUTINE RCTINT(IN,NPX,NPY,IPOL)
C
      IMPLICIT  NONE
      REAL      IN(1)
      INTEGER   NPX
      INTEGER   NPY
      REAL      IPOL(3*NPX,3*NPY)
C     
      REAL      SUB(3,3)
      INTEGER   NPERC
      INTEGER   INDXX
      INTEGER   INDXY
      CHARACTER PROGR*50
      INTEGER   ISTAT
      INTEGER   IPOS, TWO 
      INTEGER   I, J, K, L
C
 9000 FORMAT ('*** INFO: ',I3,' percent completed ...')
C
      CALL STTPUT('*** INFO: Deconvolution started',ISTAT)
      NPERC  = 0
      TWO = 2
C
C *** Begin with all pixels on the edge of the input frame
C *** lower left corner
C
      DO 20 K = 1,3
         DO 10 L = 1,3
            SUB(K,L) = IN(1)
   10    CONTINUE
   20 CONTINUE
C  
      SUB(3,1) = .667*IN(1) + .333*IN(TWO)        !   first index: x, second: y
      SUB(3,2) =  SUB(3,1)
      SUB(3,3) = .5*IN(1) + .1667* (IN(TWO)+IN(NPX+1)+IN(NPX+2))
      SUB(2,3) = .667*IN(1) + .333*IN(NPX+1)
      SUB(1,3) =  SUB(2,3)
      CALL RCTSCL(SUB,IN(1))
C
      DO 40 K = 1,3
         DO 30 L = 1,3
            INDXX  = K
            INDXY  = L
            IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L)
   30    CONTINUE
   40 CONTINUE
C
C *** lower right corner
C
      DO 60 K = 1,3
         DO 50 L = 1,3
            SUB(K,L) = IN(NPX)
   50    CONTINUE
   60 CONTINUE
C
      SUB(1,1) = .667*IN(NPX) + .333*IN(NPX-1)
      SUB(1,2) = SUB(1,1)
      SUB(1,3) = .5*IN(NPX) + .1667* (IN(NPX-1)+IN(2*NPX-1)+IN(2*NPX))
      SUB(2,3) = .667*IN(NPX) + .333*IN(2*NPX)
      SUB(3,3) = SUB(2,3)
C
      CALL RCTSCL(SUB,IN(NPX))
C
      DO 80 K = 1,3
         DO 70 L = 1,3
            INDXX  = NPX*3 - 3 + K
            INDXY  = L
            IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L)
   70    CONTINUE
   80 CONTINUE
C
C *** upper left corner
C
      IPOS  = NPX*(NPY-1) + 1
      DO 100 K = 1,3
         DO 90 L = 1,3
            SUB(K,L) = IN(IPOS)
   90    CONTINUE
  100 CONTINUE
C
      SUB(1,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX)
      SUB(2,1) = SUB(1,1)
      SUB(3,1) = .5*IN(IPOS) + .1667* (IN(IPOS-NPX)+IN(IPOS-NPX-1)+
     +           IN(IPOS+1))
      SUB(3,2) = .667*IN(IPOS) + .333*IN(IPOS+1)
      SUB(3,3) = SUB(3,2)
      CALL RCTSCL(SUB,IN(IPOS))
C
      DO 120 K = 1,3
         DO 110 L = 1,3
            INDXX  = K
            INDXY  = NPY*3 - 3 + L
            IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L)
  110    CONTINUE
  120 CONTINUE
C
C *** upper right corner
C
      IPOS   = NPX*NPY
      DO 140 K = 1,3
         DO 130 L = 1,3
            SUB(K,L) = IN(IPOS)
  130    CONTINUE
  140 CONTINUE
C
      SUB(1,3) = .667*IN(IPOS) + .333*IN(IPOS-1)
      SUB(1,2) =  SUB(1,3)
      SUB(3,1) = .5*IN(IPOS) + .1667* (IN(IPOS-1)+IN(IPOS-NPX-1)+
     +            IN(IPOS-NPX))
      SUB(2,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX)
      SUB(3,1) =  SUB(2,1)
      CALL RCTSCL(SUB,IN(IPOS))
C
      DO 160 K = 1,3
         DO 150 L = 1,3
            INDXX  = NPX*3 - 3 + K
            INDXY  = NPY*3 - 3 + L
            IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L)
  150    CONTINUE
  160 CONTINUE
C
C  *** bottom edge
C
      DO 210 I = 2,NPX - 1
         IPOS   = I
         DO 180 K = 1,3
            DO 170 L = 1,3
               SUB(K,L) = IN(IPOS)
  170       CONTINUE
  180    CONTINUE
C
         SUB(1,1) = .667*IN(IPOS) + .333*IN(IPOS-1)
         SUB(1,2) = SUB(1,1)
         SUB(1,3) = .5*IN(IPOS) + .1667*
     +              (IN(IPOS-1)+IN(IPOS+NPX-1)+IN(IPOS+NPX))
         SUB(2,3) = .667*IN(IPOS) + .333*IN(IPOS+NPX)
         SUB(3,3) = .5*IN(IPOS) + .1667*
     +              (IN(IPOS+1)+IN(IPOS+NPX)+IN(IPOS+NPX+1))
         SUB(3,2) = .667*IN(IPOS) + .333*IN(IPOS+1)
         SUB(3,1) = SUB(3,2)
         CALL RCTSCL(SUB,IN(IPOS))
C
         DO 200 K = 1,3
            DO 190 L = 1,3
               INDXX  = I*3 - 3 + K
               INDXY  = L
               IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L)
  190       CONTINUE
  200    CONTINUE
  210 CONTINUE
C
C *** top edge
C
      DO 260 I = 2,NPX - 1
         IPOS   = (NPY-1)*NPX + I
         DO 230 K = 1,3
            DO 220 L = 1,3
               SUB(K,L) = IN(IPOS)
  220       CONTINUE
  230    CONTINUE
C
         SUB(1,3) = .667*IN(IPOS) + .333*IN(IPOS-1)
         SUB(1,2) = SUB(1,3)
         SUB(1,1) = .5*IN(IPOS) + .1667*
     +              (IN(IPOS-1)+IN(IPOS-NPX-1)+IN(IPOS-NPX))
         SUB(2,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX)
         SUB(3,1) = .5*IN(IPOS) + .1667*
     +              (IN(IPOS-NPX)+IN(IPOS-NPX)+IN(IPOS-NPX+1))
         SUB(3,2) = .667*IN(IPOS) + .333*IN(IPOS+1)
         SUB(3,3) = SUB(3,2)
         CALL RCTSCL(SUB,IN(IPOS))
C
         DO 250 K = 1,3
            DO 240 L = 1,3
               INDXX  = I*3 - 3 + K
               INDXY  = NPY*3 - 3 + L
               IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L)
  240       CONTINUE
  250    CONTINUE
  260 CONTINUE
C
C *** left edge
C
      DO 310 J = 2,NPY - 1
         IPOS   = (J-1)*NPX + 1
         DO 280 K = 1,3
            DO 270 L = 1,3
               SUB(K,L) = IN(IPOS)
  270       CONTINUE
  280    CONTINUE
C
         SUB(1,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX)
         SUB(2,1) = SUB(1,1)
         SUB(3,1) = .5*IN(IPOS) + .1667*
     +              (IN(IPOS-NPX)+IN(IPOS-NPX+1)+IN(IPOS+1))
         SUB(3,2) = .667*IN(IPOS) + .333*IN(IPOS+1)
         SUB(3,3) = .5*IN(IPOS) + .1667*
     +              (IN(IPOS+1)+IN(IPOS+NPX+1)+IN(IPOS+NPX))
         SUB(2,3) = .667*IN(IPOS) + .333*IN(IPOS+NPX)
         SUB(1,3) = SUB(2,3)
         CALL RCTSCL(SUB,IN(IPOS))
C
         DO 300 K = 1,3
            DO 290 L = 1,3
               INDXX  = K
               INDXY  = J*3 - 3 + L
               IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L)
  290       CONTINUE
  300    CONTINUE
  310 CONTINUE
C
C *** right edge
C
      DO 360 J = 2,NPY - 1
         IPOS   = J*NPX
         DO 330 K = 1,3
            DO 320 L = 1,3
               SUB(K,L) = IN(IPOS)
  320       CONTINUE
  330    CONTINUE
C
         SUB(3,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX)
         SUB(2,1) = SUB(3,1)
         SUB(1,1) = .5*IN(IPOS) + .1667*
     +              (IN(IPOS-NPX)+IN(IPOS-NPX-1)+IN(IPOS-1))
         SUB(1,2) = .667*IN(IPOS) + .333*IN(IPOS-1)
         SUB(1,3) = .5*IN(IPOS) + .1667*
     +              (IN(IPOS-1)+IN(IPOS+NPX-1)+IN(IPOS+NPX))
         SUB(2,3) = .667*IN(IPOS) + .333*IN(IPOS+NPX)
         SUB(3,3) = SUB(2,3)
         CALL RCTSCL(SUB,IN(IPOS))
C
         DO 350 K = 1,3
            DO 340 L = 1,3
               INDXX  = NPX*3 - 3 + K
               INDXY  = J*3 - 3 + L
               IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L)
  340       CONTINUE
  350    CONTINUE
  360 CONTINUE
C
C *** and now the rest is straight forward
C
      DO 400 J = 2,NPY - 1
         DO 390 I = 2,NPX - 1
            IPOS   = (J-1)*NPX + I  !   CENTRAL SUBPIXEL
            SUB(2,2) = IN(IPOS)
            SUB(1,1) = .5*IN(IPOS) + .1667*
     +                 (IN(IPOS-1)+IN(IPOS-NPX-1)+IN(IPOS-NPX))
            SUB(2,1) = .667*IN(IPOS) + .333*IN(IPOS-NPX)
            SUB(3,1) = .5*IN(IPOS) + .1667*
     +                 (IN(IPOS-NPX)+IN(IPOS-NPX+1)+IN(IPOS+1))
            SUB(3,2) = .667*IN(IPOS) + .333*IN(IPOS+1)
            SUB(3,3) = .5*IN(IPOS) + .1667*
     +                (IN(IPOS+1)+IN(IPOS+NPX+1)+IN(IPOS+NPX))
            SUB(2,3) = .667*IN(IPOS) + .333*IN(IPOS+NPX)
            SUB(1,3) = .5*IN(IPOS) + .1667*
     +                 (IN(IPOS+NPX)+IN(IPOS+NPX-1)+IN(IPOS-1))
            SUB(1,2) = .667*IN(IPOS) + .333*IN(IPOS-1)
            CALL RCTSCL(SUB,IN(IPOS))
C
            DO 380 K = 1,3
               DO 370 L = 1,3
                  INDXX  = I*3 - 3 + K
                  INDXY  = J*3 - 3 + L
                  IPOL(INDXX,INDXY) = IPOL(INDXX,INDXY) + SUB(K,L)
  370          CONTINUE
  380       CONTINUE
  390    CONTINUE
C
C *** feed them progress messages ...
C
         IF (((10*J)/NPY).GE. (NPERC+1)) THEN
            NPERC  = NPERC + 1
            WRITE (PROGR,9000) 10*NPERC
            CALL STTPUT(PROGR,ISTAT)
         END IF
  400 CONTINUE

      RETURN
      END
