C @(#)rctmap.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
      SUBROUTINE RCTMAP(OUT,NPX,NPY,START,IPOL,NPXIP,NPYIP,CX,CY,
     2                  OUTSTA,SUBFAC,NREP)
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
      INTEGER      NPX
      INTEGER      NPY
      INTEGER      NPXIP
      INTEGER      NPYIP
      REAL         OUT(NPX,NPY)
      REAL         START(2)
      REAL         IPOL(NPXIP,NPYIP)
      DOUBLE PRECISION CX(20),CY(20)
      REAL         OUTSTA(2)
      REAL         SUBFAC
      INTEGER      NREP
C
      DOUBLE PRECISION X,Y
      DOUBLE PRECISION X0,Y0
      REAL         WIDTH
      REAL         OFFX,OFFY
      INTEGER      IXR, IYR 
      INTEGER      ISTAT
      INTEGER      K, L
      INTEGER      M, N
      INTEGER      NPERC
      INTEGER      OFFINT
      REAL         XS,YS
      REAL         XR,YR
      REAL         FXIN,FXOUT,FYIN,FYOUT
      REAL         OFFMAX
      REAL         FLUX
      LOGICAL      CENTX
      LOGICAL      CENTY
      LOGICAL      LOWX
      LOGICAL      HIGHX
      LOGICAL      LOWY
      LOGICAL      HIGHY
      CHARACTER*50 PROGR
C
 9000 FORMAT ('          ',I3,' percent completed ...')
C  
      XS     = START(1) - 1                                     ! permit work on
      YS     = START(2) - 1                                          ! 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 *** transformation in x
C     ------------------------------------------------------------------------
C     |  note the CONVENIENT definition of a 2D third-order regression       |
C     |  analysis with two independent variables adopted                     |
C     |  in the midas command regression/poly. conformance to his definition |
C     |  is, however, mandatory.                                             |
C     ------------------------------------------------------------------------
C
                  XR     = CX(1)+X*(CX(2) +X*(CX(3) +X*CX(4))) +
     +                    (CX(5)+X*(CX(6) +X*(CX(7) +X*CX(8))))*Y +
     +                    (CX(9)+X*(CX(10)+X*(CX(11)+X*CX(12))))*
     +                     Y*Y  +  (CX(13)+X*(CX(14)+X*(CX(15)+X*
     +                     CX(16))))*Y*Y*Y
C
C *** and the same in y
C
                  YR     = CY(1)+X*(CY(2) +X*(CY(3) +X*CY(4))) +
     +                    (CY(5)+X*(CY(6) +X*(CY(7) +X*CY(8))))*Y +
     +                    (CY(9)+X*(CY(10)+X*(CY(11)+X*CY(12))))*
     +                     Y*Y  +  (CY(13)+X*(CY(14)+X*(CY(15)+X*
     +                     CY(16))))*Y*Y*Y
C
C *** transformation from world coordinates to pixel coordinates 
C     to be used throughout the rest of this subroutine
C
                  XR     = XR - OUTSTA(1)
                  YR     = YR - OUTSTA(2)
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
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
