C @(#)rebin.for 17.1.1.1 (ESO-IPG) 01/25/02 17:41: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 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence 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 PROGRAM REBIN C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C program REBIN version 1.00 881028 C K. Banse ESO - Garching C from version 3.30 as of 871215 of VMS REBIN.FOR C C.KEYWORDS C linear rebinning C C.PURPOSE C rebin an image, i.e. calculate new gridpoints with new stepsizes C C.ALGORITHM C get the names of input + output frames from IN_A + OUT_A C scale everything down to unity, i.e. original grid will have distances = 1. C start point is 1,1 and end point nx,ny C if flux_consv = Y scale intensity of new pixels C to preserve total_signal/no._of_pixels ratio C (very crude thing - not really recommended) C C.INPUT/OUTPUT C the following keys are used: C C IN_A/C/1/80 input frame C IN_B/C/1/80 reference frame if given C OUT_A/C/1/80 output frame C INPUTD/D/1/6 new stepsize in x and y C "real" offsets in x and y (in wc) C translation values for start in x and y (in wc) C DEFAULT/C/1/5 (1) = Y for flux conservation, else = N C (2) = Y for absolute new start, C else old start + offset C (3) = qualif of command, S(pline) or L(inear) C (4) = processing option, R(ow), P(lane), all C (5) = input WC(world coordinates) or PIX(pixels) C C.VERSIONS C C 011115 last modif C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,NAXIS,MAPSIZ,MAXLIN INTEGER N,NK,NL,FELEI,FELEO,FLINI,FLINO,NOFFA,NOFFB INTEGER NLINES,NPROC,NZ,INSTA,IALL INTEGER NSAVE,MAXO,OPTIO,MXPROC,NKK INTEGER IMNOA,IMNOC,IMNOX,IMNOI INTEGER NSIZE(2),NULCNT,N1X,SWI INTEGER SIZE,SIZEA,SIZEC,STAT,IFLU INTEGER PLINCI,PLINCO,PLNSZI,PLNSZO INTEGER NPIXA(3),NPIXC(3),NPIXT(3),SUBLO(2) INTEGER RSMP,UNI(1),NULO,MADRID(1) C INTEGER*8 PNTRA,PSAVE,PNTRC,PNTRT,PNTRW,PNTRX,PNTRY,PNTRZ INTEGER*8 RESPTR C CHARACTER*80 FRAMEA,FRAMEC,REFFRA,CBUF CHARACTER CUNIT*64,IDENT*72 CHARACTER KFLUX*1,DEFAUL*5,PROCOP*1 C DOUBLE PRECISION STEPA(3),STARTA(3) ! old stepsize, old start DOUBLE PRECISION STEPC(3),STARTC(3) ! new stepsize, new start DOUBLE PRECISION INPUTD(6),OFFA(2),DDUM,ENDA(2),ZFRAC,XREST DOUBLE PRECISION STEP(2) ! stepsize in u.p.s. DOUBLE PRECISION STA(2),OFF(2),ASX,ASY C REAL RNULL,CUTS(4),CUTVAL(2),RFLU C COMMON /NULLC/ NULCNT,RNULL,CUTVAL COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA SUBLO /1,1/ DATA NPIXA /1,1,1/, NPIXC /1,1,1/ DATA STARTA /1.D0,1.D0,1.D0/, STARTC /1.D0,1.D0,1.D0/ DATA STEPA /1.D0,1.D0,1.D0/, STEPC /1.D0,1.D0,1.D0/ DATA STEP /1.0,1.0/ DATA IDENT /' '/, CUNIT /' '/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C set up MIDAS environment + enable automatic error abort CALL STSPRO('REBIN') NULCNT = 0 CUTVAL(1) = 999999.9 CUTVAL(2) = -999999.9 PLNSZI = 0 !increment between in-planes PLNSZO = 0 !same for result planes PLINCI = 0 !start increment for in-planes PLINCO = 0 !same for result planes RSMP = 0 C C get absolute_start, flux_conservation and rebin_method flag CALL STKRDC('DEFAULT',1,1,5,IAV,DEFAUL,UNI,NULO,STAT) KFLUX = DEFAUL(1:1) CALL STKRDR('NULL',2,1,IAV,RNULL,UNI,NULO,STAT) CALL STKRDI('MONITPAR',20,1,IAV,MAPSIZ,UNI,NULO,STAT) MAPSIZ = MAPSIZ * MAPSIZ !recommended buffer size PROCOP = DEFAUL(4:4) C C get input frame + open (or map) it CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT) IF (DEFAUL(3:3).NE.'S') THEN !linear rebinning CALL STFOPN(FRAMEA,D_R4_FORMAT,0,F_IMA_TYPE,IMNOA,STAT) CALL STDRDI(IMNOA,'NAXIS',1,1,IAV,NAXIS,UNI,NULO,STAT) CALL STDRDI(IMNOA,'NPIX',1,3,IAV,NPIXA,UNI,NULO,STAT) CALL STDRDD(IMNOA,'START',1,3,IAV,STARTA,UNI,NULO,STAT) CALL STDRDD(IMNOA,'STEP',1,3,IAV,STEPA,UNI,NULO,STAT) CALL STDRDC(IMNOA,'CUNIT',1,1,64,IAV,CUNIT,UNI,NULO,STAT) CALL STDRDC(IMNOA,'IDENT',1,1,72,IAV,IDENT,UNI,NULO,STAT) IF (NAXIS.GT.3) NAXIS = 3 ELSE CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXIS,NPIXA,STARTA,STEPA, + IDENT,CUNIT,PNTRA,IMNOA,STAT) IF (NAXIS.GT.3) NAXIS = 2 ENDIF C IF (NAXIS.EQ.3) THEN IF (PROCOP.EQ.'p') PROCOP = 'P' !process plane after plane IF (PROCOP.NE.'P') THEN CALL STTPUT('only 1. plane of cube is processed...',STAT) ELSE PLNSZI = NPIXA(1)*NPIXA(2) ENDIF STARTC(3) = STARTA(3) STEPC(3) = STEPA(3) NPIXC(3) = NPIXA(3) NAXIS = 2 ENDIF IF (NAXIS.EQ.2) THEN IF ((PROCOP.EQ.'R') .OR. (PROCOP.EQ.'r')) THEN PROCOP = 'R' !process row after row NAXIS = 1 !fake 1-dim input STEP(2) = 1.0D0 OFF(2) = 0.D0 STEPC(2) = STEPA(2) STARTC(2) = STARTA(2) ENDIF ELSE IF (NAXIS.EQ.1) THEN PROCOP = '+' ENDIF C C test for RESAMPLE/IMAGE IF (DEFAUL(3:3).EQ.'I') THEN CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT) STARTC(1) = STARTA(1) STARTC(2) = STARTA(2) CALL STKRDC('P3',1,1,80,IAV,REFFRA,UNI,NULO,STAT) IF ((REFFRA(1:1).EQ.'H').OR.(REFFRA(1:1).EQ.'h')) THEN RSMP = 1 STEPC(1) = STEPA(1) / 2 STEPC(2) = STEPA(2) / 2 NPIXC(1) = 2*NPIXA(1) - 1 NPIXC(2) = 2*NPIXA(2) - 1 ELSE RSMP = 2 STEPC(1) = STEPA(1) * 2.0D0 STEPC(2) = STEPA(2) * 2.0D0 NPIXC(1) = (NPIXA(1) + 1) / 2 NPIXC(2) = (NPIXA(2) + 1) / 2 ENDIF GOTO 1220 ENDIF C C look, if we use a reference frame CALL STKRDC('IN_B',1,1,80,IAV,REFFRA,UNI,NULO,STAT) IF (REFFRA(1:1).NE.'+') THEN CALL STFOPN(REFFRA,D_OLD_FORMAT,0,F_IMA_TYPE,IMNOX,STAT) CALL STDRDI(IMNOX,'NAXIS',1,1,IAV,N,UNI,NULO,STAT) IF (NAXIS.NE.N) CALL STETER + (4,'NAXIS must be the same as of input frame...') CALL STDRDD(IMNOX,'START',1,2,IAV,STARTC,UNI,NULO,STAT) CALL STDRDD(IMNOX,'STEP',1,2,IAV,STEPC,UNI,NULO,STAT) DO 500 N=1,NAXIS STEP(N) = STEPC(N)/STEPA(N) !stepsize in unit pixel space OFFA(N) = STARTC(N) - STARTA(N) OFF(N) = OFFA(N) / STEPA(N) !offset in unit pixel space 500 CONTINUE CALL STFCLO(IMNOX,STAT) ELSE C C get new stepsizes, offsets and translation values in x + y C get input like entered in MIDAS command line CALL STKRDD('INPUTD',1,6,IAV,INPUTD,UNI,NULO,STAT) C now let see, if P8 set to pixels PIX: IF (DEFAUL(5:5).EQ.'P') THEN C input in pixels instead of worlds coordinates DO 900 N=1,NAXIS C handle pixels: STEP(N) = INPUTD(N) ! new stepsize in pixels OFFA(N) = INPUTD(N+2) !"real" offset from old start ! this value would be still wc C C adjust for pixel center - new stepsize STEPC(N) = STEP(N) * STEPA(N) !stepsize is already in u.p.s. !stepc is later needed for output OFFA(N) = OFFA(N) + (0.5D0*(STEPC(N)-STEPA(N))) OFF(N) = OFFA(N) / STEPA(N) !offset in unit pixel space C IF (DEFAUL(2:2).EQ.'Y') THEN STARTC(N) = INPUTD(N+4) !absolute value for new start ELSE STARTC(N) = STARTA(N) + OFFA(N) ENDIF C ^^^^^^^ end of part with pixels 900 CONTINUE C ELSE C this is the default case: input in world coordinates DO 1000 N=1,NAXIS STEPC(N) = INPUTD(N) !new stepsize OFFA(N) = INPUTD(N+2) !"real" offset from old start C C adjust for pixel center - new stepsize OFFA(N) = OFFA(N) + (0.5D0*(STEPC(N)-STEPA(N))) OFF(N) = OFFA(N) / STEPA(N) !offset in unit pixel space STEP(N) = STEPC(N) / STEPA(N) !stepsize in u.p.s. IF (DEFAUL(2:2).EQ.'Y') THEN STARTC(N) = INPUTD(N+4) !absolute value for new start ELSE STARTC(N) = STARTA(N) + OFFA(N) ENDIF 1000 CONTINUE ENDIF C new end ENDIF C C calculate no. of points for output frame DO 1200,N=1,NAXIS IF (STEP(N) .LT. 1.0D-20) THEN WRITE(CBUF,10000) STEPC(N),N CALL STETER(1,CBUF) ENDIF C STA(N) = 0.5D0 + OFF(N) !shift 0 -> 1 origin (- half pixel) ENDA(N) = STARTA(N)+(NPIXA(N)-1)*STEPA(N) !end of input frame DDUM = ENDA(N) + 0.5*STEPA(N) - + ((STARTA(N)+OFFA(N)) - 0.5*STEPC(N)) !do not use STARTC... IF (STEP(N).GT.1.D0) THEN NPIXC(N) = INT(DDUM/STEPC(N)+0.0000001D0) ELSE NPIXC(N) = NINT(DDUM/STEPC(N)) ENDIF 1200 CONTINUE C C get output frame 1220 CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT) IF (FRAMEC.EQ.FRAMEA) !in- and output frame must be different + CALL STETER(2,'input + output frame must be different...') C C calculate different image sizes IF (PROCOP.EQ.'R') THEN NAXIS = 2 !reset NAXIS NPIXC(2) = NPIXA(2) ENDIF SIZEA = 1 SIZEC = 1 DO 2000 N=1,NAXIS IF (NPIXC(N).LT.1) THEN WRITE(CBUF,10001) N,NPIXC(N) CALL STETER(12,CBUF) ENDIF SIZEA = SIZEA * NPIXA(N) SIZEC = SIZEC * NPIXC(N) 2000 CONTINUE C IFLU = 0 RFLU = 0.0 IF (DEFAUL(3:3).NE.'I') THEN IF ((KFLUX.EQ.'Y') .OR. (KFLUX.EQ.'y')) THEN DDUM = STEP(1) * STEP(2) IF (DDUM .GT. 1.D0) THEN IFLU = 1 RFLU = (1.D0 / DDUM) ENDIF ENDIF ENDIF C C branch according to method (bilinear or cubic spline) IF (DEFAUL(3:3).NE.'S') GOTO 3900 C C C here for cubic spline rebinning C C CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT) !get old cuts C C map output file CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIXC,STARTC,STEPC, + IDENT,CUNIT,PNTRC,IMNOC,STAT) C ASX = STEP(1) !but for stepsize > 1.0 no splines...! IF (NAXIS.EQ.1) THEN ASY = ASX ELSE ASY = STEP(2) ENDIF IF ((ASX.GT.1.0D0).OR.(ASY.GT.1.0D0)) THEN CALL STTPUT + ('new step(s) larger than input => linear rebinning...',STAT) GOTO 3900 !no spline for those ENDIF C C map array for 2. derivative (row) and temporary space (row) SIZE = MAX(NPIXA(1),NPIXA(2)) CALL STFXMP(SIZE,D_R4_FORMAT,PNTRW,STAT) CALL STFXMP(SIZE,D_R4_FORMAT,PNTRT,STAT) C C if 2-dim frame, map intermed. arrays with interpolated points only rowwise IF (NPIXC(2) .GT. 1) THEN SIZE = NPIXC(1) * NPIXA(2) CALL STFXMP(SIZE,D_R4_FORMAT,PNTRZ,STAT) CALL STFXMP(SIZE,D_R4_FORMAT,PNTRX,STAT) C C we need also temporary result frame, because we have to transpose again SIZE = NPIXC(1) * NPIXC(2) CALL STFXMP(SIZE,D_R4_FORMAT,PNTRY,STAT) RESPTR = PNTRZ NPIXT(1) = NPIXC(1) NPIXT(2) = NPIXA(2) ELSE RESPTR = PNTRC ENDIF C C first do it along the rows NOFFA = 0 NOFFB = 0 C DO 3000 NL=1,NPIXA(2) C calculate the second derivative along the rows CALL DERIV2(MADRID(PNTRA),NOFFA,MADRID(PNTRW), + MADRID(PNTRT),NPIXA(1)) C and resample in x CALL RESMPX(MADRID(PNTRA),NOFFA,NPIXA(1),STA(1),STEP(1), + MADRID(PNTRW),MADRID(RESPTR),NPIXC(1),NOFFB) NOFFA = NOFFA + NPIXA(1) NOFFB = NOFFB + NPIXC(1) 3000 CONTINUE C IF (NPIXC(2).LE.1) GOTO 6000 C C if 2-dim frame we transpose and move along the former columns NSIZE(1) = 128 NSIZE(2) = 256 CALL LINCOL(MADRID(PNTRZ),NPIXT,NSIZE,MADRID(PNTRX)) C C now do it along the columns NOFFA = 0 NOFFB = 0 DO 3300 NL=1,NPIXC(1) C calculate the second derivative along the columns CALL DERIV2(MADRID(PNTRX),NOFFA,MADRID(PNTRW), + MADRID(PNTRT),NPIXA(2)) C and resample in y CALL RESMPX(MADRID(PNTRX),NOFFA,NPIXA(2),STA(2),STEP(2), + MADRID(PNTRW),MADRID(PNTRY),NPIXC(2),NOFFB) NOFFA = NOFFA + NPIXA(2) NOFFB = NOFFB + NPIXC(2) 3300 CONTINUE C C transpose back NPIXA(1) = NPIXC(2) NPIXA(2) = NPIXC(1) CALL LINCOL(MADRID(PNTRY),NPIXA,NSIZE,MADRID(PNTRC)) IF (IFLU.EQ.1) THEN CALL FLUFAC(MADRID(PNTRC),SIZEC,RFLU) DO 3500, N=1,4 CUTS(N) = RFLU * CUTS(N) 3500 CONTINUE ENDIF GOTO 6000 C C C here for bilinear rebinning C C 3900 IF (PLNSZI.GT.0) THEN NAXIS = 3 !reset input NAXIS PLNSZO = NPIXC(1)*NPIXC(2) SIZEC = SIZEC*NPIXC(3) ENDIF C CALL STFCRE(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, !create result + SIZEC,IMNOC,STAT) !frame CALL STDWRI(IMNOC,'NAXIS',NAXIS,1,1,UNI,STAT) CALL STDWRI(IMNOC,'NPIX',NPIXC,1,NAXIS,UNI,STAT) CALL STDWRD(IMNOC,'START',STARTC,1,NAXIS,UNI,STAT) CALL STDWRD(IMNOC,'STEP',STEPC,1,NAXIS,UNI,STAT) CALL STDWRC(IMNOC,'IDENT',1,IDENT,1,72,UNI,STAT) NL = (NAXIS+1)*16 CALL STDWRC(IMNOC,'CUNIT',1,CUNIT,1,NL,UNI,STAT) C IF (RSMP.GE.1) THEN CALL STFXMP(NPIXA(1),D_R4_FORMAT,PNTRA,STAT) CALL STFXMP(NPIXC(1),D_R4_FORMAT,PNTRX,STAT) !1 line of output CALL STFXMP(NPIXC(1),D_R4_FORMAT,PNTRY,STAT) CALL STFXMP(NPIXC(1),D_R4_FORMAT,PNTRZ,STAT) C CALL RESMPS(RSMP,IMNOA,IMNOC,MADRID(PNTRA),MADRID(PNTRX), + NPIXA,NPIXC,MADRID(PNTRY),MADRID(PNTRZ)) GOTO 6000 ENDIF C CALL STFXMP(NPIXC(1),D_R8_FORMAT,PNTRX,STAT) !XS array of RESMPA CALL STFXMP(NPIXC(1),D_I4_FORMAT,PNTRY,STAT) !LS array of RESMPA CALL STFXMP(NPIXC(1),D_I4_FORMAT,PNTRZ,STAT) !HS array of RESMPA CALL BUILDA(NPIXA(1),STA(1),STEP(1),NPIXC(1), + MADRID(PNTRX),MADRID(PNTRY),MADRID(PNTRZ),NKK) NULCNT = NKK * NPIXA(2) CUTS(1) = 0.0 CUTS(2) = 0.0 C C 1-dim processing C IF (NPIXA(2) .LE. 1) THEN CALL STFXMP(NPIXA(1),D_R4_FORMAT,PNTRA,STAT) CALL STFXMP(NPIXC(1),D_R4_FORMAT,PNTRC,STAT) CALL STFGET(IMNOA,1,NPIXA(1),IAV,MADRID(PNTRA),STAT) C CALL RESMPA(MADRID(PNTRA),NPIXA(1),STA(1),STEP(1),NPIXC(1), + MADRID(PNTRC),1,1,NKK, + MADRID(PNTRX),MADRID(PNTRY),MADRID(PNTRZ)) C CUTS(3) = CUTVAL(1) CUTS(4) = CUTVAL(2) !calculate also min,max CALL FLUFAK(IFLU,MADRID(PNTRC),NPIXC(1),RFLU,CUTS(3)) CALLSTFPUT(IMNOC,1,NPIXC(1),MADRID(PNTRC),STAT) C C 2-dim processing C ELSE IF (NPIXC(1).GT.NPIXA(1)) THEN !determine no. of lines N = NPIXC(1) !for all buffers ELSE !with different line length N = NPIXA(1) ENDIF NLINES = MAPSIZ / N IF (NPIXC(2).GT.NPIXA(2)) THEN !also depends on offset MAXLIN = NPIXC(2) ELSE MAXLIN = NPIXA(2) ENDIF IF (NLINES .GT. MAXLIN) NLINES = MAXLIN C INSTA = 1 IF (STEP(2).GE.1.D0) THEN DDUM = STA(2) - (STEP(2)*0.5D0) !go half a new stepsize back IF (DDUM.LE.0.D0) THEN NZ = 0 ZFRAC = 0.D0 ELSE NZ = INT(DDUM) ZFRAC = DDUM - NZ INSTA = NZ - 1 ENDIF N = (INT(STEP(2)) + 1) * 3 IF (NLINES .LT. N) NLINES = N !at least 3 steps must fit SWI = 1 ELSE INSTA = INT(STA(2)) IF (NLINES .LT. 10) NLINES = 10 SWI = 0 ENDIF IF (INSTA.LT.1) INSTA = 1 !1. row we have to work on C NSAVE = NLINES !save it for later NL = NLINES * NPIXC(1) CALL STFXMP(NL,D_R4_FORMAT,PNTRW,STAT) !out buffer CALL STFXMP(NL,D_R4_FORMAT,PNTRC,STAT) !out buffer NK = NLINES * NPIXA(1) CALL STFCRE('middummin',D_R4_FORMAT,F_X_MODE, !alloc memory + F_IMA_TYPE,NK,IMNOI,STAT) CALL STFMAP(IMNOI,F_X_MODE,1,NK,IAV,PNTRA,STAT) !in buffer C FLINI = INSTA IF ((FLINI+NLINES-1).GE.NPIXA(2)) THEN IALL = 1 !all data fits into int_buff PSAVE = PNTRW NLINES = NPIXA(2) ELSE IALL = 0 IF (PROCOP.NE.'R') THEN SIZE = NPIXC(1) * NPIXA(2) CALL STFCRE('middummbin.bdf',D_R4_FORMAT, !create tempfile in + F_O_MODE,F_IMA_TYPE,SIZE,IMNOX,STAT) !Midas binary format ENDIF ENDIF C 4000 FELEI = (FLINI-1)*NPIXA(1) + 1 FELEO = (FLINI-1)*NPIXC(1) + 1 CALL STFGET(IMNOA,(FELEI+PLINCI),NK,IAV,MADRID(PNTRA),STAT) C 4040 CALL RESMPA(MADRID(PNTRA),NPIXA(1),STA(1),STEP(1),NPIXC(1), + MADRID(PNTRW),FLINI,NLINES,NKK, + MADRID(PNTRX),MADRID(PNTRY),MADRID(PNTRZ)) C IF (PROCOP.EQ.'R') THEN !we work only on the rows IMNOX = IMNOC !write already to result frame CALL FLUFAK(IFLU,MADRID(PNTRW),NL,RFLU,CUTVAL) IF (IALL.EQ.1) + CALL STFPUT(IMNOX,FELEO,NL,MADRID(PNTRW),STAT) ENDIF C IF (IALL.EQ.0) THEN !we use intermediate disk file CALL STFPUT(IMNOX,FELEO,NL,MADRID(PNTRW),STAT) !write to disk FELEO = FELEO + NL FLINI = FLINI + NLINES C IF (FLINI .LE. NPIXA(2)) THEN FELEI = FELEI + NK CALL STFGET(IMNOA,(FELEI+PLINCI),NK,IAV, + MADRID(PNTRA),STAT) N = NPIXA(2) - FLINI + 1 IF (NLINES .GT. N) THEN NLINES = N NL = NLINES * NPIXC(1) NK = NLINES * NPIXA(1) ENDIF GOTO 4040 ENDIF ENDIF C IF (NAXIS.LT.3) + CALL STFCLO(IMNOI,STAT) !we can get rid of that buffer IF (PROCOP.EQ.'R') THEN !we work only on the rows IF (IFLU.EQ.1) THEN CUTS(3) = CUTVAL(1) * RFLU CUTS(4) = CUTVAL(2) * RFLU ELSE CUTS(3) = CUTVAL(1) CUTS(4) = CUTVAL(2) ENDIF GOTO 6000 !job already done ENDIF C C prepare processing along columns OPTIO = 0 NLINES = NSAVE !we must reset no. of new rows NL = NLINES * NPIXC(1) FLINI = INSTA !FLINI is updated in RESMPB FELEI = (FLINI-1)*NPIXC(1) + 1 FLINO = 1 !FLINO is updated here FELEO = 1 MAXO = SIZEC MXPROC = NPIXC(2) IF (IALL.EQ.0) THEN CALL STFGET(IMNOX,FELEI,NL,IAV,MADRID(PNTRW),STAT) ELSE RESPTR = PSAVE - (INSTA - 1)*NPIXC(1) !only use buffer offset PNTRW = RESPTR + FELEI - 1 ENDIF C 4400 IF (SWI .EQ. 0) THEN CALL RESMPB(OPTIO,MADRID(PNTRW),NPIXA,STA,STEP,NPIXC, + MADRID(PNTRC),FLINI,NLINES,MAXO,NPROC) ELSE CALL RESMPC(OPTIO,XREST,N1X,MADRID(PNTRW),NPIXA,STA, + STEP,NPIXC,ZFRAC,NZ,MADRID(PNTRC),FLINI, + NLINES,MXPROC,NPROC) ENDIF IF (NPROC.GT.0) THEN NK = NPROC * NPIXC(1) IF (IFLU.EQ.1) + CALL FLUFAC(MADRID(PNTRC),NK,RFLU) CALL STFPUT(IMNOC,(FELEO+PLINCO),NK,MADRID(PNTRC),STAT) C FELEO = FELEO + NK FLINO = FLINO + NPROC MAXO = MAXO - NK C IF (FLINO .LE. NPIXC(2)) THEN IF (OPTIO.NE.0) THEN FELEI = (FLINI-1)*NPIXC(1) + 1 IF (IALL.EQ.0) THEN CALL STFGET(IMNOX,FELEI,NL,IAV,MADRID(PNTRW),STAT) ELSE PNTRW = RESPTR + FELEI - 1 ENDIF ENDIF IF (STEP(2).GE.1.D0) THEN N = NPIXA(2) - FLINI + 1 !max no. of lines available IF (NLINES.GT.N) NLINES = N ENDIF GOTO 4400 ENDIF ENDIF C IF (NAXIS.EQ.3) THEN !if we work on planes of cube PLINCI = PLINCI + PLNSZI !move to next plane + loop PLINCO = PLINCO + PLNSZO IF (PLINCI.LT.(PLNSZI*NPIXA(3))) THEN FLINI = INSTA NK = NSAVE * NPIXA(1) NL = NSAVE * NPIXC(1) GOTO 4000 ENDIF ENDIF C IF (IFLU.EQ.1) THEN CUTS(3) = CUTVAL(1) * RFLU CUTS(4) = CUTVAL(2) * RFLU ELSE CUTS(3) = CUTVAL(1) CUTS(4) = CUTVAL(2) ENDIF ENDIF C C copy descriptors and update LHCUTS + HISTORY 6000 CALL DSCUPT(IMNOA,IMNOC,' ',STAT) CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNI,STAT) RNULL = NULCNT !convert NULL count to real CALL STKWRR('NULL',RNULL,1,1,UNI,STAT) CALL STSEPI C 10000 FORMAT('invalid stepsize ',F12.6,' in ',I1,'. dimension...') 10001 FORMAT('invalid NPIX(',I1,') of result frame (=',I5,') ...') C END SUBROUTINE FLUFAK(ISW,C,SIZE,FACT,RCUTS) C IMPLICIT NONE C REAL C(*),FACT,RCUTS(*) C INTEGER ISW,SIZE INTEGER N C IF (ISW.EQ.1) THEN DO 1000, N=1,SIZE !adjust result frame C(N) = C(N) * FACT IF (C(N).LT.RCUTS(1)) THEN !and check min,max RCUTS(1) = C(N) ELSE IF (C(N).GT.RCUTS(2)) THEN RCUTS(2) = C(N) ENDIF 1000 CONTINUE C ELSE DO 2000, N=1,SIZE !just check min,max IF (C(N).LT.RCUTS(1)) THEN RCUTS(1) = C(N) ELSE IF (C(N).GT.RCUTS(2)) THEN RCUTS(2) = C(N) ENDIF 2000 CONTINUE ENDIF C RETURN END SUBROUTINE FLUFAC(C,SIZE,FACT) C IMPLICIT NONE C REAL C(*),FACT C INTEGER SIZE INTEGER N C C adjust result frame DO 1000 N=1,SIZE C(N) = C(N) * FACT 1000 CONTINUE C RETURN END SUBROUTINE BUILDA(NOPI,STA,STEP,NOPC,XS,LS,HS,NKK) C DOUBLE PRECISION STA,STEP,XS(*) DOUBLE PRECISION X,XH,STHALF C INTEGER NOPI,NOPC,LS(*),HS(*),NKK INTEGER N,NX C NKK = 0 C C here for stepsize in x < 1.0 C IF (STEP.LT.1.0D0) THEN X = STA !set to center of 1.interval C 100 IF (X.LE.0.D0) THEN NKK = NKK + 1 X = X + STEP !move on GOTO 100 ENDIF C DO 200, N=NKK+1,NOPC NX = INT(X) + 1 XH = NX - 0.5 IF (X .LT. XH) THEN !X in [0,0.5] XS(N) = 1.D0 - (XH - X) NX = NX - 1 ELSE XS(N) = X - XH ENDIF LS(N) = NX IF (NX.LE.0) LS(N) = 1 IF (NX.GE.NOPI) THEN HS(N) = NOPI ELSE HS(N) = NX + 1 ENDIF X = X + STEP 200 CONTINUE C CO CO the above calculates the offsets for the original algorithm CO to interpolate between X and center before + after: CO CO NX = INT(X) + 1 CO XH = NX - 0.5 CO IF (X .LT. XH) THEN !X in [0,0.5] CO XFRAC = 1.D0 - (XH - X) CO NX = NX - 1 CO IF (NX.LE.0) THEN CO R1 = A(1+IOFFA) CO ELSE CO R1 = A(NX+IOFFA) CO ENDIF CO ELSE CO XFRAC = X - XH CO R1 = A(NX+IOFFA) CO ENDIF CO CO IF (NX.GE.NOPI) THEN CO R2 = A(NOPI+IOFFA) CO ELSE CO R2 = A(NX+1+IOFFA) CO ENDIF CO T(IOFFB) = R1 + (XFRAC*(R2-R1)) CO CO above formula is same as: A(NX)*(1.-XFRAC) + A(N1X)*XFRAC CO CO IOFFB = IOFFB + 1 CO X = X + STEP(1) CO C ELSE STHALF = STEP * 0.5D0 X = STA + STHALF !end of start sample 1100 IF (X.LE.0.D0) THEN NKK = NKK + 1 X = X + STEP !move on GOTO 1100 ENDIF C DO 1150, N=NKK+1,NOPC LS(N) = INT(X) XS(N) = X - LS(N) X = X + STEP 1150 CONTINUE ENDIF C RETURN END SUBROUTINE RESMPS(FLAG,IMA,IMC,A,C,NOPI,NOPC,D,E) C REAL A(*),C(*),D(*),E(*) C INTEGER FLAG,IMA,IMC,NOPI(*),NOPC(*) INTEGER N,NA,NY,NPXC,FELEI,FELEO,IAV,STAT C FELEI = 1 FELEO = 1 C C omit pixels + rows C IF (FLAG.EQ.2) THEN NA = 1 CALL STFGET(IMA,FELEI,NOPI(1),IAV,A,STAT) !first row DO 80, N=1,NOPC(1) C(N) = A(NA) NA = NA + 2 80 CONTINUE CALL STFPUT(IMC,FELEO,NOPC(1),C,STAT) C IF (NOPC(2) .GT. 1) THEN NPXC = NOPC(1) FELEO = FELEO + NPXC C DO 500, NY=2,NOPC(2) FELEI = FELEI + 2*NOPI(1) !omit even rows CALL STFGET(IMA,FELEI,NOPI(1),IAV,A,STAT) NA = 1 DO 330, N=1,NOPC(1) C(N) = A(NA) NA = NA + 2 330 CONTINUE CALL STFPUT(IMC,FELEO,NOPC(1),C,STAT) FELEO = FELEO + NPXC 500 CONTINUE ENDIF C C rebin between pixels + rows C ELSE CALL STFGET(IMA,FELEI,NOPI(1),IAV,A,STAT) !first row C(1) = A(1) NA = 2 DO 800, N=2,NOPC(1),2 C(N) = (A(NA-1)+A(NA))*0.5 C(N+1) = A(NA) NA = NA + 1 800 CONTINUE C(NOPC(1)) = A(NOPI(1)) CALL STFPUT(IMC,FELEO,NOPC(1),C,STAT) C IF (NOPC(2) .GT. 1) THEN NPXC = NOPC(1) FELEO = FELEO + NPXC C DO 10000, NY=2,NOPC(2),2 FELEI = FELEI + NOPI(1) CALL STFGET(IMA,FELEI,NOPI(1),IAV,A,STAT) E(1) = A(1) NA = 2 DO 1800, N=2,NPXC,2 E(N) = (A(NA-1)+A(NA))*0.5 E(N+1) = A(NA) NA = NA + 1 1800 CONTINUE E(NPXC) = A(NOPI(1)) C DO 2000, N=1,NPXC D(N) = (C(N)+E(N))*0.5 C(N) = E(N) 2000 CONTINUE C CALL STFPUT(IMC,FELEO,NOPC(1),D,STAT) FELEO = FELEO + NPXC CALL STFPUT(IMC,FELEO,NOPC(1),E,STAT) FELEO = FELEO + NPXC 10000 CONTINUE ENDIF ENDIF C RETURN END SUBROUTINE RESMPA(A,NOPI,STA,STEP,NOPC,T,STLNI,LINCNT, + NKK,XS,LS,HS) C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine RESAMPA version 1.50 840207 C K. Banse ESO - Garching C 1.80 870129 C C.KEYWORDS C bulk data frame, rebinning C C.PURPOSE C get data values at given real pixel no's. C C.ALGORITHM C work first along the lines, change lines to columns C and work along the lines of the transposed image C C.INPUT/OUTPUT C call as RESAMPA(A,NPIXA,STA,STEP,NPIXB,T,STLNI,LINCNT,NKK,XS,LS,HS) C C input par: C A: R*4 array input array of data C NOPI: I*4 no. of pixels in x of input array C STA: R*8 normed x-start for output array C STEP: R*8 normed x-stepsize for output array C NOPC: I*4 no. of pixels in x for output array C STLNI: I*4 start line of input frame C LINCNT: I*4 line count for buffer of line size NPIXB(1) C NKK: I*4 no. of X values .LE. 0 C XS: R*8 array XFRAC for each X C LS: I*4 array low index for each X C HS: I*4 array high index for each X C C output par: C T: R*4 array output array C C-------------------------------------------------- C IMPLICIT NONE C REAL A(*),T(*) REAL PIX,R1,R2,VALU,CUTVAL(2) REAL ZNULL C DOUBLE PRECISION STA,STEP,XS(*) DOUBLE PRECISION XFRAC,XREST DOUBLE PRECISION Z,ZFRAC,STHALF C INTEGER NOPI,NOPC,STLNI,LINCNT,NKK,LS(*),HS(*) INTEGER IOFFA,IOFFB,L,N,NK INTEGER NN,N1X,NX,NZ INTEGER NULCNT C COMMON /NULLC/ NULCNT,ZNULL,CUTVAL C IOFFA = 0 IOFFB = 1 C C ---------------------------- C here for stepsize in x < 1.0 C ---------------------------- C IF (STEP.LT.1.0D0) THEN NK = NKK + 1 C DO 800, L=1,LINCNT !loop over all input lines IF (NKK .GT. 0) THEN DO 550, N=1,NKK T(IOFFB) = ZNULL IOFFB = IOFFB + 1 550 CONTINUE ENDIF C DO 600, N=NK,NOPC R1 = A(LS(N)+IOFFA) R2 = A(HS(N)+IOFFA) T(IOFFB) = R1 + (XS(N)*(R2-R1)) IOFFB = IOFFB + 1 600 CONTINUE C IOFFA = IOFFA + NOPI 800 CONTINUE C C ----------------------------- C here for stepsize in x >= 1.0 C ----------------------------- C ELSE C STHALF = STEP * 0.5D0 Z = STA - STHALF !go half a new stepsize back IF (Z.LE.0.D0) THEN Z = 0.D0 NZ = 0 ELSE NZ = INT(Z) ENDIF ZFRAC = Z - NZ !from there go NK = NKK + 2 C DO 1800, L=1,LINCNT !loop over all input lines IF (NKK .GT. 0) THEN DO 1200, N=1,NKK T(IOFFB) = ZNULL IOFFB = IOFFB + 1 1200 CONTINUE ENDIF C C work on first new pixel differently... C NX = LS(NK-1) XFRAC = XS(NK-1) PIX = -(ZFRAC*A(NZ+1+IOFFA)) DO 1300, NN=1,NX-NZ !sum over full regions PIX = PIX + A(NN+NZ+IOFFA) 1300 CONTINUE N1X = NX + 1 !point to next region T(IOFFB) = PIX + XFRAC*A(N1X+IOFFA) !and add fraction of that IOFFB = IOFFB + 1 XREST = 1.D0 - XFRAC !update XREST for next iteration C C now the remaining pixels... DO 1500, N=NK,NOPC-1 NX = LS(N) XFRAC = XS(N) PIX = XREST * A(N1X+IOFFA) !init to remainder from last time N1X = N1X + 1 !point to 1. full region IF (NX .GE. N1X) THEN DO 1400, NN=N1X,NX !sum over full regions PIX = PIX + A(NN+IOFFA) 1400 CONTINUE ENDIF N1X = NX + 1 !point to next region and T(IOFFB) = PIX + (XFRAC * A(N1X+IOFFA)) !add fraction of that IOFFB = IOFFB + 1 XREST = 1.D0 - XFRAC !update XREST for next iteration 1500 CONTINUE C C take care of end of row ... NX = LS(NOPC) XFRAC = XS(NOPC) IF (N1X.LE.NOPI) THEN !init to remainder from last time PIX = XREST * A(N1X+IOFFA) ELSE PIX = XREST * A(NOPI+IOFFA) ENDIF N1X = N1X + 1 !point to 1. full region C DO 1600, NN=N1X,NX !sum over full regions IF (NN.LE.NOPI) THEN PIX = PIX + A(NN+IOFFA) ELSE PIX = PIX + A(NOPI+IOFFA) ENDIF 1600 CONTINUE C C PN 3/99 C execute this part only if no. of pixels in outima greater/equal 2 IF( NOPC .GE. 2) THEN N1X = NX + 1 !point to next region IF (N1X.LE.NOPI) THEN VALU = A(N1X+IOFFA) ELSE VALU = A(NOPI+IOFFA) ENDIF T(IOFFB) = PIX + (XFRAC*VALU) !and add fraction of that IOFFB = IOFFB + 1 C C update pointers to offset in data arrays IOFFA = IOFFA + NOPI ENDIF 1800 CONTINUE ENDIF C RETURN END SUBROUTINE RESMPB(FLAG,TA,NPIXA,STA,STEP,NPIXB,B,STLNI,LINCNT, + MXIOB,NPROC) C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine RESMPB version 1.00 970425 C K. Banse ESO - Garching C C.KEYWORDS C bulk data frame, rebinning C C.PURPOSE C get data values at given real pixel no's. C C.ALGORITHM C work along the columns of the intermediate frame C on all the rows available in the buffer C C.INPUT/OUTPUT C call as RESMPB(FLAG,TA,NPIXA,STA,STEP,NPIXB,B,STLNI,LINCNT, C MXIOB,STLNO,NPROC) C C input par: C TA: R*4 array input buffer of NPIXB(1)*LINCNT size C NPIXA: I*4 array no. of pixels in x+y of input frame C STEP: R*8 array normed stepsize for output array C NPIXB: I*4 no. of pixels in x+y for output frame C LINCNT: I*4 no. of lines in input buffer C MXIOB: I*4 max. no. of values to be produced C C output par: C B: R*4 array output buffer of data C NPROC: I*4 no. of output lines produced C C in/output par: C FLAG: I*4 = 0, at start + set to 1, if zero data C finished C STA: R*8 array normed start for output buffer C STLNI: I*4 starting line of input buffer C C-------------------------------------------------- C IMPLICIT NONE C REAL TA(*),B(*) REAL R1,VALU REAL ZNULL,CUTVAL(2) C DOUBLE PRECISION STA(*),STEP(*) DOUBLE PRECISION XFRAC,X,XH C INTEGER FLAG,NPIXA(*),NPIXB(*),STLNI,LINCNT,MXIOB,NPROC INTEGER LSTLN1,L,KOFF1,KOFF2 INTEGER NX,IOB,NOPI INTEGER NULCNT,KX C COMMON /NULLC/ NULCNT,ZNULL,CUTVAL C C use TA as input - TA is of size NPIXB(1) * LINCNT C NOPI = NPIXB(1) NPROC = 0 !no. of processed lines IOB = 1 !index for output buffer B LSTLN1 = NPIXA(2) - STLNI C C test FLAG C IF (FLAG .EQ. 0) THEN !first time, so check if outside IF (STA(2) .LT. 0.0D0) THEN C C determine no. of NULL lines X = STA(2) 2800 IF (X.LT.0.0D0) THEN NPROC = NPROC + 1 X = X + STEP(2) IF (NPROC.LT.LINCNT) GOTO 2800 ENDIF C STA(2) = X DO 3300, IOB=1,NPROC*NOPI B(IOB) = ZNULL 3300 CONTINUE C NULCNT = NULCNT + (NPROC*NOPI) IF (CUTVAL(1).GT.ZNULL) THEN CUTVAL(1) = ZNULL ELSE IF (CUTVAL(2).LT.ZNULL) THEN CUTVAL(2) = ZNULL ENDIF RETURN ENDIF C FLAG = 1 ENDIF X = STA(2) C C here we always have a line `before' C 6600 NX = INT(X) + 1 XH = NX - 0.5 KX = NX - STLNI + 1 IF (X .LT. XH) THEN !X in [0,0.5[ XFRAC = 1.D0 - (XH - X) NX = NX - 1 !NX = NX - 1 KX = KX - 1 IF (KX.GE.LINCNT) THEN NX = NX - 1 !for security one more line at bottom GOTO 6900 ENDIF IF (NX.LT.STLNI) THEN KOFF1 = 1 ELSE KOFF1 = (KX-1)*NOPI + 1 ENDIF ELSE !X in [0.5,1.0] IF (KX.GE.LINCNT) THEN NX = NX - 1 !for security one more line at bottom GOTO 6900 ENDIF XFRAC = X - XH KOFF1 = (KX-1)*NOPI + 1 ENDIF C IF (KX.GE.LSTLN1) THEN KOFF2 = (LSTLN1 * NOPI) + 1 ELSE KOFF2 = (KX * NOPI) + 1 ENDIF C DO 6800, L=IOB,IOB+NOPI-1 !loop over all columns R1 = TA(KOFF1) VALU = R1 + (XFRAC*(TA(KOFF2)-R1)) B(L) = VALU C C above formula is same as: A(NX)*(1.-XFRAC) + A(N1X)*XFRAC C IF (CUTVAL(1).GT.VALU) THEN CUTVAL(1) = VALU ELSE IF (CUTVAL(2).LT.VALU) THEN CUTVAL(2) = VALU ENDIF C KOFF1 = KOFF1 + 1 KOFF2 = KOFF2 + 1 6800 CONTINUE C IOB = IOB + NOPI NPROC = NPROC + 1 IF (IOB .GT. MXIOB) RETURN !end_of_job check C X = X + STEP(2) IF (NPROC.LT.LINCNT) GOTO 6600 C 6900 STA(2) = X STLNI = NX - 1 IF (STLNI.LT.1) STLNI = 1 C RETURN END SUBROUTINE RESMPC(FLAG,XREST,N1X,TA,NPIXA,STA,STEP,NPIXB,ZFRAC,NZ, + B,STLNI,LINCNT,MXPROC,NPROC) C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine RESMPC version 1.00 970425 C K. Banse ESO - Garching C C.KEYWORDS C bulk data frame, rebinning C C.PURPOSE C get data values at given real pixel no's. C C.ALGORITHM C work along the lines of the transposed image C C.INPUT/OUTPUT C call as RESMPC(FLAG,TA,NPIXA,STA,STEP,NPIXB,B,STLNI,LINCNT, C MXPROC,STLNO,NPROC) C C input par: C TA: R*4 array input buffer of NPIXB(1)*NPIXA(2) temp frame C NPIXA: I*4 array no. of pixels in x+y of input frame C STEP: R*8 array normed stepsize for output array C NPIXB: I*4 no. of pixels in x+y for output frame C MXPROC: I*4 max. no. of rowd to be produced C C output par: C B: R*4 array output buffer of data C NPROC: I*4 no. of output lines produced C C in/output par: C FLAG: I*4 = 0, at start + set to 1, if zero data C finished C STA: R*8 array normed start for output buffer C STLNI: I*4 starting line of input buffer C C-------------------------------------------------- C IMPLICIT NONE C REAL TA(*),B(*) REAL PIXVAL,VALU REAL ZNULL,CUTVAL(2) C INTEGER FLAG,NPIXA(*),NPIXB(*),STLNI,LINCNT,MXPROC,NPROC,NZ INTEGER LASTIN,L,NX,IOB,KL,ML,JL INTEGER MZ,NSAV,NN,N1X,NOLI,NOPI,STLNIM INTEGER NULCNT,N2X C DOUBLE PRECISION STA(*),STEP(*) DOUBLE PRECISION XFRAC,X,XREST DOUBLE PRECISION ZFRAC,STHALF C COMMON /NULLC/ NULCNT,ZNULL,CUTVAL C C use TA as input - TA is of size NPIXB(1) * chunk_lines C NOPI = NPIXB(1) NOLI = NPIXA(2) NPROC = 0 STLNIM = STLNI - 1 LASTIN = LINCNT + STLNIM IOB = 1 C C here for stepsize >= 1.0 C STHALF = STEP(2) * 0.5D0 X = STA(2) + STHALF !end of start sample C IF (FLAG .EQ. 0) THEN IF (X .GT. 0.0D0) THEN !determine no. of NULL lines FLAG = 1 ELSE 2800 IF (X .LE. 0.0D0) THEN NPROC = NPROC + 1 X = X + STEP(2) IF (NPROC.LT.LINCNT) GOTO 2800 ENDIF C STA(2) = X - STHALF DO 3300, IOB=1,NPROC*NOPI B(IOB) = ZNULL 3300 CONTINUE C NULCNT = NULCNT + (NPROC*NOPI) MXPROC = MXPROC - NPROC IF (CUTVAL(1).GT.ZNULL) THEN CUTVAL(1) = ZNULL ELSE IF (CUTVAL(2).LT.ZNULL) THEN CUTVAL(2) = ZNULL ENDIF RETURN ENDIF ENDIF C C if first real line, build up NX, N1X, KL, XFRAC, XREST C IF (FLAG .EQ. 1) THEN !work on first new pixel differently... IF (NZ .GT. 0) THEN MZ = NZ - STLNIM ELSE MZ = 0 ENDIF X = STA(2) + STHALF NX = INT(X) N1X = NX + 1 !point to next region XFRAC = X - NX XREST = 1.D0 - XFRAC !update XREST for next iteration KL = MZ*NOPI + 1 JL = (NX - STLNIM)*NOPI + 1 NSAV = MZ * NOPI !save for loop C IF (NX.GT.NZ) THEN DO 5100, L=1,NOPI !loop over all columns PIXVAL = -(ZFRAC*TA(KL)) ML = NSAV DO 5000, NN=1,NX-NZ !sum over full regions PIXVAL = PIXVAL + TA(ML+L) ML = ML + NOPI 5000 CONTINUE VALU = PIXVAL + (XFRAC * TA(JL)) !add fraction of that B(IOB) = VALU IOB = IOB + 1 IF (CUTVAL(1).GT.VALU) THEN CUTVAL(1) = VALU ELSE IF (CUTVAL(2).LT.VALU) THEN CUTVAL(2) = VALU ENDIF KL = KL + 1 JL = JL + 1 5100 CONTINUE C ELSE DO 5200, L=IOB,IOB+NOPI-1 !loop over all columns VALU = -(ZFRAC*TA(KL)) + (XFRAC*TA(JL)) B(L) = VALU IF (CUTVAL(1).GT.VALU) THEN CUTVAL(1) = VALU ELSE IF (CUTVAL(2).LT.VALU) THEN CUTVAL(2) = VALU ENDIF KL = KL + 1 JL = JL + 1 5200 CONTINUE IOB = IOB + NOPI ENDIF C FLAG = 2 NPROC = 1 MXPROC = MXPROC - 1 ENDIF C C now the remaining rows... C IF (FLAG.EQ.2) THEN 6000 X = X + STEP(2) NX = INT(X) C IF (NX .GE. LASTIN) THEN !we have to check NX... IF((NPROC.EQ.0).OR.(LASTIN.EQ.NOLI)) GOTO 9090 X = X - STEP(2) STLNI = INT(X) N1X = STLNI + 1 !point to next region XFRAC = X - STLNI XREST = 1.D0 - XFRAC !update XREST for next iteration STA(2) = X - STHALF RETURN ENDIF C XFRAC = X - NX N2X = N1X + 1 !point to 1. full region NSAV = (N2X - STLNI) * NOPI KL = (N1X - STLNI) * NOPI + 1 JL = (NX - STLNIM) * NOPI + 1 C IF (NX .GE. N2X) THEN DO 7500, L=1,NOPI PIXVAL = XREST * TA(KL) !init to remainder from last time ML = NSAV DO 7400, NN=N2X,NX !sum over full regions PIXVAL = PIXVAL + TA(ML+L) ML = ML + NOPI 7400 CONTINUE VALU = PIXVAL + (XFRAC * TA(JL)) !add fraction of that B(IOB) = VALU IOB = IOB + 1 IF (CUTVAL(1).GT.VALU) THEN CUTVAL(1) = VALU ELSE IF (CUTVAL(2).LT.VALU) THEN CUTVAL(2) = VALU ENDIF KL = KL + 1 !instead of TA(KL+L) JL = JL + 1 7500 CONTINUE C ELSE DO 7700, L=IOB,IOB+NOPI-1 VALU = XREST*TA(KL) + XFRAC*TA(JL) B(L) = VALU IF (CUTVAL(1).GT.VALU) THEN CUTVAL(1) = VALU ELSE IF (CUTVAL(2).LT.VALU) THEN CUTVAL(2) = VALU ENDIF KL = KL + 1 !instead of TA(KL+L) JL = JL + 1 7700 CONTINUE IOB = IOB + NOPI ENDIF C N1X = NX + 1 XREST = 1.D0 - XFRAC !update XREST for next iteration NPROC = NPROC + 1 MXPROC = MXPROC - 1 C C check, if we reached last row IF (MXPROC.EQ.1) THEN FLAG = 3 IF (NPROC.LT.LINCNT) THEN GOTO 9000 !move to final section right away ELSE GOTO 8000 ENDIF ENDIF IF (NPROC.LT.LINCNT) GOTO 6000 C 8000 STLNI = NX STA(2) = X - STHALF RETURN ENDIF C C take care of last row (FLAG = 3) C 9000 X = X + STEP(2) NX = INT(X) C 9090 IF (NX.GE.LASTIN) THEN !we have to check NX... IF ((NX.EQ.LASTIN) .AND. (LASTIN.EQ.NOLI)) GOTO 9099 X = X - STEP(2) STLNI = INT(X) N1X = STLNI + 1 !point to next region XFRAC = X - STLNI XREST = 1.D0 - XFRAC !update XREST for next iteration STA(2) = X - STHALF RETURN ENDIF C 9099 XFRAC = X - NX N2X = N1X + 1 !point to 1. full region IF (N1X.LE.NOLI) THEN !init to remainder from last time ML = (N1X - STLNI)*NOPI + 1 ELSE ML = (NOLI - STLNI)*NOPI + 1 ENDIF N1X = NX + 1 !point to next region IF (N1X.LE.NOLI) THEN KL = (N1X - STLNI)*NOPI + 1 ELSE KL = (NOLI - STLNI)*NOPI + 1 ENDIF C IF (NX .GE. N2X) THEN DO 9700, L=1,NOPI !loop over all columns PIXVAL = XREST * TA(ML) DO 9600, NN=N2X,NX !sum over full regions IF (NN.LE.NOLI) THEN JL = NN - STLNI ELSE JL = NOLI - STLNI ENDIF PIXVAL = PIXVAL + TA(JL*NOPI+L) 9600 CONTINUE VALU = PIXVAL + XFRAC*TA(KL) !and add fraction of that B(IOB) = VALU IOB = IOB + 1 IF (CUTVAL(1).GT.VALU) THEN CUTVAL(1) = VALU ELSE IF (CUTVAL(2).LT.VALU) THEN CUTVAL(2) = VALU ENDIF KL = KL + 1 !instead of TA(KL+L) ML = ML + 1 9700 CONTINUE C ELSE DO 9770, L=IOB,IOB+NOPI-1 !loop over all columns VALU = XREST*TA(ML) + XFRAC*TA(KL) B(L) = VALU IF (CUTVAL(1).GT.VALU) THEN CUTVAL(1) = VALU ELSE IF (CUTVAL(2).LT.VALU) THEN CUTVAL(2) = VALU ENDIF KL = KL + 1 !instead of TA(KL+L) ML = ML + 1 9770 CONTINUE IOB = IOB + NOPI ENDIF MXPROC = 0 NPROC = NPROC + 1 C RETURN END SUBROUTINE RESMPX(A,OFFA,NDIMA,STAX,STEPX,DERIV, + B,NDIMB,OFFB) C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine RESMPX version 1.00 900713 C K. Banse ESO - Garching C 1.10 920109 C C.KEYWORDS C bulk data frame, rebinning C C.PURPOSE C get data values at given real pixel no's. of a single row from an image C C.ALGORITHM C C.INPUT/OUTPUT C call as RESMPX(A,OFFA,NDIMA,STAX,STEPX,DERIV,B,NDIMB,OFFB,ZNULL) C C input par: C A: R*4 array row of input data frame C OFFA: I*4 offset within total input array C NDIMA: I*4 size of row of A C STAX: R*8 start in x of output array C STEPX: R*8 normalized stepsize in x of output array C DERIV: R*4 array second derivative of row C NDIMB: I*4 size of output row C OFFB: I*4 offset for output array C C output par: C B: R*4 array output array of data C C.VERSIONS C 1.00 creation C 1.10 clean + speed up C C-------------------------------------------------- C IMPLICIT NONE C REAL A(1),DERIV(1),B(1) REAL ZNULL,CUTVAL(2) C DOUBLE PRECISION STAX,STEPX DOUBLE PRECISION X,GX,RF,QQ,RR C INTEGER NDIMA,NDIMB,OFFA,OFFB INTEGER N,NO1,NO2,KOFF,NULCNT C COMMON /NULLC/ NULCNT,ZNULL,CUTVAL C GX = NDIMA X = STAX + (0.5D0*STEPX) RF = 1.0D0 / 6.0D0 KOFF = OFFB + 1 C DO 1800 N=1,NDIMB IF ((X.LT.1.0D0) .OR. (X.GT.GX)) THEN B(KOFF) = ZNULL NULCNT = NULCNT + 1 ELSE C NO1 = INT(X) IF (NO1.GE.NDIMA) THEN B(KOFF) = A(NO1+OFFA) ELSE NO2 = NO1 + 1 !get next pixel in x RR = NO2 - X QQ = 1.D0 - RR C C original formula: PIX = RR*A(NO1+OFFA) + QQ*A(NO2+OFFA) + C ((RR**3-RR)*DERIV(NO1) + (QQ**3-QQ)*DERIV(NO2)) * RF C B(KOFF) = RR*A(NO1+OFFA) + QQ*A(NO2+OFFA) + + ( (RR*(RR*RR-1)*DERIV(NO1)) + + (QQ*(QQ*QQ-1.)*DERIV(NO2)) ) * RF ENDIF ENDIF C KOFF = KOFF + 1 X = X + STEPX 1800 CONTINUE C RETURN END SUBROUTINE DERIV2(A,IOFF,DERV2,TEMP,NDIM) C IMPLICIT NONE C INTEGER NDIM,IOFF,K,I,IM1 C REAL A(NDIM),DERV2(NDIM),TEMP(NDIM) REAL P C C init DERV2(1) = 0. !use "natural" spline TEMP(1) = 0. C C compute array DERV2 (second derivative of A) C C original loop: C DO 2000 I=2,NDIM-1 C K = IOFF + I C P = 0.5*DERV2(I-1) + 2.0 !interval size = 1.0 C DERV2(I) = -0.5 / P C TEMP(I) = ( 3.0*(A(K+1) - A(K) - A(K) + A(K-1)) C + - 0.5*TEMP(I-1) ) / P C2000 CONTINUE C DO 2000 I=2,NDIM-1 IM1 = I - 1 K = IOFF + I P = 1.0 / (DERV2(IM1) + 4.0) !this is 0.5*(1./P) of above DERV2(I) = - P TEMP(I) = ( 6.0*(A(K+1)-A(K)-A(K)+A(K-1)) - TEMP(IM1) ) * P 2000 CONTINUE C DERV2(NDIM) = 0.0 C DO 3000, I=NDIM-1,1,-1 DERV2(I) = DERV2(I)*DERV2(I+1) + TEMP(I) 3000 CONTINUE C RETURN END