C @(#)rectspec.for 17.1.1.1 (ES0-DMD) 01/25/02 17:54:00 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 17.1.1.1 (ESO) 01/25/02 17:54:00 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(1) 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