C @(#)rctmap.for 17.1.1.1 (ES0-DMD) 01/25/02 17:19:01 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