← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to stdred/spec/src/fripple.for

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
stdred/spec/src/fripple.for
 
1
C===========================================================================
 
2
C Copyright (C) 1995-2008 European Southern Observatory (ESO)
 
3
C
 
4
C This program is free software; you can redistribute it and/or 
 
5
C modify it under the terms of the GNU General Public License as 
 
6
C published by the Free Software Foundation; either version 2 of 
 
7
C the License, or (at your option) any later version.
 
8
C
 
9
C This program is distributed in the hope that it will be useful,
 
10
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
C GNU General Public License for more details.
 
13
C
 
14
C You should have received a copy of the GNU General Public 
 
15
C License along with this program; if not, write to the Free 
 
16
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
 
17
C MA 02139, USA.
 
18
C
 
19
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
C       Internet e-mail: midas@eso.org
 
21
C       Postal address: European Southern Observatory
 
22
C                       Data Management Division 
 
23
C                       Karl-Schwarzschild-Strasse 2
 
24
C                       D 85748 Garching bei Muenchen 
 
25
C                       GERMANY
 
26
C===========================================================================
 
27
C
 
28
      PROGRAM FRIPLE
 
29
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
C
 
31
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  11:28 - 7 DEC 1987
 
32
C
 
33
C.LANGUAGE: F77+ESOext
 
34
C
 
35
C.AUTHOR: J.D.Ponz, but I did not code it!
 
36
C
 
37
C.IDENTIFICATION
 
38
C  program    FRIPPLE              version 1.0          851030
 
39
C                                  correct some bugs    900619  MP
 
40
C
 
41
C.Keywords
 
42
C  filtering
 
43
C
 
44
C.PURPOSE
 
45
C  correct 1-D images for periodic modulation (FRIPLE)
 
46
C
 
47
C.ALGORITHM
 
48
C  Fold flux with period, calculate mean flux per phase bin, determine
 
49
C  deviation from total mean, adjust to mean level.
 
50
C
 
51
C.INPUT/OUTPUT
 
52
C the following keywords are used:
 
53
C IN_A/C/1/8           name of input image
 
54
C INPUTR/R/1/1         period
 
55
C INPUTI/I/1/2         first and last pixel to be considered
 
56
C OUT_A/C/1/8          name of ripple corrected image
 
57
C
 
58
C.VERSIOn
 
59
C 080717        last modif
 
60
C 
 
61
C --------------------------------------------------------------------
 
62
C
 
63
      IMPLICIT NONE
 
64
C
 
65
      INTEGER      MADRID
 
66
      CHARACTER*60 INIMA,OUTIMA
 
67
      CHARACTER*80 OUTEXT,HISTOR
 
68
      CHARACTER    IDENT*72,CUNIT*64
 
69
      INTEGER      NAXIS,NPIX,ISTAT,IAV
 
70
      INTEGER*8    PNTRIN,PNTRO
 
71
      INTEGER      RANGE(2),ERROR(2),NWORK
 
72
      INTEGER      NBIN,I,SCNPIX
 
73
      INTEGER*8    SCRPTR
 
74
      INTEGER      KUN,KNUL,IMIN,IMOUT
 
75
      REAL         PERIOD,AMP(500),TEST
 
76
C
 
77
      DOUBLE PRECISION DSTART(3),DSTEP(3)
 
78
C
 
79
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
80
 
 
81
      COMMON /VMR/MADRID(1)
 
82
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
83
 
 
84
C
 
85
C ... set up MIDAS environment
 
86
C
 
87
      CALL STSPRO('FRIPLE')
 
88
C
 
89
C ... get name of input frame, map it
 
90
C
 
91
      CALL STKRDC('IN_A',1,1,60,IAV,INIMA,KUN,KNUL,ISTAT)
 
92
      IF (ISTAT.NE.0) GO TO 20
 
93
      CALL STIGET(INIMA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,1,NAXIS,NPIX,
 
94
     +     DSTART,DSTEP,IDENT,CUNIT,PNTRIN,IMIN,ISTAT)
 
95
      IF (ISTAT.NE.0) GO TO 20
 
96
      IF (NAXIS.NE.1) THEN
 
97
          CALL STTPUT('ERROR: input frame must be one-dimensional',
 
98
     +                ISTAT)
 
99
          ERROR(1) = 1
 
100
          GO TO 20
 
101
      END IF
 
102
C
 
103
C ... get name of output frame, map it
 
104
C
 
105
      CALL STKRDC('OUT_A',1,1,60,IAV,OUTIMA,KUN,KNUL,ISTAT)
 
106
      CALL STIPUT(OUTIMA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,NAXIS,NPIX,
 
107
     +            DSTART,DSTEP,IDENT,CUNIT,PNTRO,IMOUT,ISTAT)
 
108
C
 
109
C  ...  get period, branch according to type (integer or real) of
 
110
C  ...  period to subroutines IAMPLT or RAMPLT to find 
 
111
C  ...  amplitude of modulation
 
112
C
 
113
      CALL STKRDR('INPUTR',1,1,IAV,PERIOD,KUN,KNUL,ISTAT)
 
114
C
 
115
      CALL STKRDI('INPUTI',1,2,IAV,RANGE,KUN,KNUL,ISTAT)
 
116
C
 
117
      TEST   = PERIOD - INT(PERIOD)
 
118
C  !   period near-integer?
 
119
      IF (TEST*NPIX/PERIOD.LT.0.05) THEN
 
120
          NBIN   = PERIOD
 
121
          CALL IAMPLT(MADRID(PNTRIN),PERIOD,RANGE,AMP)  !   no, real
 
122
      ELSE  !   use a fixed number of 20 phase bins
 
123
          NBIN   = 20
 
124
C
 
125
C ... in this case, map SCRAT file to hold intermediate results
 
126
C ... first, determine number of pixels needed
 
127
C
 
128
          SCNPIX = INT(20*NPIX/PERIOD)
 
129
          NWORK  = SCNPIX*4
 
130
          CALL TDMGET(NWORK,SCRPTR,ISTAT)
 
131
          CALL RAMPLT(MADRID(PNTRIN),SCNPIX,PERIOD,RANGE,AMP,
 
132
     +                    MADRID(SCRPTR))
 
133
      END IF
 
134
C
 
135
C ... communicate results:
 
136
C
 
137
      CALL STTPUT('Results:',ISTAT)
 
138
      DO 10 I = 1,NBIN,4
 
139
          WRITE (OUTEXT,9000) I,I + 3,AMP(I),AMP(I+1),AMP(I+2),
 
140
     +      AMP(I+3)
 
141
          CALL STTPUT(OUTEXT,ISTAT)
 
142
   10 CONTINUE
 
143
C
 
144
C ... the ripple correction is, then, done by subroutines
 
145
C ... IRIPC and RRIPC, respectively
 
146
C
 
147
      IF (TEST*NPIX/PERIOD.LT.0.05) THEN
 
148
          CALL IRIPC(MADRID(PNTRIN),NPIX,PERIOD,RANGE,AMP,
 
149
     +                 MADRID(PNTRO))
 
150
      ELSE
 
151
          CALL RRIPC(MADRID(PNTRIN),NPIX,PERIOD,RANGE,AMP,
 
152
     +                 MADRID(SCRPTR),SCNPIX,MADRID(PNTRO))
 
153
      END IF
 
154
C
 
155
C ... copy descriptor LHCUTS
 
156
C
 
157
      CALL STDCOP(IMIN,IMOUT,4,'LHCUTS',ISTAT)
 
158
C
 
159
C ... copy and update descriptor HISTOR
 
160
C
 
161
      CALL STDCOP(IMIN,IMOUT,4,'HISTORY',ISTAT)
 
162
      CALL STKRDC('HISTORY',1,1,80,IAV,HISTOR,KUN,KNUL,ISTAT)
 
163
      CALL STDWRC(IMOUT,'HISTOR',1,HISTOR,-1,80,KUN,ISTAT)
 
164
C
 
165
C ... the end
 
166
C
 
167
      CALL STSEPI
 
168
      STOP 
 
169
C
 
170
C ... direct trouble makers here:
 
171
C
 
172
   20 ERROR(1) = ISTAT
 
173
      CALL STKWRI('PROGSTAT',ERROR,1,2,KUN,ISTAT)
 
174
      CALL STSEPI
 
175
      STOP 
 
176
C
 
177
 9000 FORMAT ('Phase bins No.',I3,' to',I3,' have ','amplitudes ',
 
178
     +       3 (F7.4,1X),F7.4)
 
179
      END
 
180
C
 
181
C
 
182
C
 
183
      SUBROUTINE IAMPLT(IN,PERIOD,RANGE,AMP)
 
184
C
 
185
C
 
186
      INTEGER NCYCL,I,J,RANGE(2),IPER, IR1, IR2
 
187
      REAL PERIOD,IN(1),AMP(500),SUM(500),MEAN
 
188
C  !   consider only integer number of cycles
 
189
      NCYCL  = INT((RANGE(2)-RANGE(1)+1)/PERIOD)
 
190
C
 
191
      IPER = PERIOD
 
192
      DO 10 I = 1,IPER
 
193
          SUM(I) = 0
 
194
   10 CONTINUE
 
195
C
 
196
C ... build up flux in phase bins
 
197
C ... to compensate for slopes in input data, refer to mean 
 
198
C ... flux per cycle
 
199
C
 
200
      IR1 = RANGE(1)
 
201
      IR2 = RANGE(1) + (NCYCL-1)*PERIOD
 
202
      DO 40 I = IR1, IR2,IPER
 
203
          MEAN   = 0
 
204
          DO 20 J = 1,IPER
 
205
              MEAN   = MEAN + IN(I+J-1)
 
206
   20     CONTINUE  !   average flux in one cycle
 
207
          MEAN   = MEAN/PERIOD
 
208
          DO 30 J = 1, IPER
 
209
              SUM(J) = SUM(J) + IN(I+J-1)/MEAN
 
210
   30     CONTINUE
 
211
   40 CONTINUE
 
212
C
 
213
C  ... division into mean flux yields the relative amplitude of
 
214
C  ... each phase bin
 
215
C
 
216
      MEAN   = 0
 
217
      DO 50 I = 1,IPER
 
218
          MEAN   = MEAN + SUM(I)
 
219
   50 CONTINUE
 
220
      MEAN   = MEAN/PERIOD
 
221
C
 
222
      DO 60 I = 1,IPER
 
223
          AMP(I) = MEAN/SUM(I)
 
224
   60 CONTINUE
 
225
C
 
226
      RETURN
 
227
 
 
228
      END
 
229
C
 
230
C
 
231
C
 
232
      SUBROUTINE IRIPC(IN,NPIX,PERIOD,RANGE,AMP,OUT)
 
233
C
 
234
C
 
235
      INTEGER RANGE(2),NPIX,OFF,PHASE,I,J,INDX,IPER
 
236
      REAL IN(1),OUT(1),PERIOD,AMP(500)
 
237
C
 
238
C ... take into account possible phase offset introduced if RANGE
 
239
C ... was not defaulted to [1,NPIX]
 
240
C
 
241
      OFF    = NINT(((RANGE(1)-1.)/PERIOD-INT((RANGE(1)-1.)/PERIOD))*
 
242
     +         PERIOD)
 
243
C
 
244
      IPER = PERIOD
 
245
      DO 20 I = 1,NPIX,IPER
 
246
          DO 10 J = 1,IPER
 
247
              INDX  = I + J - 1
 
248
              IF (INDX.GT.NPIX) GO TO 30
 
249
              PHASE  = OFF + J
 
250
              IF (PHASE.GT.PERIOD) PHASE  = PHASE - PERIOD  !   do the FRIPLE compensation
 
251
              OUT(INDX) = IN(INDX)*AMP(PHASE)
 
252
   10     CONTINUE
 
253
   20 CONTINUE
 
254
C
 
255
   30 RETURN
 
256
      END
 
257
C
 
258
C
 
259
C
 
260
      SUBROUTINE RAMPLT(IN,SCNPIX,PERIOD,RANGE,AMP,SCRAT)
 
261
C
 
262
      INTEGER RANGE(2),PIXLO,PIXHI,I,J,SCNPIX,BIN,START,OFF,END
 
263
      REAL PERIOD,IN(1),AMP(500),SUM(20),FRACLO,FRACHI,XLO,XHI,MEAN,
 
264
     +     SCRAT(1),STEP
 
265
C
 
266
      DO 10 I = 1,20
 
267
          SUM(I) = 0.
 
268
   10 CONTINUE
 
269
C
 
270
C ... rebin to stepsize corresponding to one phasebin
 
271
C ... the algorithm is a pure projection, no interpolation is done
 
272
C ... (since the origin and nature of the ripple are not known)
 
273
C
 
274
      STEP   = 0.05*PERIOD
 
275
      DO 30 I = 1,SCNPIX  !   ith pixel extends from i to i+1 in both spaces
 
276
          XLO    = (I-1)*STEP + 1  !   lower and upper limit of projected output pixel
 
277
          XHI    = XLO + STEP
 
278
          PIXLO  = INT(XLO)
 
279
          PIXHI  = INT(XHI)  !   all in same output pixel
 
280
          IF ((PIXHI-PIXLO).EQ.0) THEN
 
281
              FRACLO = XHI - XLO
 
282
              SCRAT(I) = FRACLO*IN(PIXLO)  !   at least two output pixels hit
 
283
          ELSE
 
284
              FRACLO = PIXLO + 1 - XLO
 
285
              FRACHI = XHI - PIXHI
 
286
              SCRAT(I) = FRACLO*IN(PIXLO) + FRACHI*IN(PIXHI)  !   if in fact more than two
 
287
              DO 20 J = PIXLO + 1,PIXHI - 1
 
288
                  SCRAT(I) = SCRAT(I) + IN(J)
 
289
   20         CONTINUE
 
290
          END IF
 
291
 
 
292
   30 CONTINUE
 
293
C
 
294
C ... evaluate the periodicity over the given range
 
295
C ... do so over an integer number of periods
 
296
C ... compensate for a slope in the input data by averaging over
 
297
C ... one cycle of 20 pixels each
 
298
C
 
299
      START  = 20*INT((RANGE(1)-1.)/PERIOD) + 1  !   phase offset due to zeropoint used
 
300
      OFF    = 20*INT((RANGE(1)-1.)/PERIOD) + 1 - START
 
301
      END    = 20*INT((RANGE(2)-1.)/PERIOD) + 1
 
302
      DO 60 I = START,END,20
 
303
          MEAN   = 0
 
304
          DO 40 J = 1,20
 
305
              MEAN   = MEAN + SCRAT(I+J-1)
 
306
   40     CONTINUE
 
307
          MEAN   = MEAN/20
 
308
          DO 50 J = 1,20
 
309
C              BIN    = JMOD(I+J-1+OFF,20) + 1
 
310
              BIN    = MOD(I+J-1+OFF,20) + 1
 
311
              SUM(BIN) = SUM(BIN) + SCRAT(I+J-1)/MEAN
 
312
   50     CONTINUE
 
313
   60 CONTINUE
 
314
C
 
315
C  ...  division into mean value yields the relative amplitudes of 
 
316
C  ...  each each phas
 
317
C
 
318
      MEAN   = 0
 
319
      DO 70 I = 1,20
 
320
          MEAN   = MEAN + SUM(I)
 
321
   70 CONTINUE
 
322
      MEAN   = MEAN/20.
 
323
C
 
324
      DO 80 I = 1,20
 
325
          AMP(I) = MEAN/SUM(I)
 
326
   80 CONTINUE
 
327
C
 
328
      RETURN
 
329
 
 
330
      END
 
331
C
 
332
C
 
333
C
 
334
      SUBROUTINE RRIPC(IN,NPIX,PERIOD,RANGE,AMP,SCRAT,SCNPIX,OUT)
 
335
C
 
336
C
 
337
      INTEGER RANGE(2),NPIX,SCNPIX,I,J,BIN,PIXLO,PIXHI
 
338
      REAL IN(1),OUT(1),PERIOD,AMP(500),XLO,XHI,FRACLO,FRACHI,
 
339
     +     SCRAT(1),STEP
 
340
C
 
341
      DO 10 I = 1,NPIX  !   initialize output frame to zero
 
342
          OUT(I) = 0
 
343
   10 CONTINUE
 
344
C
 
345
C ... here the amplitude correction is done
 
346
C
 
347
      DO 20 I = 1,SCNPIX
 
348
C          BIN    = JMOD(I,20) + 1
 
349
          BIN    = MOD(I,20) + 1
 
350
          SCRAT(I) = SCRAT(I)*AMP(BIN)
 
351
   20 CONTINUE
 
352
C
 
353
C ... finally rebin to step = 1 (as in the input frame)
 
354
C ... the algorithm is a pure projection, no interpolation is done
 
355
C ... (since the origin and nature of the ripple are not known)
 
356
C
 
357
      STEP   = 20/PERIOD
 
358
      DO 40 I = 1,NPIX
 
359
          XLO    = (I-1)*STEP + 1
 
360
          XHI    = XLO + STEP
 
361
          PIXLO  = INT(XLO)
 
362
          PIXHI  = INT(XHI)  !   all into one output pixel
 
363
          IF ((PIXHI-PIXLO).EQ.0) THEN
 
364
              FRACLO = XHI - XLO
 
365
              OUT(I) = OUT(I) + FRACLO*SCRAT(PIXLO)  !   at least two output pixels hit
 
366
          ELSE
 
367
              FRACLO = PIXLO + 1 - XLO
 
368
              FRACHI = XHI - PIXHI
 
369
              OUT(I) = OUT(I) + FRACLO*SCRAT(PIXLO) +
 
370
     +                 FRACHI*SCRAT(PIXHI)  !   if in fact more than two hit
 
371
              DO 30 J = PIXLO + 1,PIXHI - 1
 
372
                  OUT(I) = OUT(I) + SCRAT(J)
 
373
   30         CONTINUE
 
374
          END IF
 
375
 
 
376
   40 CONTINUE
 
377
C
 
378
      RETURN
 
379
 
 
380
      END