← Back to branch summary

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

« back to all changes in this revision

Viewing changes to prim/general/src/avwndw.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:
prim/general/src/avwndw.for
 
1
C  @(#)avwndw.for       19.1 (ESO-DMD) 02/25/03 14:01:40 
 
2
C===========================================================================
 
3
C Copyright (C) 1995 European Southern Observatory (ESO)
 
4
C
 
5
C This program is free software; you can redistribute it and/or 
 
6
C modify it under the terms of the GNU General Public License as 
 
7
C published by the Free Software Foundation; either version 2 of 
 
8
C the License, or (at your option) any later version.
 
9
C
 
10
C This program is distributed in the hope that it will be useful,
 
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
C GNU General Public License for more details.
 
14
C
 
15
C You should have received a copy of the GNU General Public 
 
16
C License along with this program; if not, write to the Free 
 
17
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
21
C       Internet e-mail: midas@eso.org
 
22
C       Postal address: European Southern Observatory
 
23
C                       Data Management Division 
 
24
C                       Karl-Schwarzschild-Strasse 2
 
25
C                       D 85748 Garching bei Muenchen 
 
26
C                       GERMANY
 
27
C===========================================================================
 
28
C
 
29
      PROGRAM AVWNDW
 
30
C
 
31
C++++++++++++++++++++++++++++++++++++++++++++++++++
 
32
C
 
33
C.IDENTIFICATION
 
34
C  program AVWNDW                       version 1.00    890123
 
35
C  K. Banse                             ESO - Garching
 
36
C  1.20       890529    1.30    900130
 
37
C
 
38
C.KEYWORDS
 
39
C  bulk data frames, average
 
40
C
 
41
C.PURPOSE
 
42
C  average a no. of image frames, take only values which don't differ too much
 
43
C  from the median value (at each pixel), the result will either be a frame with
 
44
C  size equal to common area of all frames involved or size equal to union of
 
45
C  all frame areas
 
46
C
 
47
C.ALGORITHM
 
48
C  extract input frames either from given parameter string or from catalogue,
 
49
C  add up all frames + divide in the end
 
50
C  max. 20 frames may be averaged at one time !
 
51
C
 
52
C.INPUT/OUTPUT
 
53
C  the following keys are used:
 
54
C
 
55
C  ACTION/C/1/3                 = WIN, MED, MIN or MAX
 
56
C  OUT_A/C/1/60                 result frame
 
57
C  P3/C/1/80                    list of frames to be averaged
 
58
C  INPUTR/R/1/2                 BGERR,SNOISE - Poisson + Gaussian errors
 
59
C                               for WINDOW option
 
60
C
 
61
C  the following descriptors are used with WINDOW option:
 
62
C  LHCUTS/R/5/2                 low, high cuts for calculations
 
63
C  FLAT_BKG/R/1/1               background value
 
64
C  O_TIME/D/7/1                 exposure time
 
65
C       
 
66
C.VERSIONS
 
67
C 
 
68
C 010802                last modif
 
69
C       
 
70
C--------------------------------------------------
 
71
C
 
72
      IMPLICIT NONE
 
73
C
 
74
      INTEGER     BEGIN,FRMCNT,IAV,ISEQ
 
75
      INTEGER     L,LL,LP,M,IOFF
 
76
      INTEGER     N,NAXIS,NAXISC,NCHECK,NIN,NP,NPERC,NTL
 
77
      INTEGER*8   PNTRC,PNTRT
 
78
      INTEGER     SIZEC,STAT
 
79
      INTEGER     IMNOC,IMN(20)
 
80
      INTEGER     NPIX(3),NPIXC(3)
 
81
      INTEGER     OFX(20),OFY(20),POFF(20),NPNTS(20)
 
82
      INTEGER     IORD(20),NORD(20),REJ(20),TNP(4)
 
83
      INTEGER     EC,EL,ED
 
84
      INTEGER     UNI(1),NULO,MADRID(1)
 
85
      INTEGER     LEN,INDEX            !system functions
 
86
C
 
87
      CHARACTER   FRAME(20)*60,FRAMEC*60,CATFIL*60
 
88
      CHARACTER   CUNIT*64,IDENT*72,OUTPUT*80
 
89
      CHARACTER   NOCUTS*60,LINE*80,ACTION*3
 
90
      CHARACTER   ERROR1*60,ERROR2*50,ERROR3*40,ERROR4*40,ERROR6*40
 
91
      CHARACTER   ERROR7*40,ERROR9*60
 
92
C
 
93
      DOUBLE PRECISION STEP(2,20),STEPC(2)
 
94
      DOUBLE PRECISION START(2,20),STARTC(2),OSTART(2),OEND(2)
 
95
C
 
96
      REAL        EPS(3),CUTS(4)
 
97
      REAL        VAL(20),BGV(20),RMS(20),RFC(20)
 
98
      REAL        AJCUL(20),AJCUH(20),AJCU(2),SNOISE,BGERR
 
99
      REAL        INPUTR(2),EXPTI(20)
 
100
C
 
101
      DOUBLE PRECISION      EXPTIM
 
102
C
 
103
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
104
C
 
105
      COMMON   /VMR/     MADRID
 
106
      COMMON   /MYREAL/  VAL,BGV,RMS,RFC,AJCUL,AJCUH,SNOISE,BGERR,EXPTI
 
107
      COMMON   /MYINT/   IORD,NORD,REJ,TNP
 
108
C
 
109
      DATA    IDENT    /'average frame '/
 
110
      DATA    ERROR1   /'operands do not match in stepsize... '/
 
111
      DATA    ERROR2   /'operands do not overlap...'/
 
112
      DATA    ERROR3   /'stepsizes of different sign...'/
 
113
      DATA    ERROR4   /'catalogue empty...'/
 
114
      DATA    ERROR6   /'no. of input frames < 2 ...'/
 
115
      DATA    ERROR7   /'more than 20 entries in catalogue...'/
 
116
      DATA    ERROR9   /'no. of axes do not match... '/
 
117
      DATA    NOCUTS   /'LHCUTS(5,6) not there - use LHCUTS(3,4)... '/
 
118
      DATA    STARTC   /2*0.D0/,    STEPC  /2*1.D0/
 
119
C
 
120
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
121
C
 
122
C  set up MIDAS environment + enable automatic error abort
 
123
      CALL STSPRO('AVWNDW')
 
124
C 
 
125
C  init COMMON variables
 
126
      DO 50 N=1,20
 
127
         REJ(N) = 0
 
128
50    CONTINUE
 
129
      DO 60 N=1,4
 
130
         TNP(N) = 0
 
131
60    CONTINUE
 
132
C
 
133
C  get necessary keys
 
134
      CALL STKRDC('ACTION',1,1,3,IAV,ACTION,UNI,NULO,STAT)    !get option
 
135
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,UNI,NULO,STAT)    !get result frame
 
136
      CALL STKRDC('P3',1,1,80,IAV,LINE,UNI,NULO,STAT)         !list of frames
 
137
      CALL UPCAS(ACTION,ACTION)
 
138
      IF (ACTION(1:1).EQ.'W') THEN                        !only for real window option
 
139
         CALL STKRDR('INPUTR',1,2,IAV,INPUTR,UNI,NULO,STAT)
 
140
         BGERR = INPUTR(1) * INPUTR(1)                   !take square of errors
 
141
         SNOISE = INPUTR(2) * INPUTR(2)
 
142
      ENDIF
 
143
C
 
144
C  test, if we have list of frames or catalogue
 
145
      LL = INDEX(LINE,'.cat')
 
146
      IF (LL.LE.0) THEN
 
147
         LL = INDEX(LINE,'.CAT')
 
148
         IF (LL.LE.0) GOTO 500                     !No. Treat list input
 
149
      ENDIF
 
150
C
 
151
C  here we handle input from a catalog
 
152
C  get name of first image file in catalog
 
153
C
 
154
      ISEQ = 0
 
155
      CATFIL(1:) = LINE(1:)
 
156
      CALL STCGET(CATFIL,0,FRAME(1),OUTPUT,ISEQ,STAT)
 
157
      IF (FRAME(1)(1:1).EQ.' ')                 !there must be at least 1 image
 
158
     +   CALL STETER(4,ERROR4)
 
159
C
 
160
C  now loop through catalogue
 
161
      DO 200 N=2,20                                !max. 20 frames
 
162
         CALL STCGET(CATFIL,0,FRAME(N),OUTPUT,ISEQ,STAT)
 
163
         IF (FRAME(N)(1:1).EQ.' ') THEN            !indicates the end ...
 
164
            FRMCNT = N - 1
 
165
            GOTO 1000
 
166
         ENDIF
 
167
200   CONTINUE
 
168
C
 
169
      CALL STCGET(CATFIL,0,LINE,OUTPUT,ISEQ,STAT)
 
170
      IF (LINE(1:1).NE.' ')
 
171
     +   CALL STTPUT
 
172
     +   ('Max. 20 frames are used, following frames ignored...',STAT)
 
173
      FRMCNT = 20
 
174
      GOTO 1000
 
175
C
 
176
C  here we handle input from single inout line
 
177
C  pull out operands      (max. 20)
 
178
C
 
179
500   BEGIN = 1
 
180
      DO 800 N=1,21                             !max. of 20 frames...!
 
181
         CALL EXTRSS(LINE,',',BEGIN,FRAME(N),LL)
 
182
         IF (LL.LE.0) THEN
 
183
            FRMCNT = N - 1
 
184
            GOTO 1000
 
185
         ELSE
 
186
            CALL CLNFRA(FRAME(N),FRAME(N),0)
 
187
         ENDIF
 
188
800   CONTINUE
 
189
C
 
190
C  check no. of input frames
 
191
1000  IF (FRMCNT.LT.2) CALL STETER(6,ERROR6)
 
192
C
 
193
      DO 1200 N=1,FRMCNT
 
194
         DO 1100 M=1,2                        !init values...
 
195
            START(M,N) = 0.D0
 
196
            STEP(M,N) = 1.D0
 
197
            NPIX(M) = 1
 
198
1100     CONTINUE
 
199
1200  CONTINUE
 
200
C
 
201
C  get initial area from 1. frame
 
202
      DO 1400 M=1,2                  !init values...
 
203
         OSTART(M) = 0.D0
 
204
         OEND(M) = 0.D0
 
205
         NPIX(M) = 1
 
206
1400  CONTINUE
 
207
C
 
208
      CALL STFOPN(FRAME(1),D_R4_FORMAT,0,F_IMA_TYPE,IMN(1),STAT)
 
209
      CALL STDRDI(IMN(1),'NAXIS',1,1,IAV,NAXISC,UNI,NULO,STAT)
 
210
      IF (NAXISC.GT.2) THEN                !we only work on max 2-dim. frames
 
211
         NAXISC = 2
 
212
         OUTPUT(1:) = 'currently only 1 or 2-dim frames supported... '
 
213
         CALL STTPUT(OUTPUT,STAT)
 
214
      ENDIF
 
215
C
 
216
      CALL STDRDI(IMN(1),'NPIX',1,2,IAV,NPIX,UNI,NULO,STAT)
 
217
      NPNTS(1) = NPIX(1)
 
218
      CALL STDRDD(IMN(1),'START',1,2,IAV,START(1,1),UNI,NULO,STAT)
 
219
      OSTART(1) = START(1,1)
 
220
      OSTART(2) = START(2,1)
 
221
      CALL STDRDD(IMN(1),'STEP',1,2,IAV,STEP(1,1),UNI,NULO,STAT)
 
222
      STEPC(1) = STEP(1,1)
 
223
      STEPC(2) = STEP(2,1)
 
224
      CALL STDRDC(IMN(1),'CUNIT',1,1,64,IAV,CUNIT,UNI,NULO,STAT)
 
225
C
 
226
C  for "real" window option get exposure-time, low + high cuts and background
 
227
      IF (ACTION(1:1).EQ.'W') THEN
 
228
         CALL STECNT('GET',EC,EL,ED)            !disable error abort
 
229
         CALL STECNT('PUT',1,0,0)
 
230
C
 
231
         CALL STDRDD(IMN(1),'O_TIME',7,1,IAV,EXPTIM,UNI,NULO,STAT)
 
232
         IF (STAT.NE.0) THEN
 
233
            EXPTI(1) = 1.0                      !default to 1 second
 
234
         ELSE
 
235
            EXPTI(1) = EXPTIM
 
236
         ENDIF
 
237
C
 
238
C  get LHCUTS(5,6) - if not there use LHCUTS(3,4)
 
239
         CALL STDRDR(IMN(1),'LHCUTS',5,2,IAV,AJCU,UNI,NULO,STAT)
 
240
         IF (STAT.NE.0) THEN
 
241
            CALL STTPUT(NOCUTS,STAT)
 
242
            CALL STDRDR(IMN(1),'LHCUTS',3,2,IAV,AJCU,UNI,NULO,STAT)
 
243
         ENDIF
 
244
         IF (AJCU(2).LE.AJCU(1))
 
245
     +      CALL STETER(43,'invalid LHCUTS: LHCUTS(6).LE.LHCUTS(5)...')
 
246
         AJCUL(1) = AJCU(1)
 
247
         AJCUH(1) = AJCU(2)
 
248
C
 
249
         CALL STDRDR(IMN(1),'FLAT_BKG',1,1,IAV,BGV(1),UNI,NULO,STAT)
 
250
         IF (STAT.NE.0) BGV(1) = 0.             !default to background = 0.
 
251
         CALL STECNT('PUT',EC,EL,ED)
 
252
      ENDIF
 
253
C
 
254
      DO 2000 M=1,2
 
255
         EPS(M) = 0.0001 * ABS(STEPC(M))        !0.01% of stepsize as epsilon
 
256
         OEND(M) = OSTART(M) + (NPIX(M)-1)*STEPC(M)
 
257
2000  CONTINUE
 
258
C
 
259
C  now loop through the other input frames
 
260
      DO 3000 N=2,FRMCNT
 
261
         CALL STFOPN(FRAME(N),D_R4_FORMAT,0,F_IMA_TYPE,IMN(N),STAT)
 
262
C
 
263
         CALL STDRDI(IMN(N),'NAXIS',1,1,IAV,NAXIS,UNI,NULO,STAT)
 
264
         IF (NAXISC.NE.NAXIS) CALL STETER(1,ERROR9)
 
265
C
 
266
         CALL STDRDI(IMN(N),'NPIX',1,2,IAV,NPIX,UNI,NULO,STAT)
 
267
         NPNTS(N) = NPIX(1)
 
268
         CALL STDRDD(IMN(N),'START',1,2,IAV,START(1,N),UNI,NULO,STAT)
 
269
         CALL STDRDD(IMN(N),'STEP',1,2,IAV,STEP(1,N),UNI,NULO,STAT)
 
270
C
 
271
         IF (ACTION(1:1).EQ.'W') THEN
 
272
            CALL STECNT('GET',EC,EL,ED)            !disable error abort
 
273
            CALL STECNT('PUT',1,0,0)
 
274
C
 
275
            CALL STDRDD(IMN(N),'O_TIME',7,1,IAV,EXPTIM,UNI,NULO,STAT)
 
276
            IF (STAT.NE.0) THEN
 
277
               EXPTI(N) = 1.0                      !default to 1 second
 
278
            ELSE
 
279
               EXPTI(N) = EXPTIM
 
280
            ENDIF
 
281
C
 
282
            CALL STDRDR(IMN(N),'LHCUTS',5,2,IAV,AJCU,UNI,NULO,STAT)
 
283
            IF (STAT.NE.0)
 
284
     +         CALL STDRDR(IMN(N),'LHCUTS',3,2,IAV,AJCU,UNI,NULO,STAT)
 
285
            AJCUL(N) = AJCU(1)
 
286
            AJCUH(N) = AJCU(2)
 
287
C
 
288
            CALL STDRDR(IMN(N),'FLAT_BKG',1,1,IAV,BGV(N),UNI,NULO,STAT)
 
289
            IF (STAT.NE.0) BGV(N) = 0.            !default to background = 0.
 
290
            CALL STECNT('PUT',EC,EL,ED)
 
291
         ENDIF
 
292
C
 
293
C  create new result frame with dimension as intersection of input frames
 
294
C  stepsizes should have same sign and not differ too much...
 
295
         DO 2400 M=1,NAXISC
 
296
            IF (STEPC(M)*STEP(M,N).LE.0.) CALL STETER(1,ERROR3)
 
297
            IF (ABS(STEP(M,N)-STEPC(M)).GT.EPS(M))
 
298
     +         CALL STTPUT(ERROR1,STAT)                 !just warning message
 
299
2400     CONTINUE
 
300
C
 
301
C  get intersection of image areas
 
302
         DO 2600 M=1,NAXISC
 
303
            IF (STEPC(M).GT.0.) THEN
 
304
               OSTART(M) = MAX(OSTART(M),START(M,N))
 
305
               OEND(M) = MIN(OEND(M),START(M,N) + (NPIX(M)-1)*STEP(M,N))
 
306
            ELSE
 
307
               OSTART(M) = MIN(OSTART(M),START(M,N))
 
308
               OEND(M) = MAX(OEND(M),START(M,N) + (NPIX(M)-1)*STEP(M,N))
 
309
            ENDIF
 
310
2600     CONTINUE
 
311
C
 
312
C  compute individual offsets
 
313
         DO 2800 NIN=1,N
 
314
            OFX(NIN) = NINT((OSTART(1)-START(1,NIN))/STEP(1,NIN)) + 1
 
315
            OFY(NIN) = NINT((OSTART(2)-START(2,NIN))/STEP(2,NIN)) + 1
 
316
2800     CONTINUE
 
317
C
 
318
3000  CONTINUE
 
319
C
 
320
C  test, if something is left...
 
321
C
 
322
      DO 3300 M=1,NAXISC
 
323
         IF (STEPC(M)*(OEND(M)-OSTART(M)).LT.0.) CALL STETER(2,ERROR2)
 
324
3300  CONTINUE
 
325
C
 
326
C  set up standard stuff for result frame
 
327
      SIZEC = 1
 
328
      NPIXC(2) = 1
 
329
      DO 3600 M=1,NAXISC
 
330
         STARTC(M) = OSTART(M)
 
331
         NPIXC(M) = NINT((OEND(M)-OSTART(M))/STEPC(M)) + 1
 
332
         SIZEC = SIZEC * NPIXC(M)
 
333
3600  CONTINUE
 
334
      LP = NPIXC(2)
 
335
      NP = NPIXC(1)
 
336
C
 
337
C  map result frame
 
338
      CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
339
     +            NAXISC,NPIXC,STARTC,STEPC,
 
340
     +            IDENT,CUNIT,PNTRC,IMNOC,STAT)
 
341
C
 
342
C  initialize some working variables
 
343
      DO 3900 N=1,FRMCNT
 
344
         IF (ACTION(1:1).EQ.'W') RFC(N) = 1./EXPTI(N)
 
345
         POFF(N) = (OFY(N)-1)*NPNTS(N) + OFX(N)      !determine 1.pixel no.
 
346
3900  CONTINUE
 
347
C
 
348
C  and allocate working space
 
349
      NTL = NP * FRMCNT
 
350
      IF (NTL.LT.512) NTL = 512
 
351
      CALL STFXMP(NTL,D_R4_FORMAT,PNTRT,STAT)
 
352
C
 
353
C  set up counters for calculation of progress
 
354
      NTL = INT(FLOAT(LP)*.2)                  !20% of total no. of lines
 
355
      IF (NTL.LE.2) THEN
 
356
         NCHECK = LP + 1                       !not worth it...
 
357
      ELSE
 
358
         NPERC = 20
 
359
         NCHECK = NTL
 
360
      ENDIF
 
361
C
 
362
C  now, loop through lines
 
363
      DO 5000 L=1,LP
 
364
C
 
365
C  and loop through frames
 
366
         IOFF = 1
 
367
         DO 4200 N=1,FRMCNT
 
368
            CALL COPULA(IOFF,IMN(N),POFF(N),MADRID(PNTRT),NP)
 
369
            POFF(N) = POFF(N) + NPNTS(N)                  !move to next line
 
370
            IOFF = IOFF + NP
 
371
4200     CONTINUE
 
372
C
 
373
C  and do it ...
 
374
         IF (ACTION(1:1).EQ.'W') THEN
 
375
            CALL AVWIN1(MADRID(PNTRT),FRMCNT,NP)
 
376
         ELSE
 
377
            CALL AVWIN2(ACTION,MADRID(PNTRT),FRMCNT,NP)
 
378
         ENDIF
 
379
         CALL COPUL(MADRID(PNTRT),MADRID(PNTRC),L,NP)    !save resulting line
 
380
C
 
381
C  display progress
 
382
         IF (L.GE.NCHECK) CALL PROGRS(20,NTL,NPERC,NCHECK)
 
383
5000  CONTINUE
 
384
C
 
385
C  display more info
 
386
      IF (ACTION(1:1).EQ.'W') THEN
 
387
         DO 5500 N=1,FRMCNT
 
388
            WRITE(OUTPUT,10000) AJCUL(N),AJCUH(N)
 
389
            L = INDEX(FRAME(N),' ') - 1
 
390
            IF (L.LE.0) L = LEN(FRAME(N))
 
391
            LINE(1:) = FRAME(N)(1:L)//OUTPUT(1:)
 
392
            CALL STTPUT(LINE,STAT)
 
393
            WRITE(OUTPUT,10001) EXPTI(N),BGV(N),REJ(N)
 
394
            CALL STTPUT(OUTPUT,STAT)
 
395
5500     CONTINUE
 
396
         WRITE(OUTPUT,10002) TNP(1),TNP(2),TNP(3),TNP(4)
 
397
         CALL STTPUT(OUTPUT,STAT)
 
398
      ENDIF
 
399
C
 
400
C  copy descriptors + calculate new dynamic range of result frame
 
401
      CALL DSCUPT(IMN(1),IMNOC,' ',STAT)
 
402
      CALL MYMX(MADRID(PNTRC),SIZEC,CUTS(3))
 
403
      CUTS(1) = CUTS(3)
 
404
      CUTS(2) = CUTS(4)
 
405
      CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNI,STAT)
 
406
C
 
407
      CALL STSEPI
 
408
C
 
409
C  Formats
 
410
10000 FORMAT(' low and high cuts',G12.5,' ,',G12.5)
 
411
10001 FORMAT('         EXPTIM',G12.5,', BGV',G12.5,', rej_pix',I8)
 
412
10002 FORMAT('pixels with N=0:',I8,', N=1:',I8,', N=2:',I8,', N>2:',I8)
 
413
      END
 
414
 
 
415
      SUBROUTINE COPULA(IOFF,IMNO,OFFSET,A,NP)
 
416
C
 
417
      IMPLICIT NONE
 
418
C
 
419
      INTEGER     IOFF,IMNO,OFFSET,NP
 
420
      INTEGER     IAV,STAT
 
421
C
 
422
      REAL        A(1)
 
423
C
 
424
      CALL STFGET(IMNO,OFFSET,NP,IAV,A(IOFF),STAT)
 
425
      RETURN
 
426
      END
 
427
 
 
428
      SUBROUTINE COPUL(A,B,NLINE,NP)
 
429
C
 
430
      IMPLICIT NONE
 
431
C
 
432
      INTEGER     IOFF,N,NLINE,NP
 
433
C
 
434
      REAL        A(1),B(1)
 
435
C
 
436
      IOFF = (NLINE-1) * NP
 
437
      DO 800 N=1,NP
 
438
         B(N+IOFF) = A(N)
 
439
800   CONTINUE
 
440
C
 
441
      RETURN
 
442
      END
 
443
 
 
444
      SUBROUTINE AVWIN1(T,TOTAL,NP)
 
445
C
 
446
      IMPLICIT NONE
 
447
C
 
448
      REAL        T(1)
 
449
      REAL        V,VM,ER,ERR,WG,AM,ALP,F
 
450
      REAL        VAL(20),BGV(20),RMS(20),RFC(20)
 
451
      REAL        AJCUL(20),AJCUH(20),SNOISE,BGERR
 
452
      REAL        EXPTI(20)
 
453
      REAL        SQRT                        !system function
 
454
C
 
455
      LOGICAL     NOK(20)
 
456
C
 
457
      INTEGER     TOTAL,NP
 
458
      INTEGER     IORD(20),NORD(20),REJ(20),TNP(4)
 
459
      INTEGER     I,IDX,II,IM1,IO
 
460
      INTEGER     N1,N2,NH,NN,NX
 
461
      INTEGER     TNP1,TNP2,TNP3,TNP4
 
462
C
 
463
      COMMON  /MYREAL/  VAL,BGV,RMS,RFC,AJCUL,AJCUH,SNOISE,BGERR,EXPTI
 
464
      COMMON  /MYINT/   IORD,NORD,REJ,TNP
 
465
C
 
466
      EQUIVALENCE    (NORD(1),N1),(NORD(2),N2)
 
467
      EQUIVALENCE    (TNP(1),TNP1),(TNP(2),TNP2),
 
468
     +               (TNP(3),TNP3),(TNP(4),TNP4)
 
469
C
 
470
      ALP = 0.
 
471
C
 
472
C  loop through pixels in each line
 
473
      DO 5000 NX=1,NP
 
474
         NN = 0                          !init no. of valid pixels
 
475
         IDX = NX                        !init offset in working area
 
476
C
 
477
         DO 400 I=1,TOTAL
 
478
            V = T(IDX)
 
479
            NOK(I) = V.LT.AJCUL(I) .OR. V.GT.AJCUH(I)      !set NOK, if outside
 
480
C
 
481
            IF (.NOT.NOK(I)) THEN
 
482
               V = V - BGV(I)
 
483
               VAL(I) = RFC(I) * V       !fill value stack
 
484
             NN = NN + 1                 !increment no. of valid pixels
 
485
             IF (V.LT.0.) THEN
 
486
                F = 0.
 
487
             ELSE
 
488
                F = SNOISE * V
 
489
             ENDIF
 
490
             RMS(I) = RFC(I) * SQRT(BGERR + F)      !estimate deviation
 
491
            ENDIF
 
492
C
 
493
            IDX = IDX + NP                          !increase offset
 
494
400      CONTINUE
 
495
C
 
496
C  order pixel values + find median
 
497
         IF (NN.GE.3) THEN
 
498
            GOTO 2000                       !more than 2 pixels o.k.
 
499
         ELSE
 
500
            GOTO (1400,1410,1450),NN+1      !branch acc. to no. of good pixels
 
501
         ENDIF
 
502
C
 
503
C  no good pixels, take 1. frame
 
504
1400     NOK(1) = .FALSE.
 
505
         TNP1 = TNP1 + 1
 
506
         GOTO 3900
 
507
C
 
508
C  only one good pixel, take it
 
509
1410     TNP2 = TNP2 + 1
 
510
         GOTO 3900
 
511
C
 
512
C  two good pixels, test them
 
513
1450     II = 1
 
514
         DO 1600 I=1,TOTAL
 
515
            IF (.NOT.NOK(I)) THEN
 
516
               NORD(II) = I
 
517
               II = II + 1
 
518
            ENDIF
 
519
1600     CONTINUE
 
520
C
 
521
         IF (ABS(VAL(N1)-VAL(N2)).LE.RMS(N1)+RMS(N2)) THEN
 
522
            TNP3 = TNP3 + 1
 
523
         ELSE
 
524
            TNP2 = TNP2 + 1
 
525
            IF (ABS(VAL(N1)-ALP).LE.ABS(VAL(N2)-ALP)) THEN
 
526
               NOK(N2) = .TRUE.
 
527
            ELSE
 
528
               NOK(N1) = .TRUE.
 
529
            ENDIF
 
530
         ENDIF
 
531
         GOTO 3900
 
532
C
 
533
C  more than 2 pixels, compare them
 
534
2000     DO 2500 I=1,TOTAL
 
535
            IF (.NOT.NOK(I)) THEN
 
536
               V = VAL(I)
 
537
               IO = 0
 
538
               IF (I.EQ.1) GOTO 2400
 
539
               IM1 = I - 1
 
540
               DO 2200 II=1,IM1
 
541
                  IF (.NOT.NOK(II)) THEN
 
542
                     IF (V.GT.VAL(II)) THEN
 
543
                        IO = IO + 1
 
544
                     ELSE
 
545
                        IORD(II) = IORD(II) + 1
 
546
                     ENDIF
 
547
                  ENDIF
 
548
2200           CONTINUE
 
549
2400           IORD(I) = IO
 
550
            ENDIF
 
551
2500     CONTINUE
 
552
C
 
553
C  update NORD
 
554
         DO 2600 I=1,TOTAL
 
555
            IF (.NOT.NOK(I)) THEN
 
556
               II = IORD(I) + 1
 
557
               NORD(II) = I
 
558
            ENDIF
 
559
2600     CONTINUE
 
560
C
 
561
C  and get median
 
562
         NH = (NN+1) / 2
 
563
         IM1 = NH - 1
 
564
         I = NORD(NH)
 
565
         VM = VAL(I)
 
566
         ER = RMS(I)
 
567
         NH = NN
 
568
C
 
569
C  now remove all pixels with high residuals (not close to median...)
 
570
         DO 3000 I=1,IM1
 
571
            II = NORD(I)                                !start at bottom...
 
572
            ERR = RMS(II) + ER
 
573
            IF (ABS(VAL(II)-VM).GT.ERR) THEN    !discard, if value differs much
 
574
               NN = NN - 1
 
575
               NOK(II) = .TRUE.
 
576
            ENDIF
 
577
            II = NORD(NH)                               !and move to top...
 
578
            NH = NH - 1
 
579
            ERR = RMS(II) + ER
 
580
            IF (ABS(VAL(II)-VM).GT.ERR) THEN
 
581
               NN = NN - 1
 
582
               NOK(II) = .TRUE.
 
583
            ENDIF
 
584
3000     CONTINUE
 
585
C
 
586
         NN = MIN(4,NN+1)
 
587
         TNP(NN) = TNP(NN) + 1
 
588
C
 
589
C  calculate final flux
 
590
3900     WG = 0.
 
591
         AM = 0.
 
592
         DO 4000 I=1,TOTAL
 
593
            IF (.NOT.NOK(I)) THEN
 
594
               WG = WG + EXPTI(I)
 
595
               AM = AM + EXPTI(I)*VAL(I)
 
596
            ELSE
 
597
               REJ(I) = REJ(I) + 1
 
598
            ENDIF
 
599
4000     CONTINUE
 
600
         ALP = AM / WG
 
601
         T(NX) = ALP
 
602
C
 
603
5000  CONTINUE                               !end of loop through pixels
 
604
      RETURN
 
605
      END
 
606
 
 
607
      SUBROUTINE AVWIN2(ACTION,T,TOTAL,NP)
 
608
C
 
609
      IMPLICIT NONE
 
610
C
 
611
      REAL        T(1)
 
612
      REAL        VAL(20),V,O(11),RVAL
 
613
      REAL        VVAL(122)                         !just to pad the COMMON ...
 
614
C
 
615
      INTEGER     TOTAL,NP
 
616
      INTEGER     I,IDX,NN
 
617
      INTEGER     NH,NX,NH1
 
618
C
 
619
      CHARACTER   ACTION*3
 
620
C
 
621
      COMMON   /MYREAL/  VAL,VVAL
 
622
C
 
623
C ***
 
624
C
 
625
C  handle maximum search
 
626
C
 
627
C ***
 
628
C
 
629
      IF (ACTION(1:2).EQ.'MA') THEN
 
630
C
 
631
C  loop through pixels in each line
 
632
         DO 800 NX=1,NP
 
633
C
 
634
            IDX = NX                           !init offset in working area
 
635
            DO 200 I=1,TOTAL
 
636
               VAL(I) = T(IDX)
 
637
               IDX = IDX + NP                  !increase offset
 
638
200         CONTINUE
 
639
C
 
640
            V = VAL(1)
 
641
            DO 400 I=2,TOTAL
 
642
               IF (V.LT.VAL(I)) V = VAL(I)
 
643
400         CONTINUE
 
644
            T(NX) = V
 
645
800      CONTINUE
 
646
C
 
647
C ***
 
648
C
 
649
C  handle minimum search
 
650
C
 
651
C ***
 
652
C
 
653
      ELSE IF (ACTION(1:2).EQ.'MI') THEN
 
654
C
 
655
C  loop through pixels in each line
 
656
         DO 1800 NX=1,NP
 
657
C
 
658
            IDX = NX                           !init offset in working area
 
659
            DO 1200 I=1,TOTAL
 
660
               VAL(I) = T(IDX)
 
661
               IDX = IDX + NP                  !increase offset
 
662
1200        CONTINUE
 
663
C
 
664
            V = VAL(1)
 
665
            DO 1400 I=2,TOTAL
 
666
               IF (V.GT.VAL(I)) V = VAL(I)
 
667
1400        CONTINUE
 
668
            T(NX) = V
 
669
1800     CONTINUE
 
670
C
 
671
C ***
 
672
C
 
673
C  handle median search
 
674
C
 
675
C ***
 
676
C
 
677
      ELSE
 
678
C
 
679
C  loop through pixels in each line
 
680
         NH = (TOTAL+1) / 2
 
681
         NH1 = NH + 1
 
682
         DO 5000 NX=1,NP
 
683
C
 
684
            IDX = NX                           !init offset in working area
 
685
            DO 2200 I=1,TOTAL
 
686
               VAL(I) = T(IDX)
 
687
               IDX = IDX + NP                  !increase offset
 
688
2200        CONTINUE
 
689
C
 
690
C  order pixel values
 
691
            O(1) = VAL(1)
 
692
C
 
693
C  build up lower half of ordered array (including median index)
 
694
            DO 3600 I=2,NH
 
695
               NN = I
 
696
               RVAL = VAL(NN)
 
697
C
 
698
3100           IF ((NN.LE.1) .OR. (O(NN-1).LE.RVAL) ) GOTO 3200
 
699
               O(NN) = O(NN-1)
 
700
               NN = NN - 1
 
701
               GOTO 3100
 
702
C
 
703
3200           O(NN) = RVAL
 
704
3600        CONTINUE
 
705
C
 
706
C  now only merge upper half of array VAL into ordered array O if necessary
 
707
            DO 4000 I=NH1,TOTAL
 
708
               RVAL = VAL(I)
 
709
               NN = NH
 
710
               IF (RVAL.GE.O(NN)) GOTO 4000
 
711
C
 
712
3700           IF ( (NN.LE.1) .OR. (O(NN-1).LE.RVAL) ) GOTO 3800
 
713
               O(NN) = O(NN-1)
 
714
               NN = NN - 1
 
715
               GOTO 3700
 
716
C
 
717
3800           O(NN) = RVAL
 
718
4000        CONTINUE
 
719
C
 
720
            T(NX) = O(NH)
 
721
5000  CONTINUE
 
722
      ENDIF
 
723
C
 
724
      RETURN
 
725
      END