C @(#)rectspec.for	19.1 (ES0-DMD) 02/25/03 14:25: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 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 @(#)rectspec.for	19.1  (ESO)  02/25/03  14:25:23 
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C.IDENTIFICATION RECTSPEC
C.PURPOSE:       Geometrically rectificy distorted 2-D spectra, rebin to 
C                constant step size in wavelength maintaining the number of 
C                pixels, take care that flux is conserved
C.AUTHOR:        D. Baade, J.D.Ponz
C.KEYWORDS:      geometrical rectification of 2-D spectra, rebinning, flux
C                conservation, wavelength calibration
C.LANGUAGE:      F77+ESOext
C.ALGORITHM:     a) rectification: The parameters of a bipolynomial regression
C                   of variable degree have to be calculated first by 
C                   REGRESSION/POLY on TABLE data (e.g. with positions of 
C                   reseaux marks). They are stored in the keyword COEFXj 
C                   (x-transformation) and COEFYj (y-transformation), where
C                   j = I, C, R, D to keep all the regression information.
C                b) rebinning:     Each pixel is routinely subdivided into 9 
C                   subpixels (unless SUBF=1.). Assuming that
C                   the coordinates of the input pixel refer to its center
C                   and that its contents is the average flux integrated over
C                   its area (no dead space), a simple linear interpolation
C                   with nearest neighboring original pixels is done for each
C                   of the 9 subpixels. After this interpolation the flux
C                   summed over the subpixels is scaled to the flux of the
C                   original pixel so that flux conservation is achieved on a
C                   pixel basis. Finally, each subpixel is imaged to the
C                   rectified coordinate frame and its contents added to the
C                   appropriate pixel(s) (of about the same size as in the input
C                   frame). Optionally (NREP.GT.1) each subpixel can first be
C                   subdivided into still smaller units. Sub-subpixels
C                   originating from the same parent-subpixel all contain the
C                   same flux, i.e. no further interpolation is done.
C.RESTRICTIONS:  All STEP sizes are implicitly assumed to be positive.
C.INPUT/OUTPUT:  the following keywords are used:
C                P1    name of input frame
C                P2    name of output frame
C                P3    reference frame or sampling domain 
C                COEFXj  regression coefficients for X-coordinate
C                COEFYj  regression coefficients for Y-coordinate
C                INPUTI/I/1/1   repetition factor (for the further 
C                                 subdivision of subpixel
C                INPUTR/R/1/1     SUBF
C.VERSION:       861010  D. Baade                  ESO
C.VERSION:       871123  R.H. Warmels              ESO-FORTRAN Conversion
C.VERSION:       881214  J.D.Ponz                  New functionality
C --------------------------------------------------------------------
      PROGRAM     RCTSPC
C
      IMPLICIT    NONE
C
      INTEGER      NLEN
      PARAMETER    (NLEN=80)
      INTEGER      MADRID
      INTEGER      I,IAV,IAC,KUN,KNUL,COEFXI(13),COEFYI(13)
      INTEGER      NAXIS,NPIX(2),ISTAT
      INTEGER*8    IPNTRA,IPNTRC
      INTEGER      NREP,OUTNPX(2)
      INTEGER      IMA,IMC,IMD
      INTEGER*8    IPNTRD
      INTEGER      NPIXIP(2)
C 
      REAL         SUBF
C 
      DOUBLE PRECISION OUTSTR(2),OUTSTP(2)
      DOUBLE PRECISION INSTR(2),INSTP(2)
      DOUBLE PRECISION COEFXD(100),COEFYD(100),START(2),STEP(2)
C 
      CHARACTER*80 FRAMEA,FRAMEC,COEFXC,COEFYC
      CHARACTER*60 FRAMER
      CHARACTER    IDENT*72,CUNIT*64
C 
      INCLUDE     'MID_INCLUDE:ST_DEF.INC'
      COMMON      /VMR/MADRID(1)
      INCLUDE     'MID_INCLUDE:ST_DAT.INC'
C
C ... get into MIDAS
C
      CALL STSPRO('RECTIFY')  
C
C ... read parameters
C
      CALL STKRDC('P1',1,1,NLEN,IAC,FRAMEA,KUN,KNUL,ISTAT)
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .       2,NAXIS,NPIX,START,STEP,IDENT,CUNIT,IPNTRA,IMA,ISTAT)
C
      IF (NAXIS.NE.2) THEN                ! insist that frame be two-dimensional
         CALL STTPUT('*** FATAL: only 2-d frames allowed ... ',ISTAT)
         GO TO 1111
      END IF
      CALL STKRDC('P2',1,1,NLEN,IAC,FRAMEC,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_B',1,1,60,IAC,FRAMER,KUN,KNUL,ISTAT)
C
C ... read coefficients
C
      CALL STKRDC('KEYLONGC',1,1,NLEN,IAC,COEFXC,KUN,KNUL,ISTAT)
      CALL STKRDI('KEYLONGI',1,13,IAC,COEFXI,KUN,KNUL,ISTAT)
      CALL STKRDD('KEYLONGD',1,50,IAC,COEFXD,KUN,KNUL,ISTAT)

      CALL STKRDC('COEFYC',1,1,NLEN,IAC,COEFYC,KUN,KNUL,ISTAT)
      CALL STKRDI('COEFYI',1,13,IAC,COEFYI,KUN,KNUL,ISTAT)
      CALL STKRDD('COEFYD',1,50,IAC,COEFYD,KUN,KNUL,ISTAT)

       DO 66 I = 1,13
          COEFYI(I) = 0
66     CONTINUE

       DO 67 I = 1, 50
          COEFYD(I) = 0.D0
67     CONTINUE
C
C ... repetition factor and substepping flag
C
      CALL STKRDI('INPUTI',1,1,IAV,NREP,KUN,KNUL,ISTAT)
      CALL STKRDR('INPUTR',1,1,IAV,SUBF,KUN,KNUL,ISTAT)
      CALL STKRDD('INPUTD',1,2,IAV,INSTR,KUN,KNUL,ISTAT)
      CALL STKRDD('INPUTD',3,2,IAV,INSTP,KUN,KNUL,ISTAT)
C
      IF (NREP.GT.5) THEN                                      ! limit nrep to 5
         NREP   = 5
         CALL STTPUT('*** WARNING: Substepping factor set to 5 !',
     .               ISTAT)
      END IF
C
C ... calculate values of outstr and outnpx for output frame
C
      CALL OUTDSP(FRAMER,START,STEP,NPIX,
     .            COEFXI(6),COEFXI(7),COEFXD,
     .            OUTSTR,OUTSTP,OUTNPX)
C
      CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .   NAXIS,OUTNPX,OUTSTR,OUTSTP,IDENT,CUNIT,IPNTRC,IMC,ISTAT)
      CALL INTNUL(OUTNPX(1)*OUTNPX(2),MADRID(IPNTRC))  ! initialize its contents
      CALL STDCOP(IMA,IMC,4,'LHCUTS',ISTAT)       ! copy descriptor lhcuts
C
C ... test if deconvolution required
C
      IF (SUBF.GT.1.1) THEN                           
         DO 10 I = 1,2         
            NPIXIP(I) = NPIX(I)*SUBF
   10    CONTINUE
         CALL STIPUT('IPOLDUM',D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .    NAXIS,NPIXIP,START,STEP,IDENT,CUNIT,IPNTRD,IMD,ISTAT)
         CALL INTNUL(NPIXIP(1)*NPIXIP(2),MADRID(IPNTRD))
         CALL RCTSIN(MADRID(IPNTRA),NPIX(1),NPIX(2),
     .               MADRID(IPNTRD))
         CALL RCTSPE(MADRID(IPNTRC),OUTNPX(1),OUTNPX(2),START,
     .               MADRID(IPNTRD),NPIXIP(1),NPIXIP(2), 
     .               COEFXI(6),COEFXI(7),COEFXD,
     .               COEFYI(6),COEFYI(7),COEFYD,
     .               OUTSTR,OUTSTP,SUBF,NREP,INSTR,INSTP)
  
      ELSE                                                
C
C ... no deconvolution required
C
         CALL RCTSPE(MADRID(IPNTRC),OUTNPX(1),OUTNPX(2),START,
     .               MADRID(IPNTRA),NPIX(1),NPIX(2),
     .               COEFXI(6),COEFXI(7),COEFXD,
     .               COEFYI(6),COEFYI(7),COEFYD,
     .               OUTSTR,OUTSTP,SUBF,NREP,INSTR,INSTP)
      END IF
C
C ...  release files and update keywords
C
1111  CALL STSEPI
      END

      SUBROUTINE INTNUL(N, X)
C
C ... initialize to zero
C
      IMPLICIT NONE
      INTEGER N
      REAL    X(N)
C
      INTEGER I
C
      DO 10 I = 1, N
         X(I) = 0.
10    CONTINUE
      RETURN
      END


      SUBROUTINE RCTSPE(OUT,NPX,NPY,START,IPOL,NPXIP,NPYIP,
     .          KX,LX,CX,KY,LY,CY,OUTSTA,OUTSTP,SUBFAC,NREP,INSTR,INSTP)
C +++
C.PURPOSE: This subroutine maps each set of subpixels into the rectified
C          coordinate frame and adds the contents of each subpixel into
C          the output pixel concerned.
C ---
      IMPLICIT NONE
C
      INTEGER      NPXIP, NPYIP, NPX, NPY, KX, KY, LX, LY
      REAL         IPOL(NPXIP,NPYIP)
      REAL         OUT(NPX,NPY)
      DOUBLE PRECISION XS,YS,XR,YR
      DOUBLE PRECISION INSTR(2),INSTP(2)
      DOUBLE PRECISION START(2),WIDTH,OFFX,OFFY
      DOUBLE PRECISION FXIN,FXOUT,FYIN,FYOUT,OFFMAX,FLUX
      DOUBLE PRECISION OFFINT,OUTSTA(2),OUTSTP(2)
      REAL         SUBFAC
      DOUBLE PRECISION CX(20),CY(20),X,Y,X0,Y0, XR1, YR1
      LOGICAL      LOWX,CENTX,HIGHX,LOWY,CENTY,HIGHY
      INTEGER      IXR, IYR, NREP, NPERC, K, L, M, N, ISTAT
      CHARACTER*50 PROGR
C
 9000 FORMAT ('          ',I3,' percent completed ...')
C  
      XS     = 0.                                     ! permit work on
      YS     = 0.                                          ! subframes
      WIDTH  = 1./(SUBFAC*NREP)                        ! width of (sub-)subpixel
      OFFMAX = 0.5* (1-WIDTH)   ! max. offs. to be still in nearest output pixel
      OFFINT = INT(SUBFAC/2.)
      NPERC  = 0
      CALL STTPUT('*** INFO: Mapping started',ISTAT)
C
C *** map coordinates of subpixels into rectified grid
C
      DO 40 L = 1,NPYIP
         DO 30 K = 1,NPXIP
            X0     = XS + (K+OFFINT)/SUBFAC           ! coordinates of subpixels
            Y0     = YS + (L+OFFINT)/SUBFAC
            FLUX   = IPOL(K,L)/(NREP*NREP)
            DO 20 M = 1,NREP           ! if requested, go to still smaller units
               X      = X0 + (M-0.5-NREP*0.5)*WIDTH
               DO 10 N = 1,NREP
                  Y  = Y0 + (N-0.5-NREP*0.5)*WIDTH

C
C ... coefficients apply to wolrd coordinates
C
                  X = (X-1.D0)*INSTP(1)+INSTR(1)
                  Y = (Y-1.D0)*INSTP(2)+INSTR(2)
C
C ... transformation in x
C
                  CALL COMPOL(KX, LX, CX, X, Y, XR1)
C
C *** and the same in y
C
                  IF (KY.LE.0 .AND. LY.LE.0) THEN
                     YR1 = Y
                  ELSE
                     CALL COMPOL(KY, LY, CY, X, Y, YR1)
                  ENDIF
C
C *** transformation from world coordinates to pixel coordinates 
C     to be used throughout the rest of this subroutine
C
                  XR = (XR1 - OUTSTA(1))/OUTSTP(1) + 1
                  YR = (YR1 - OUTSTA(2))/OUTSTP(2) + 1
C
C *** drop flux into respective output pixel
C     if necessary, divide flux in one subpixel among several output
C     pixels assuming that the flux does not vary over one subpixel
C  
                  IXR    = NINT(XR)                     !   nearest output pixel
                  IYR    = NINT(YR)        !   (should receive most of the flux)
                  OFFX   = XR - IXR           !   offset of subpixel center from
                  OFFY   = YR - IYR                      !   output pixel center
                  IF (IXR.LE.1 .OR. IXR.GE.NPX .OR.
     .                IYR.LE.1 .OR. IYR.GE.NPY) GOTO 10
C
C *** determine position in x of subpixel relative to nearest output pixel
C
                  IF ((ABS(OFFX)-OFFMAX).LE.0) THEN     ! within cen. pixel in x
                     LOWX   = .FALSE.  
                     CENTX  = .TRUE.
                     HIGHX  = .FALSE.
                     FXIN   = 1.
                  ELSE IF ((XR-IXR-OFFMAX).GT.0) THEN  ! beyond cent. pixel in x
                     LOWX   = .FALSE.
                     CENTX  = .FALSE.  
                     HIGHX  = .TRUE.
                     FXIN   = (IXR+0.5-XR)/WIDTH + 0.5
                  ELSE                                ! below central pixel in x
                     LOWX   = .TRUE.
                     CENTX  = .FALSE.
                     HIGHX  = .FALSE.
                     FXIN   = (XR-IXR+0.5)/WIDTH + 0.5
                  END IF
                  FXOUT  = 1. - FXIN
C
C *** the same in y
C
                  IF ((ABS(OFFY)-OFFMAX).LE.0) THEN
                      LOWY   = .FALSE.
                      CENTY  = .TRUE.
                      HIGHY  = .FALSE.
                      FYIN   = 1.
                  ELSE IF ((YR-IYR-OFFMAX).GT.0) THEN
                      LOWY   = .FALSE.
                      CENTY  = .FALSE.
                      HIGHY  = .TRUE.
                      FYIN   = (IYR+0.5-YR)/WIDTH + 0.5
                  ELSE
                      LOWY   = .TRUE.
                      CENTY  = .FALSE.
                      HIGHY  = .FALSE.
                      FYIN   = (YR-IYR+0.5)/WIDTH + 0.5
                  END IF
                  FYOUT  = 1. - FYIN
C
C ***  First case, then: subpixel completely within one output pixel
                  IF (CENTX .AND. CENTY) THEN
                     OUT(IXR,IYR) = OUT(IXR,IYR) + FLUX
                  ELSE IF (HIGHX) THEN          ! test if subpixel on right edge
                     IF (HIGHY) THEN               ! pixel on upper right corner
C
C *** direct flux of subpixel into the four output pixels concerned
                        OUT(IXR+1,IYR+1) = OUT(IXR+1,IYR+1) +
     +                                     FXOUT*FYOUT*FLUX
                        OUT(IXR+1,IYR)   = OUT(IXR+1,IYR) +
     +                                     FXOUT*FYIN*FLUX
                        OUT(IXR,IYR+1)   = OUT(IXR,IYR+1) +
     +                                     FXIN*FYOUT*FLUX
                        OUT(IXR,IYR)     = OUT(IXR,IYR) +
     +                                     FXIN*FYIN*FLUX
                     ELSE IF (LOWY) THEN     ! three subp. on lower right corner
                        OUT(IXR+1,IYR-1) = OUT(IXR+1,IYR-1) +
     +                                     FXOUT*FYOUT*FLUX
                        OUT(IXR+1,IYR) = OUT(IXR+1,IYR) +
     +                                     FXOUT*FYIN*FLUX
                        OUT(IXR,IYR-1) = OUT(IXR,IYR-1) +
     +                                     FXIN*FYOUT*FLUX
                        OUT(IXR,IYR) = OUT(IXR,IYR) +
     +                                     FXIN*FYIN*FLUX
C
                     ELSE          ! four: subpixel on center right edge
C                                  ! (need to consider now only 2 output pixels)
                        OUT(IXR+1,IYR) = OUT(IXR+1,IYR) + FXOUT*FLUX
                        OUT(IXR,IYR)   = OUT(IXR,IYR) + FXIN*FLUX
                     END IF
                  ELSE IF (LOWX) THEN              ! subpixel falls on left edge
                     IF (HIGHY) THEN       ! five: subpixel on upper left corner
                        OUT(IXR-1,IYR+1) = OUT(IXR-1,IYR+1) +
     +                                     FXOUT*FYOUT*FLUX
                        OUT(IXR-1,IYR)   = OUT(IXR-1,IYR) +
     +                                     FXOUT*FYIN*FLUX
                        OUT(IXR,IYR+1)   = OUT(IXR,IYR+1) +
     +                                     FXIN*FYOUT*FLUX
                        OUT(IXR,IYR)     = OUT(IXR,IYR) +
     +                                     FXIN*FYIN*FLUX
C
                     ELSE IF (LOWY) THEN    ! six: subpixel on lower left corner
                        OUT(IXR-1,IYR-1) = OUT(IXR-1,IYR-1) +
     +                                     FXOUT*FYOUT*FLUX
                        OUT(IXR-1,IYR)   = OUT(IXR-1,IYR) +
     +                                     FXOUT*FYIN*FLUX
                        OUT(IXR,IYR-1)   = OUT(IXR,IYR-1) +
     +                                     FXIN*FYOUT*FLUX
                        OUT(IXR,IYR)     = OUT(IXR,IYR) +
     +                                      FXIN*FYIN*FLUX
C
                      ELSE          ! seventh case: subpixel on center left edge
                         OUT(IXR-1,IYR) = OUT(IXR-1,IYR) + FXOUT*FLUX
                         OUT(IXR,IYR)   = OUT(IXR,IYR) + FXIN*FLUX
                      END IF
C
                  ELSE      ! now only the center upper and lower edges are left
                     IF (LOWY) THEN       ! eigth: subpixel on center lower edge
                        OUT(IXR,IYR-1) = OUT(IXR,IYR-1) + FYOUT*FLUX
                        OUT(IXR,IYR) = OUT(IXR,IYR) + FYIN*FLUX

                     ELSE                  !nine: subpixel on center upper edge 
                        OUT(IXR,IYR+1) = OUT(IXR,IYR+1) + FYOUT*FLUX
                        OUT(IXR,IYR) = OUT(IXR,IYR) + FYIN*FLUX
                     END IF
                  END IF  
   10          CONTINUE                       !   this row of sub-subpixels done
   20       CONTINUE                       !   this column of sub-subpixels done
   30    CONTINUE                             !   this row of (sub-) pixels done
C
C *** feed them progress messages ...
C
         IF (((10*L)/NPYIP).GE. (NPERC+1)) THEN
            NPERC  = NPERC + 1
            WRITE (PROGR,9000) 10*NPERC
            CALL STTPUT(PROGR,ISTAT)
         END IF
C  
   40 CONTINUE                             !   this column of (sub-) pixels done
C
      RETURN
      END


      SUBROUTINE RCTSIN(IN,NPX,NPY,IPOL)
C +++
C.PURPOSE:     Interpolation
C.ALGORITHM:   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 ---
      IMPLICIT  NONE
C
      INTEGER   NPX, NPY
      REAL      IPOL(3*NPX,3*NPY)
      REAL      SUB(3,3),IN(*)
C
      INTEGER   NPERC,INDXX,INDXY
      INTEGER   ISTAT, K, L, I, J, IPOS
      CHARACTER PROGR*50

 9000 FORMAT ('*** INFO: ',I3,' percent completed ...')
C
      CALL STTPUT('*** INFO: Deconvolution started',ISTAT)
      NPERC  = 0
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(2)          !   first index: x, second: y
      SUB(3,2) =  SUB(3,1)
      SUB(3,3) = .5*IN(1) + .1667* (IN(2)+IN(NPX+1)+IN(NPX+2))
      SUB(2,3) = .667*IN(1) + .333*IN(NPX+1)
      SUB(1,3) =  SUB(2,3)
      CALL RCTSSC(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 RCTSSC(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 RCTSSC(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 RCTSSC(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 RCTSSC(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 RCTSSC(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 RCTSSC(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 RCTSSC(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 RCTSSC(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


      SUBROUTINE RCTSSC(SUB,IN)
C +++
C.PURPOSE: This subroutine compares the sum of the fluxes in the individual
C          subpixels to the flux in the original input pixel. The fluxes in
C          the subpixels are then scaled such that their sum equals the
C          original flux. Hence, flux conservation should be achieved on the
C          pixel scale.
C ---
      IMPLICIT  NONE
C
      REAL      SUB(3,3),IN,SUM,SCALE
C
      INTEGER   L, K
C
      SCALE  = 1
      SUM    = SUB(1,1) + SUB(1,2) + SUB(1,3) + SUB(2,1) + SUB(2,2) +
     +         SUB(2,3) + SUB(3,1) + SUB(3,2) + SUB(3,3)
C  
      IF (SUM.NE.0) THEN     ! assume sum does not vanish due to neg. fluxes ...
         SCALE  = IN/SUM 
      ENDIF
C
      DO 20 K = 1,3
         DO 10 L = 1,3
            SUB(K,L) = SUB(K,L)*SCALE
   10    CONTINUE
   20 CONTINUE
C
      RETURN
      END

      SUBROUTINE OUTDSP(REFIMA,START,STEP,NPIX,
     .                  K,L,CX,OUTSTA,OUTSTP,OUTPIX)
C +++
C.PURPOSE: Determine starting coordinates and stepsizes of output frame
C.ALGORITHM
C          Use the coordinates in REFIMA if specified, otherwise
C          use default values for the wavelength domain.
C          
C ---
      IMPLICIT NONE
C
      CHARACTER*60     REFIMA
      DOUBLE PRECISION START(2),STEP(2),OUTSTA(2),END(2),OUTSTP(2)
      DOUBLE PRECISION OUTEND(2),CORNER(4,2)
      DOUBLE PRECISION CX(20),X,Y,OUTX,WORK(2)
      INTEGER          NPIX(2),OUTPIX(2),K,L
      INTEGER          IMF,KUN,KNUL
C
      INTEGER  MADRID, I, IAC, ISTAT
      DOUBLE PRECISION RESULT
      DOUBLE PRECISION OUTXMN,OUTXMX
C
      INCLUDE    'MID_INCLUDE:ST_DEF.INC/NOLIST'
      COMMON     /VMR/MADRID(1)
      INCLUDE    'MID_INCLUDE:ST_DAT.INC/NOLIST'
C
      DATA     OUTXMN/1.E8/
      DATA     OUTXMX/-1.E8/
C
      IF (REFIMA(1:1) .EQ. '?') THEN
        DO 10 I = 1,2
           END(I) = START(I) + (NPIX(I)-1)*STEP(I)
   10   CONTINUE
C
        CORNER(1,1) = START(1)
        CORNER(1,2) = START(2)
        CORNER(2,1) = START(1)
        CORNER(2,2) = END(2)
        CORNER(3,1) = END(1)
        CORNER(3,2) = START(2)
        CORNER(4,1) = END(1)
        CORNER(4,2) = END(2)
C
C ... image corners to output space
C
        DO 20 I = 1,4
            X      = CORNER(I,1)
            Y      = CORNER(I,2)
            CALL COMPOL(K,L,CX,X,Y,RESULT)
            OUTX = RESULT
            IF (OUTX.GT.OUTXMX) OUTXMX = OUTX
            IF (OUTX.LT.OUTXMN) OUTXMN = OUTX
   20   CONTINUE
C            OUTX   = CX(1) + X * (CX(2)+X  * (CX(3)+X*CX(4))) +
C     +             (CX(5)+X  * (CX(6)+X  * (CX(7)+X*CX(8))))*Y +
C     +             (CX(9)+X  * (CX(10)+X * (CX(11)+X*CX(12))))*Y*Y +
C     +             (CX(13)+X * (CX(14)+X * (CX(15)+X*CX(16))))*Y*Y*Y
C
C ... find new steps
C
        OUTSTP(1) = (OUTXMX-OUTXMN)/NPIX(1)
        OUTSTP(2) = STEP(2)
C
C ... Deduce maximal dimensions of output frame
C
        OUTSTA(1) = OUTXMN - 5*OUTSTP(1)               ! safety margin
        OUTEND(1) = OUTXMX + 5*OUTSTP(1)               ! safety margin
        OUTSTA(2) = START(2)
        OUTPIX(1) = (OUTEND(1)-OUTSTA(1))/OUTSTP(1) + 1
        OUTPIX(2) = NPIX(2)
      ELSE
C
C ... get descriptors from reference image
C
        CALL STFOPN(REFIMA,D_R4_FORMAT,0,F_IMA_TYPE,IMF,ISTAT)
        CALL STDRDI(IMF,'NPIX',1,2,IAC,OUTPIX,KUN,KNUL,ISTAT)
        CALL STDRDD(IMF,'START',1,2,IAC,WORK,KUN,KNUL,ISTAT)
        OUTSTA(1) = WORK(1)
        OUTSTA(2) = WORK(2)
        CALL STDRDD(IMF,'STEP',1,2,IAC,WORK,KUN,KNUL,ISTAT)
        OUTSTP(1) = WORK(1)
        OUTSTP(2) = WORK(2)
      ENDIF
      RETURN
      END

      SUBROUTINE COMPOL(K,L,CX,X,Y,RESULT)
C +++
C.PURPOSE: Compute polynomial value
C ---
      IMPLICIT NONE
C
      DOUBLE PRECISION CX(20),X,Y,WORK(30)
      INTEGER  K,L
C
      INTEGER  IP, K1, L1
      DOUBLE PRECISION DC, RESULT
C
      IP     = 0
      DC     = 1.D0
      RESULT = 0.D0
      DO 20 L1 = 0, L
         IP       = IP + 1
         WORK(IP) = DC
         RESULT   = RESULT + WORK(IP)*CX(IP)
         DO 10 K1 = 1, K
            IP = IP + 1
            WORK(IP) = WORK(IP-1)*X
            RESULT   = RESULT + WORK(IP)*CX(IP)
10       CONTINUE
         DC = DC*Y
20    CONTINUE 
      RETURN
      END
