C @(#)fripple.for 17.1.1.1 (ES0-DMD) 01/25/02 17:56:13 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 PROGRAM FRIPLE C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 11:28 - 7 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: J.D.Ponz, but I did not code it! C C.IDENTIFICATION C program FRIPPLE version 1.0 851030 C correct some bugs 900619 MP C C.Keywords C filtering C C.PURPOSE C correct 1-D images for periodic modulation (FRIPLE) C C.ALGORITHM C Fold flux with period, calculate mean flux per phase bin, determine C deviation from total mean, adjust to mean level. C C.INPUT/OUTPUT C the following keywords are used: C IN_A/C/1/8 name of input image C INPUTR/R/1/1 period C INPUTI/I/1/2 first and last pixel to be considered C OUT_A/C/1/8 name of ripple corrected image C C -------------------------------------------------------------------- C IMPLICIT NONE C INTEGER MADRID CHARACTER*60 INIMA,OUTIMA CHARACTER*80 OUTEXT,HISTOR CHARACTER IDENT*72,CUNIT*64 INTEGER NAXIS,NPIX,ISTAT,IAV INTEGER*8 PNTRIN,PNTRO INTEGER RANGE(2),ERROR(2),NWORK INTEGER NBIN,I,SCNPIX INTEGER*8 SCRPTR INTEGER KUN,KNUL,IMIN,IMOUT REAL PERIOD,AMP(500),TEST C DOUBLE PRECISION DSTART(3),DSTEP(3) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C ... set up MIDAS environment C CALL STSPRO('FRIPLE') C C ... get name of input frame, map it C CALL STKRDC('IN_A',1,1,60,IAV,INIMA,KUN,KNUL,ISTAT) IF (ISTAT.NE.0) GO TO 20 CALL STIGET(INIMA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,1,NAXIS,NPIX, + DSTART,DSTEP,IDENT,CUNIT,PNTRIN,IMIN,ISTAT) IF (ISTAT.NE.0) GO TO 20 IF (NAXIS.NE.1) THEN CALL STTPUT('ERROR: input frame must be one-dimensional', + ISTAT) ERROR(1) = 1 GO TO 20 END IF C C ... get name of output frame, map it C CALL STKRDC('OUT_A',1,1,60,IAV,OUTIMA,KUN,KNUL,ISTAT) CALL STIPUT(OUTIMA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,NAXIS,NPIX, + DSTART,DSTEP,IDENT,CUNIT,PNTRO,IMOUT,ISTAT) C C ... get period, branch according to type (integer or real) of C ... period to subroutines IAMPLT or RAMPLT to find C ... amplitude of modulation C CALL STKRDR('INPUTR',1,1,IAV,PERIOD,KUN,KNUL,ISTAT) C CALL STKRDI('INPUTI',1,2,IAV,RANGE,KUN,KNUL,ISTAT) C TEST = PERIOD - INT(PERIOD) C ! period near-integer? IF (TEST*NPIX/PERIOD.LT.0.05) THEN NBIN = PERIOD CALL IAMPLT(MADRID(PNTRIN),PERIOD,RANGE,AMP) ! no, real ELSE ! use a fixed number of 20 phase bins NBIN = 20 C C ... in this case, map SCRAT file to hold intermediate results C ... first, determine number of pixels needed C SCNPIX = INT(20*NPIX/PERIOD) NWORK = SCNPIX*4 CALL TDMGET(NWORK,SCRPTR,ISTAT) C CALL SXIPUT('SCRAT','X',1,SCNPIX,START,STEP,IDENT,CUNIT, C + SCRPTR,ISTAT) CALL RAMPLT(MADRID(PNTRIN),SCNPIX,PERIOD,RANGE,AMP, + MADRID(SCRPTR)) END IF C C ... communicate results: C CALL STTPUT('Results:',ISTAT) DO 10 I = 1,NBIN,4 WRITE (OUTEXT,9000) I,I + 3,AMP(I),AMP(I+1),AMP(I+2), + AMP(I+3) CALL STTPUT(OUTEXT,ISTAT) 10 CONTINUE C C ... the ripple correction is, then, done by subroutines C ... IRIPC and RRIPC, respectively C IF (TEST*NPIX/PERIOD.LT.0.05) THEN CALL IRIPC(MADRID(PNTRIN),NPIX,PERIOD,RANGE,AMP, + MADRID(PNTRO)) ELSE CALL RRIPC(MADRID(PNTRIN),NPIX,PERIOD,RANGE,AMP, + MADRID(SCRPTR),SCNPIX,MADRID(PNTRO)) END IF C C ... copy descriptor LHCUTS C CALL STDCOP(IMIN,IMOUT,4,'LHCUTS',ISTAT) C C ... copy and update descriptor HISTOR C CALL STDCOP(IMIN,IMOUT,4,'HISTORY',ISTAT) CALL STKRDC('HISTORY',1,1,80,IAV,HISTOR,KUN,KNUL,ISTAT) CALL STDWRC(IMOUT,'HISTOR',1,HISTOR,-1,80,KUN,ISTAT) C C ... the end C CALL STSEPI STOP C C ... direct trouble makers here: C 20 ERROR(1) = ISTAT CALL STKWRI('PROGSTAT',ERROR,1,2,KUN,ISTAT) CALL STSEPI STOP C 9000 FORMAT ('Phase bins No.',I3,' to',I3,' have ','amplitudes ', + 3 (F7.4,1X),F7.4) END C C C SUBROUTINE IAMPLT(IN,PERIOD,RANGE,AMP) C C INTEGER NCYCL,I,J,RANGE(2),IPER, IR1, IR2 REAL PERIOD,IN(1),AMP(500),SUM(500),MEAN C ! consider only integer number of cycles NCYCL = INT((RANGE(2)-RANGE(1)+1)/PERIOD) C IPER = PERIOD DO 10 I = 1,IPER SUM(I) = 0 10 CONTINUE C C ... build up flux in phase bins C ... to compensate for slopes in input data, refer to mean C ... flux per cycle C IR1 = RANGE(1) IR2 = RANGE(1) + (NCYCL-1)*PERIOD DO 40 I = IR1, IR2,IPER MEAN = 0 DO 20 J = 1,IPER MEAN = MEAN + IN(I+J-1) 20 CONTINUE ! average flux in one cycle MEAN = MEAN/PERIOD DO 30 J = 1, IPER SUM(J) = SUM(J) + IN(I+J-1)/MEAN 30 CONTINUE 40 CONTINUE C C ... division into mean flux yields the relative amplitude of C ... each phase bin C MEAN = 0 DO 50 I = 1,IPER MEAN = MEAN + SUM(I) 50 CONTINUE MEAN = MEAN/PERIOD C DO 60 I = 1,IPER AMP(I) = MEAN/SUM(I) 60 CONTINUE C RETURN END C C C SUBROUTINE IRIPC(IN,NPIX,PERIOD,RANGE,AMP,OUT) C C INTEGER RANGE(2),NPIX,OFF,PHASE,I,J,INDEX REAL IN(1),OUT(1),PERIOD,AMP(500) C C ... take into account possible phase offset introduced if RANGE C ... was not defaulted to [1,NPIX] C OFF = NINT(((RANGE(1)-1.)/PERIOD-INT((RANGE(1)-1.)/PERIOD))* + PERIOD) C DO 20 I = 1,NPIX,PERIOD DO 10 J = 1,PERIOD INDEX = I + J - 1 IF (INDEX.GT.NPIX) GO TO 30 PHASE = OFF + J IF (PHASE.GT.PERIOD) PHASE = PHASE - PERIOD ! do the FRIPLE compensation OUT(INDEX) = IN(INDEX)*AMP(PHASE) 10 CONTINUE 20 CONTINUE C 30 RETURN END C C C SUBROUTINE RAMPLT(IN,SCNPIX,PERIOD,RANGE,AMP,SCRAT) C INTEGER RANGE(2),PIXLO,PIXHI,I,J,SCNPIX,BIN,START,OFF,END REAL PERIOD,IN(1),AMP(500),SUM(20),FRACLO,FRACHI,XLO,XHI,MEAN, + SCRAT(1),STEP C DO 10 I = 1,20 SUM(I) = 0. 10 CONTINUE C C ... rebin to stepsize corresponding to one phasebin C ... the algorithm is a pure projection, no interpolation is done C ... (since the origin and nature of the ripple are not known) C STEP = 0.05*PERIOD DO 30 I = 1,SCNPIX ! ith pixel extends from i to i+1 in both spaces XLO = (I-1)*STEP + 1 ! lower and upper limit of projected output pixel XHI = XLO + STEP PIXLO = INT(XLO) PIXHI = INT(XHI) ! all in same output pixel IF ((PIXHI-PIXLO).EQ.0) THEN FRACLO = XHI - XLO SCRAT(I) = FRACLO*IN(PIXLO) ! at least two output pixels hit ELSE FRACLO = PIXLO + 1 - XLO FRACHI = XHI - PIXHI SCRAT(I) = FRACLO*IN(PIXLO) + FRACHI*IN(PIXHI) ! if in fact more than two DO 20 J = PIXLO + 1,PIXHI - 1 SCRAT(I) = SCRAT(I) + IN(J) 20 CONTINUE END IF 30 CONTINUE C C ... evaluate the periodicity over the given range C ... do so over an integer number of periods C ... compensate for a slope in the input data by averaging over C ... one cycle of 20 pixels each C START = 20*INT((RANGE(1)-1.)/PERIOD) + 1 ! phase offset due to zeropoint used OFF = 20*INT((RANGE(1)-1.)/PERIOD) + 1 - START END = 20*INT((RANGE(2)-1.)/PERIOD) + 1 DO 60 I = START,END,20 MEAN = 0 DO 40 J = 1,20 MEAN = MEAN + SCRAT(I+J-1) 40 CONTINUE MEAN = MEAN/20 DO 50 J = 1,20 C BIN = JMOD(I+J-1+OFF,20) + 1 BIN = MOD(I+J-1+OFF,20) + 1 SUM(BIN) = SUM(BIN) + SCRAT(I+J-1)/MEAN 50 CONTINUE 60 CONTINUE C C ... division into mean value yields the relative amplitudes of C ... each each phas C MEAN = 0 DO 70 I = 1,20 MEAN = MEAN + SUM(I) 70 CONTINUE MEAN = MEAN/20. C DO 80 I = 1,20 AMP(I) = MEAN/SUM(I) 80 CONTINUE C RETURN END C C C SUBROUTINE RRIPC(IN,NPIX,PERIOD,RANGE,AMP,SCRAT,SCNPIX,OUT) C C INTEGER RANGE(2),NPIX,SCNPIX,I,J,BIN,PIXLO,PIXHI REAL IN(1),OUT(1),PERIOD,AMP(500),XLO,XHI,FRACLO,FRACHI, + SCRAT(1),STEP C DO 10 I = 1,NPIX ! initialize output frame to zero OUT(I) = 0 10 CONTINUE C C ... here the amplitude correction is done C DO 20 I = 1,SCNPIX C BIN = JMOD(I,20) + 1 BIN = MOD(I,20) + 1 SCRAT(I) = SCRAT(I)*AMP(BIN) 20 CONTINUE C C ... finally rebin to step = 1 (as in the input frame) C ... the algorithm is a pure projection, no interpolation is done C ... (since the origin and nature of the ripple are not known) C STEP = 20/PERIOD DO 40 I = 1,NPIX XLO = (I-1)*STEP + 1 XHI = XLO + STEP PIXLO = INT(XLO) PIXHI = INT(XHI) ! all into one output pixel IF ((PIXHI-PIXLO).EQ.0) THEN FRACLO = XHI - XLO OUT(I) = OUT(I) + FRACLO*SCRAT(PIXLO) ! at least two output pixels hit ELSE FRACLO = PIXLO + 1 - XLO FRACHI = XHI - PIXHI OUT(I) = OUT(I) + FRACLO*SCRAT(PIXLO) + + FRACHI*SCRAT(PIXHI) ! if in fact more than two hit DO 30 J = PIXLO + 1,PIXHI - 1 OUT(I) = OUT(I) + SCRAT(J) 30 CONTINUE END IF 40 CONTINUE C RETURN END