← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/ccdred/src/ccdillum.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/ccdred/src/ccdillum.for
 
1
C @(#)ccdillum.for      19.1 (ES0-DMD) 02/25/03 14:16:38
 
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 Massachusetss Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Corresponding 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 MKILLU
 
30
C++++++++++++++++++++++++++++++++++++++++++++++++++
 
31
C.IDENTIFICATION: ccdillum.FOR
 
32
C.KEYWORDS:       box smoothing
 
33
C.PURPOSE:        box smoothing with growing boxcar size
 
34
C.VERSION:        930810 RHW Creation
 
35
C--------------------------------------------------
 
36
      IMPLICIT     NONE
 
37
 
 
38
      INTEGER      IAV
 
39
      INTEGER      STATUS,MADRID
 
40
      INTEGER      NAXIS,NPIX(2)
 
41
      INTEGER      IMNOC,IMNOA
 
42
      INTEGER*8    PNTRA,PNTRC
 
43
      INTEGER      KUN,KNUL
 
44
      INTEGER      XBMIN, XBMAX
 
45
      INTEGER      YBMIN, YBMAX
 
46
C
 
47
      CHARACTER*60 FRAMEA,FRAMEC
 
48
      CHARACTER*72 IDENT
 
49
      CHARACTER*48 CUNIT
 
50
      CHARACTER*3  ACTION
 
51
C
 
52
      DOUBLE PRECISION START(2),STEP(2)
 
53
C
 
54
      REAL         AVER
 
55
      REAL         XBOX(2), YBOX(2)
 
56
      REAL         XBMINR, XBMAXR
 
57
      REAL         YBMINR, YBMAXR
 
58
      REAL         CLIP(2)
 
59
      REAL         LOWSIG,HIGSIG
 
60
C
 
61
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
62
      COMMON/VMR/MADRID(1)
 
63
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
64
C
 
65
C *** do the work
 
66
      CALL STSPRO('ILLUM')
 
67
 
 
68
      CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,KUN,KNUL,STATUS)
 
69
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
70
     1                2,NAXIS,NPIX,START,STEP,IDENT,CUNIT,
 
71
     2                PNTRA,IMNOA,STATUS)
 
72
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,KUN,KNUL,STATUS)
 
73
      CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
74
     1            NAXIS,NPIX,START,STEP,IDENT,CUNIT,PNTRC,
 
75
     2            IMNOC,STATUS)
 
76
C
 
77
C *** smoothing boxes
 
78
      CALL STKRDR('XBOX',1,2,IAV,XBOX,KUN,KNUL,STATUS)
 
79
      CALL STKRDR('YBOX',1,2,IAV,YBOX,KUN,KNUL,STATUS)
 
80
      XBMINR = XBOX(1)
 
81
      XBMAXR = XBOX(2)
 
82
      YBMINR = YBOX(1)
 
83
      YBMAXR = YBOX(2)
 
84
C
 
85
C *** check the smoothing boxes
 
86
      IF (XBMINR.LT.1.) THEN
 
87
         XBMINR = XBMINR*NPIX(1)
 
88
      ENDIF 
 
89
      IF (XBMAXR.LT.1.) THEN
 
90
         XBMAXR = XBMAXR*NPIX(1)
 
91
      ENDIF      
 
92
      IF (YBMINR.LT.1.) THEN
 
93
         YBMINR = YBMINR*NPIX(2)
 
94
      ENDIF      
 
95
      IF (YBMAXR.LT.1.) THEN
 
96
         YBMAXR = YBMAXR*NPIX(2)
 
97
      ENDIF
 
98
C
 
99
      XBMIN = INT(MAX(2,MIN(NPIX(1),NINT(MIN(XBMINR,XBMAXR)))))
 
100
      XBMAX = INT(MAX(2,MIN(NPIX(1),NINT(MAX(XBMINR,XBMAXR)))))
 
101
      YBMIN = INT(MAX(2,MIN(NPIX(2),NINT(MIN(YBMINR,YBMAXR)))))
 
102
      YBMAX = INT(MAX(2,MIN(NPIX(2),NINT(MAX(YBMINR,YBMAXR)))))
 
103
C
 
104
C *** to clip or not to clip
 
105
      CALL STKRDC('CLIP',1,1,3,IAV,ACTION,KUN,KNUL,STATUS)
 
106
      CALL UPCAS(ACTION,ACTION)
 
107
      IF (ACTION(1:1).EQ.'Y') THEN
 
108
         CALL STKRDR('SIGMA',1,2,IAV,CLIP,KUN,KNUL,STATUS)
 
109
         LOWSIG = CLIP(1)
 
110
         HIGSIG = CLIP(2)
 
111
      ENDIF
 
112
C
 
113
C *** test on clipping
 
114
      IF (ACTION(1:1).EQ.'Y') THEN
 
115
         CALL ILLUM(PNTRA,PNTRC,NPIX,XBMIN,XBMAX,YBMIN,YBMAX,
 
116
     2              LOWSIG,HIGSIG,AVER)
 
117
      ELSE
 
118
         CALL QILLUM(PNTRA,PNTRC,NPIX,XBMIN,XBMAX,YBMIN,YBMAX,AVER)
 
119
      ENDIF
 
120
C
 
121
      CALL STDCOP(IMNOA,IMNOC,1,' ',STATUS) 
 
122
      CALL STDWRR(IMNOC,'CCDMEAN',AVER,1,1,KUN,STATUS)
 
123
      CALL STSEPI
 
124
      END
 
125
 
 
126
 
 
127
      SUBROUTINE ILLUM(IN,OUT,NPIX,XBMIN,XBMAX,YBMIN,YBMAX,
 
128
     2                 LOW,HIGH,MEAN)
 
129
C
 
130
      IMPLICIT NONE
 
131
 
 
132
      INTEGER       IN, OUT
 
133
      INTEGER       NPIX(2)
 
134
      INTEGER       XBMIN,XBMAX,YBMIN,YBMAX
 
135
      REAL          LOW, HIGH
 
136
      REAL          MEAN
 
137
C
 
138
      INTEGER       LINEIN, LINEOUT
 
139
      INTEGER       IMNOS, IMNOG, IMNOP
 
140
      INTEGER*8     IPOINT, PTRS, SUM, AVG, OUTPUT
 
141
      INTEGER*8     PTR
 
142
      INTEGER       I, NITER, NREJ
 
143
      INTEGER       MADRID(1)
 
144
      INTEGER       NCOLS,NLINES
 
145
      INTEGER       STAT, IAV
 
146
      INTEGER       YBOX2
 
147
 
 
148
      REAL          SCALE
 
149
      REAL          ILLSUM
 
150
      EXTERNAL      ILLSUM
 
151
      INTEGER       BXCLN
 
152
      EXTERNAL      BXCLN
 
153
      INCLUDE       'MID_INCLUDE:ST_DEF.INC'
 
154
      COMMON        /VMR/MADRID
 
155
      INCLUDE       'MID_INCLUDE:ST_DAT.INC'
 
156
      DATA          NITER/10/
 
157
C
 
158
C *** here create vitual memory
 
159
      NCOLS  = NPIX(1)
 
160
      NLINES = NPIX(2)
 
161
C      
 
162
      CALL STFCRE('dummsumm',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
 
163
     2            NCOLS,IMNOS,STAT)                          ! create sum array
 
164
      CALL STFMAP(IMNOS,F_X_MODE,1,NCOLS,IAV,SUM,STAT)
 
165
C
 
166
      CALL STFCRE('dummavg',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
 
167
     2            NCOLS,IMNOG,STAT)                        !create output array
 
168
      CALL STFMAP(IMNOG,F_X_MODE,1,NCOLS,IAV,AVG,STAT)
 
169
      OUTPUT = AVG
 
170
C
 
171
      CALL STFCRE('dummpntr',D_I4_FORMAT,F_X_MODE,F_IMA_TYPE,
 
172
     2            YBMAX,IMNOP,STAT)                      ! create pointer array
 
173
      CALL STFMAP(IMNOP,F_X_MODE,1,YBMAX,IAV,PTRS,STAT)
 
174
C
 
175
C *** now start the work
 
176
C      IF (YBMAX.LT.NLINES) THEN
 
177
C         YBMAX = NLINES
 
178
C      ENDIF
 
179
C
 
180
C *** accumulate the minimum y box
 
181
      CALL ILLZER(MADRID(SUM),NCOLS)                   ! initialize sum array
 
182
C
 
183
C *** acumulate the minimum y box
 
184
      LINEIN = 0
 
185
   10 CONTINUE
 
186
      IF (LINEIN.LT.YBMIN) THEN
 
187
         LINEIN = LINEIN + 1
 
188
         IPOINT = IN + (LINEIN-1)*NCOLS
 
189
         CALL ILLADD(MADRID(IPOINT),MADRID(SUM),MADRID(SUM),NCOLS)
 
190
         PTR = PTRS + MOD(LINEIN,YBMAX) 
 
191
         MADRID(PTR) = IPOINT
 
192
         GOTO 10
 
193
      ENDIF
 
194
C
 
195
C *** output the minimum y box
 
196
      YBOX2 = YBMIN
 
197
      SCALE = FLOAT(YBMIN)
 
198
      CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
 
199
C
 
200
C *** Iteratively clean the initial lines
 
201
      PTR = PTRS
 
202
      IF (YBOX2 .NE. YBMAX) THEN
 
203
         PTR = PTR+1
 
204
      ENDIF
 
205
C
 
206
      DO I = 1, NITER
 
207
         NREJ = 0
 
208
         DO LINEOUT = 1,LINEIN
 
209
            IPOINT = MADRID(PTR+LINEOUT-1)
 
210
            NREJ   = NREJ + BXCLN(MADRID(IPOINT),MADRID(AVG),
 
211
     2                            MADRID(SUM),NCOLS,LOW,HIGH)
 
212
         ENDDO
 
213
         IF (NREJ.GT.0) THEN
 
214
            CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
 
215
         ELSE
 
216
            GOTO 20
 
217
         ENDIF
 
218
      ENDDO
 
219
C
 
220
   20 CONTINUE
 
221
      YBOX2 = (YBMIN+1)/2
 
222
      LINEOUT = 0
 
223
          
 
224
   30 CONTINUE
 
225
      IF (LINEOUT.LT.YBOX2) THEN
 
226
         LINEOUT = LINEOUT + 1
 
227
         IPOINT  = OUT + (LINEOUT-1)*NCOLS
 
228
         CALL COPYF(MADRID(OUTPUT),MADRID(IPOINT),NCOLS)
 
229
         GO TO 30
 
230
      ENDIF
 
231
      MEAN = YBOX2*ILLSUM(MADRID(OUTPUT),NCOLS)            ! calculate the MEAN
 
232
C
 
233
C *** increase the y box size by steps of 2 until the maximum szie.
 
234
   40 CONTINUE
 
235
      IF (LINEIN.LT.YBMAX) THEN
 
236
         LINEIN = LINEIN + 1
 
237
         IPOINT = IN + (LINEIN-1)*NCOLS
 
238
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
239
         PTR = PTRS + MOD(LINEIN,YBMAX)
 
240
         MADRID(PTR) = IPOINT
 
241
         SCALE = SCALE + 1
 
242
C
 
243
         NREJ  = BXCLN(MADRID(IPOINT),MADRID(AVG),
 
244
     2                 MADRID(SUM),NCOLS,LOW,HIGH)
 
245
         CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
 
246
C
 
247
         LINEIN = LINEIN + 1
 
248
         IPOINT = IN + (LINEIN-1)*NCOLS
 
249
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
250
         PTR    = PTRS + MOD(LINEIN,YBMAX)
 
251
         MADRID(PTR) = IPOINT
 
252
C
 
253
         NREJ  = BXCLN(MADRID(IPOINT),MADRID(AVG),MADRID(SUM),
 
254
     2           NCOLS,LOW,HIGH)
 
255
         SCALE = SCALE + 1
 
256
         CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
 
257
C
 
258
         LINEOUT = LINEOUT + 1
 
259
         IPOINT = OUT + (LINEOUT-1)*NCOLS
 
260
         CALL COPYF(MADRID(AVG),MADRID(IPOINT),NCOLS)
 
261
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)            ! calculate MEAN
 
262
         GOTO 40
 
263
      ENDIF
 
264
C
 
265
C *** for each line subtract the last line from the sum, add the 
 
266
C     next line to the sum, and output a line
 
267
   50 CONTINUE
 
268
      IF (LINEIN.LT.NLINES) THEN
 
269
         LINEIN = LINEIN + 1
 
270
         PTR = PTRS + MOD(LINEIN,YBMAX)
 
271
         IPOINT = MADRID(PTR)
 
272
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
273
         IPOINT = IN + (LINEIN-1)*NCOLS
 
274
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
275
         MADRID(PTR) = IPOINT
 
276
C
 
277
         NREJ  = BXCLN(MADRID(IPOINT),MADRID(AVG),MADRID(SUM),
 
278
     2           NCOLS,LOW,HIGH)
 
279
C
 
280
         LINEOUT = LINEOUT + 1
 
281
         IPOINT = OUT + (LINEOUT-1)*NCOLS
 
282
         CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
 
283
         CALL COPYF(MADRID(AVG),MADRID(IPOINT),NCOLS)
 
284
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)        ! calculate the MEAN
 
285
         GOTO 50
 
286
      ENDIF
 
287
C
 
288
C *** decrease the y box in steps of 2 until mimimum y box
 
289
   60 CONTINUE
 
290
      IF (LINEOUT.LT.(NLINES-YBOX2)) THEN
 
291
         LINEIN = LINEIN + 1
 
292
         PTR = PTRS + MOD(LINEIN,YBMAX)
 
293
         IPOINT = MADRID(PTR)
 
294
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
295
         LINEIN = LINEIN + 1
 
296
         PTR = PTRS + MOD(LINEIN,YBMAX)
 
297
         IPOINT = MADRID(PTR)
 
298
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
299
         SCALE = SCALE-2
 
300
C
 
301
         LINEOUT  = LINEOUT + 1
 
302
         IPOINT = OUT + (LINEOUT-1)*NCOLS
 
303
         CALL AGBCAR(MADRID(SUM),MADRID(IPOINT),NCOLS,XBMIN,XBMAX,SCALE)
 
304
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)        ! calculate the MEAN
 
305
         GOTO 60
 
306
      ENDIF
 
307
C
 
308
C *** output the last liens of the minimum y box size
 
309
      CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
 
310
   70 CONTINUE
 
311
      IF (LINEOUT.LT.NLINES) THEN
 
312
         LINEOUT = LINEOUT + 1
 
313
         IPOINT = OUT + (LINEOUT-1)*NCOLS
 
314
         CALL COPYF(MADRID(OUTPUT),MADRID(IPOINT),NCOLS)
 
315
         GOTO 70
 
316
      ENDIF
 
317
C
 
318
C *** write the scaling vector
 
319
      MEAN = MEAN + (NLINES-LINEOUT)*ILLSUM(MADRID(AVG),NCOLS)
 
320
      MEAN = MEAN/(NCOLS*NLINES)
 
321
C
 
322
      RETURN
 
323
      END 
 
324
 
 
325
 
 
326
 
 
327
      SUBROUTINE QILLUM(IN,OUT,NPIX,XBMIN,XBMAX,YBMIN,YBMAX,MEAN)
 
328
C
 
329
      IMPLICIT NONE
 
330
 
 
331
      INTEGER       IN, OUT
 
332
      INTEGER       NPIX(2)
 
333
      INTEGER       XBMIN,XBMAX,YBMIN,YBMAX
 
334
      REAL          MEAN
 
335
C
 
336
      INTEGER       LINEIN, LINEOUT
 
337
      INTEGER       IMNOS, IMNOO, IMNOP
 
338
      INTEGER       IPOINT, PTRS, SUM, OUTPUT
 
339
      INTEGER       PTR
 
340
      INTEGER       MADRID(1)
 
341
      INTEGER       NCOLS,NLINES
 
342
      INTEGER       STAT, IAV
 
343
      INTEGER       YBOX1
 
344
 
 
345
      REAL          SCALE
 
346
      REAL          ILLSUM
 
347
      EXTERNAL      ILLSUM
 
348
      INCLUDE       'MID_INCLUDE:ST_DEF.INC'
 
349
      COMMON        /VMR/MADRID
 
350
      INCLUDE       'MID_INCLUDE:ST_DAT.INC'
 
351
C
 
352
C *** here create vitual memory
 
353
      NCOLS  = NPIX(1)
 
354
      NLINES = NPIX(2)
 
355
C      
 
356
      CALL STFCRE('dummsumm',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
 
357
     2            NCOLS,IMNOS,STAT)                          ! create sum array
 
358
      CALL STFMAP(IMNOS,F_X_MODE,1,NCOLS,IAV,SUM,STAT)
 
359
C
 
360
      CALL STFCRE('dummoutp',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
 
361
     2            NCOLS,IMNOO,STAT)                        !create output array
 
362
      CALL STFMAP(IMNOO,F_X_MODE,1,NCOLS,IAV,OUTPUT,STAT)
 
363
C
 
364
      CALL STFCRE('dummpntr',D_I4_FORMAT,F_X_MODE,F_IMA_TYPE,
 
365
     2            YBMAX,IMNOP,STAT)                      ! create pointer array
 
366
      CALL STFMAP(IMNOP,F_X_MODE,1,YBMAX,IAV,PTRS,STAT)
 
367
C
 
368
C *** now start the work
 
369
C      IF (YBMAX.LT.NLINES) THEN
 
370
C         YBMAX = NLINES
 
371
C      ENDIF
 
372
C
 
373
C *** accumulate the minimum y box
 
374
      CALL ILLZER(MADRID(SUM),NCOLS)                   ! initialize sum array
 
375
C
 
376
C *** acumulate the minimum y box
 
377
      LINEIN = 0
 
378
   10 CONTINUE
 
379
      IF (LINEIN.LT.YBMIN) THEN
 
380
         LINEIN = LINEIN + 1
 
381
         IPOINT = IN + (LINEIN-1)*NCOLS
 
382
         CALL ILLADD(MADRID(IPOINT),MADRID(SUM),MADRID(SUM),NCOLS)
 
383
         PTR = PTRS + MOD(LINEIN,YBMAX) 
 
384
         MADRID(PTR) = IPOINT
 
385
         GOTO 10
 
386
      ENDIF
 
387
C
 
388
C *** output the minimum y box
 
389
      YBOX1 = (YBMIN+1)/2
 
390
      SCALE = FLOAT(YBMIN)
 
391
      CALL AGBCAR(MADRID(SUM),MADRID(OUTPUT),NCOLS,XBMIN,XBMAX,SCALE)
 
392
      LINEOUT = 0
 
393
   20 CONTINUE
 
394
      IF (LINEOUT.LT.YBOX1) THEN
 
395
         LINEOUT = LINEOUT + 1
 
396
         IPOINT  = OUT + (LINEOUT-1)*NCOLS
 
397
         CALL COPYF(MADRID(OUTPUT),MADRID(IPOINT),NCOLS)
 
398
         GO TO 20
 
399
      ENDIF
 
400
      MEAN = YBOX1*ILLSUM(MADRID(OUTPUT),NCOLS)            ! calculate the MEAN
 
401
C
 
402
C *** increase the y box size by steps of 2 until the maximum szie.
 
403
   30 CONTINUE
 
404
      IF (LINEIN.LT.YBMAX) THEN
 
405
         LINEIN = LINEIN + 1
 
406
         IPOINT = IN + (LINEIN-1)*NCOLS
 
407
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
408
         PTR    = PTRS + MOD(LINEIN,YBMAX)
 
409
         MADRID(PTR) = IPOINT
 
410
         LINEIN = LINEIN + 1
 
411
         IPOINT = IN + (LINEIN-1)*NCOLS
 
412
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
413
         PTR    = PTRS + MOD(LINEIN,YBMAX)
 
414
         MADRID(PTR) = IPOINT
 
415
C
 
416
         SCALE = SCALE + 2
 
417
         LINEOUT = LINEOUT + 1
 
418
         IPOINT = OUT + (LINEOUT-1)*NCOLS
 
419
         CALL AGBCAR(MADRID(SUM),MADRID(IPOINT),NCOLS,XBMIN,XBMAX,SCALE)
 
420
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)            ! calculate MEAN
 
421
         GOTO 30
 
422
      ENDIF
 
423
C
 
424
C *** for each line subtract the last line from the sum, add the 
 
425
C     next line to the sum, and output a line
 
426
   40 CONTINUE
 
427
      IF (LINEIN.LT.NLINES) THEN
 
428
         LINEIN = LINEIN + 1
 
429
         PTR = PTRS + MOD(LINEIN,YBMAX)
 
430
         IPOINT = MADRID(PTR)
 
431
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
432
         IPOINT = IN + (LINEIN-1)*NCOLS
 
433
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
434
         MADRID(PTR) = IPOINT
 
435
C
 
436
         LINEOUT = LINEOUT + 1
 
437
         IPOINT = OUT + (LINEOUT-1)*NCOLS
 
438
         CALL AGBCAR(MADRID(SUM),MADRID(IPOINT),NCOLS,XBMIN,XBMAX,SCALE)
 
439
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)        ! calculate the MEAN
 
440
         GOTO 40
 
441
      ENDIF
 
442
C
 
443
C *** decrease the y box in steps of 2 until mimimum y box
 
444
   50 CONTINUE
 
445
      IF (LINEOUT.LT.(NLINES-YBOX1)) THEN
 
446
         LINEIN = LINEIN + 1
 
447
         PTR = PTRS + MOD(LINEIN,YBMAX)
 
448
         IPOINT = MADRID(PTR)
 
449
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
450
         LINEIN = LINEIN + 1
 
451
         PTR = PTRS + MOD(LINEIN,YBMAX)
 
452
         IPOINT = MADRID(PTR)
 
453
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
 
454
C
 
455
         LINEOUT  = LINEOUT + 1
 
456
         SCALE = SCALE -2
 
457
         IPOINT = OUT + (LINEOUT-1)*NCOLS
 
458
         CALL AGBCAR(MADRID(SUM),MADRID(IPOINT),NCOLS,XBMIN,XBMAX,SCALE)
 
459
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)        ! calculate the MEAN
 
460
         GOTO 50
 
461
      ENDIF
 
462
C
 
463
C *** output the last liens of the minimum y box size
 
464
      CALL AGBCAR(MADRID(SUM),MADRID(OUTPUT),NCOLS,XBMIN,XBMAX,SCALE)
 
465
   60 CONTINUE
 
466
      IF (LINEOUT.LT.NLINES) THEN
 
467
         LINEOUT = LINEOUT + 1
 
468
         IPOINT = OUT + (LINEOUT-1)*NCOLS
 
469
         CALL COPYF(MADRID(OUTPUT),MADRID(IPOINT),NCOLS)
 
470
         GOTO 60
 
471
      ENDIF
 
472
C
 
473
C *** write the scaling vector
 
474
      MEAN = MEAN + (NLINES-LINEOUT)*ILLSUM(MADRID(OUTPUT),NCOLS)
 
475
      MEAN = MEAN/(NCOLS*NLINES)
 
476
C
 
477
      RETURN
 
478
      END 
 
479
C
 
480
C
 
481
C
 
482
      INTEGER FUNCTION BXCLN(DATA,BXAVG,SUM,NCOLS,LOW,HIGH)
 
483
C
 
484
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
485
C   Reject data values from thr sum for the next boxcar average which
 
486
C   exceed the minimum and maximum residual values from the current
 
487
C   boxcar average. This excludes data from the moving avarge before it 
 
488
C   enters the average.
 
489
C-----------------------------------------------------------------------
 
490
C
 
491
      IMPLICIT NONE
 
492
   
 
493
      REAL     DATA(1)
 
494
      REAL     BXAVG(1)
 
495
      REAL     SUM(1)
 
496
      INTEGER  NCOLS
 
497
      REAL     LOW
 
498
      REAL     HIGH
 
499
C
 
500
      INTEGER  I,NREJ
 
501
      REAL     RMS,RESID,MINRES,MAXRES
 
502
C
 
503
      RMS    = 0.0
 
504
      DO I = 1,NCOLS
 
505
         RMS = RMS + (DATA(I)-BXAVG(1))**2
 
506
      ENDDO
 
507
      RMS    = SQRT(RMS/NCOLS)
 
508
      MINRES = -LOW*RMS
 
509
      MAXRES = HIGH*RMS
 
510
C
 
511
      NREJ = 0
 
512
      DO I = 1,NCOLS
 
513
         RESID = DATA(I)-BXAVG(I)
 
514
         IF ((RESID.LT.MINRES) .OR. (RESID.GT.MAXRES)) THEN
 
515
            DATA(I) = BXAVG(I)
 
516
            SUM(I)  = SUM(I)-RESID
 
517
            NREJ    = NREJ+1
 
518
         ENDIF
 
519
      ENDDO
 
520
C
 
521
      BXCLN = NREJ
 
522
      RETURN
 
523
      END
 
524
C
 
525
C
 
526
C
 
527
      SUBROUTINE ILLZER(A,NC)
 
528
C
 
529
      IMPLICIT NONE
 
530
      REAL     A(1)
 
531
      INTEGER  NC
 
532
      INTEGER  IC
 
533
C
 
534
      DO IC = 1,NC
 
535
         A(IC) = 0.0
 
536
      ENDDO
 
537
      RETURN
 
538
      END
 
539
C
 
540
C
 
541
C
 
542
      SUBROUTINE ILLADD(A,B,C,NC)
 
543
C
 
544
      IMPLICIT NONE
 
545
      REAL     A(1),B(1),C(1)
 
546
      INTEGER  NC
 
547
      INTEGER  IC
 
548
C
 
549
      DO IC = 1,NC
 
550
         C(IC) = A(IC)+B(IC)
 
551
      ENDDO
 
552
      RETURN
 
553
      END
 
554
C
 
555
C
 
556
C
 
557
      SUBROUTINE ILLSUB(A,B,C,NC)
 
558
C
 
559
      IMPLICIT NONE
 
560
      REAL     A(1),B(1),C(1)
 
561
      INTEGER  NC
 
562
      INTEGER  IC
 
563
C
 
564
      DO IC = 1,NC
 
565
         C(IC) = A(IC)-B(IC)
 
566
      ENDDO
 
567
      RETURN
 
568
      END
 
569
C
 
570
C
 
571
C
 
572
      REAL FUNCTION ILLSUM(A,NC)
 
573
C
 
574
      IMPLICIT NONE
 
575
      REAL     A(1)
 
576
      INTEGER  NC
 
577
      INTEGER  IC
 
578
C
 
579
      ILLSUM = 0.0       
 
580
      DO IC = 1,NC
 
581
         ILLSUM = ILLSUM + A(IC)
 
582
      ENDDO
 
583
 
 
584
      RETURN
 
585
      END
 
586
C
 
587
C
 
588
C
 
589
 
 
590
      SUBROUTINE AGBCAR(INA,OUTA,NCOLS,XBMIN,XBMAX,YBOX)
 
591
C 
 
592
C +++ Vector growing boxcar smooth. Taken from the IRAF CCD package
 
593
C
 
594
      REAL     INA(1)
 
595
      REAL     OUTA(1)
 
596
      INTEGER  NCOLS
 
597
      INTEGER  XBMIN,XBMAX
 
598
      REAL     YBOX
 
599
C
 
600
      INTEGER  COLIN, COLOUT, LASTCOL
 
601
      INTEGER  NPIX, XBMIN2
 
602
      REAL     SUM,OUTPUT
 
603
C
 
604
      XBMIN2 = (XBMIN+1)/2
 
605
      COLIN  = 0 
 
606
      SUM    = 0.0
 
607
   10 CONTINUE
 
608
      IF (COLIN.LT.XBMIN) THEN
 
609
         COLIN = COLIN + 1
 
610
         SUM   =  SUM + INA(COLIN)
 
611
         GOTO 10
 
612
      ENDIF
 
613
C
 
614
      NPIX   = XBMIN * YBOX
 
615
      OUTPUT = SUM/NPIX
 
616
      COLOUT = 0
 
617
   20 CONTINUE
 
618
      IF (COLOUT.LT.XBMIN2) THEN
 
619
         COLOUT       = COLOUT + 1
 
620
         OUTA(COLOUT) = OUTPUT
 
621
         GO TO 20
 
622
      ENDIF 
 
623
C
 
624
   30 CONTINUE
 
625
      IF (COLIN.LT.XBMAX) THEN
 
626
         COLIN  = COLIN + 1
 
627
         SUM    = SUM + INA(COLIN)
 
628
         COLIN  = COLIN + 1
 
629
         SUM    = SUM + INA(COLIN)
 
630
         NPIX   = NPIX + 2*YBOX
 
631
         COLOUT = COLOUT + 1
 
632
         OUTA(COLOUT) = SUM/NPIX
 
633
         GO TO 30
 
634
      ENDIF
 
635
C
 
636
      LASTCOL = 0
 
637
   40 CONTINUE
 
638
      IF (COLIN.LT.NCOLS) THEN
 
639
         COLIN   = COLIN + 1
 
640
         LASTCOL = LASTCOL +1
 
641
         SUM     = SUM + INA(COLIN) - INA(LASTCOL)
 
642
         COLOUT  = COLOUT + 1
 
643
         OUTA(COLOUT) = SUM/NPIX
 
644
         GO TO 40
 
645
      ENDIF
 
646
C
 
647
   50 CONTINUE
 
648
      IF (COLOUT.LT. NCOLS - XBMIN2) THEN
 
649
         LASTCOL     = LASTCOL + 1
 
650
         SUM         = SUM - INA(LASTCOL)
 
651
         LASTCOL     = LASTCOL + 1
 
652
         SUM         = SUM - INA(LASTCOL)
 
653
         NPIX        = NPIX - 2*YBOX
 
654
         COLOUT      = COLOUT + 1
 
655
         OUTA(COLOUT)= SUM/NPIX
 
656
         GO TO 50
 
657
      ENDIF
 
658
C
 
659
      OUTPUT = SUM/NPIX
 
660
   70 CONTINUE
 
661
      IF (COLOUT.LT.NCOLS) THEN
 
662
         COLOUT = COLOUT + 1
 
663
         OUTA(COLOUT) = OUTPUT
 
664
         GO TO 70
 
665
      ENDIF
 
666
 
 
667
      RETURN
 
668
      END
 
669
 
 
670
 
 
671
 
 
672