← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/echelle/src/necavmed.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/echelle/src/necavmed.for
 
1
C @(#)necavmed.for      19.1 (ES0-DMD) 02/25/03 14:20:21
 
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
C @(#)necavmed.for      19.1 (ESO)   02/25/03    14:20:21
 
30
C
 
31
C+
 
32
C.COPYRIGHT: Copyright (c) 1991 European Southern Observatory,
 
33
C                                         all rights reserved
 
34
C
 
35
C.VERSION: 1.0          23-JULY-1991
 
36
C
 
37
C.LANGUAGE: F77+ESOext
 
38
C
 
39
C.AUTHOR: P.BALLESTER
 
40
C
 
41
C.IDENTIFICATION
 
42
C
 
43
C.KEYWORDS
 
44
C
 
45
C  ECHELLE, CASPEC, BLAZE FUNCTION
 
46
C
 
47
C.PURPOSE
 
48
C
 
49
C  COMPARE SUCCESSIVE ORDERS TO GET THE MEDIAN AND AVERAGE BLAZE FUNCTION
 
50
C  
 
51
C.ALGORITHM
 
52
C
 
53
C  IN DEVELOPMENT ...
 
54
C
 
55
C.INPUT/OUTPUT
 
56
C
 
57
C-
 
58
      PROGRAM AVEMED
 
59
 
 
60
      IMPLICIT  NONE
 
61
 
 
62
      INTEGER   NAXISA,NPIXA(2),IAV,STAT,ACTVAL,MAXORD
 
63
      INTEGER*8   PNTRA,PNTRB
 
64
      INTEGER   IMNOA,IMNOB
 
65
      INTEGER   KNULL,KUNIT(1)
 
66
      INTEGER   MADRID(1)
 
67
 
 
68
      PARAMETER (MAXORD=500)
 
69
 
 
70
      CHARACTER FRAMEA*60,FRAMEB*60
 
71
      CHARACTER CUNIT*64,IDENTA*72
 
72
 
 
73
      INTEGER          ORDSTA(MAXORD),ORDEND(MAXORD),DELTAY
 
74
      REAL             THRES
 
75
      DOUBLE PRECISION STEPA(2),STARTA(2)
 
76
 
 
77
      INCLUDE 'MID_INCLUDE:st_def.inc'
 
78
 
 
79
      COMMON /VMR/ MADRID
 
80
 
 
81
      INCLUDE 'MID_INCLUDE:st_dat.inc'
 
82
 
 
83
      CALL STSPRO('AVEMED')
 
84
 
 
85
 
 
86
      CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,KUNIT,KNULL,STAT)
 
87
 
 
88
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
89
     +            2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNIT,
 
90
     +            PNTRA,IMNOA,STAT)
 
91
 
 
92
      IF (NPIXA(2).LT.MAXORD) THEN
 
93
 
 
94
          CALL STDRDI(IMNOA,'ORDSTA',1,NPIXA(2),ACTVAL,ORDSTA,
 
95
     +                                       KUNIT,KNULL,STAT)
 
96
          CALL STDRDI(IMNOA,'ORDEND',1,NPIXA(2),ACTVAL,ORDEND,
 
97
     +                                       KUNIT,KNULL,STAT)
 
98
 
 
99
      ELSE
 
100
 
 
101
          CALL STETER(10,'Buffer overflow in AVEMED.')
 
102
 
 
103
      ENDIF
 
104
 
 
105
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEB,KUNIT,KNULL,STAT)
 
106
 
 
107
      CALL STIPUT(FRAMEB,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
 
108
     +            NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNIT,
 
109
     +            PNTRB,IMNOB,STAT)
 
110
 
 
111
      CALL STKRDI('INPUTI',1,1,IAV,DELTAY,KUNIT,KNULL,STAT)
 
112
 
 
113
 
 
114
      IF (DELTAY.GT.0) THEN
 
115
 
 
116
         CALL STKRDR('INPUTR',1,1,IAV,THRES,KUNIT,KNULL,STAT)
 
117
 
 
118
         CALL MEDIAN(MADRID(PNTRA),NPIXA(1),NPIXA(2),MADRID(PNTRB),
 
119
     +            ORDSTA,ORDEND,DELTAY,THRES)
 
120
 
 
121
      ELSE 
 
122
 
 
123
         CALL AVERAGE(MADRID(PNTRA),NPIXA(1),NPIXA(2),MADRID(PNTRB),
 
124
     +            ORDSTA,ORDEND)
 
125
 
 
126
      ENDIF
 
127
 
 
128
      CALL STSEPI
 
129
 
 
130
      END
 
131
 
 
132
C ======================Median Filtering====================
 
133
 
 
134
      SUBROUTINE MEDIAN(INPFRAM,NX,NY,OUTFRAM,ORDSTA,ORDEND,
 
135
     +                  DELY,THRES)
 
136
 
 
137
      IMPLICIT   NONE
 
138
 
 
139
      INTEGER    NX,NY,DELY,ROW,COL
 
140
      INTEGER    ORDSTA(NY),ORDEND(NY)
 
141
 
 
142
      REAL       INPFRAM(NX,NY),OUTFRAM(NX,NY),THRES,RATIO
 
143
 
 
144
      INTEGER    STAMIN,STAMAX,ENDMIN,ENDMAX
 
145
      INTEGER    ROWINF,ROWSUP,FLAG,ROWCUR,Y
 
146
      REAL       TRANS,BUFFER(-50:50)
 
147
 
 
148
      IF (DELY.GT.50) THEN
 
149
 
 
150
         CALL STETER(9,'Width too big in AVEMED. Limited to 50.')
 
151
 
 
152
      ENDIF
 
153
 
 
154
C --- Initialisations
 
155
 
 
156
      STAMIN = ORDSTA(1)
 
157
      STAMAX = ORDSTA(1)
 
158
      ENDMIN = ORDEND(1)
 
159
      ENDMAX = ORDEND(1)
 
160
      ROWINF = 0
 
161
      ROWSUP = 0
 
162
 
 
163
      DO 5  ROW = 1,NY
 
164
 
 
165
             IF (ORDSTA(ROW).LT.STAMIN) STAMIN = ORDSTA(ROW)
 
166
             IF (ORDSTA(ROW).GT.STAMAX) STAMAX = ORDSTA(ROW)
 
167
             IF (ORDEND(ROW).LT.ENDMIN) ENDMIN = ORDEND(ROW)
 
168
             IF (ORDEND(ROW).GT.ENDMAX) ENDMAX = ORDEND(ROW)
 
169
 
 
170
5     CONTINUE
 
171
 
 
172
      DO 10  COL = 1,NX
 
173
 
 
174
C --- Find out what is ROWINF and ROWSUP for column COL
 
175
C --- This part is easier to program with a DO WHILE. Consider C language.
 
176
 
 
177
      IF (COL.GE.STAMAX.AND.COL.LE.ENDMIN) THEN
 
178
 
 
179
          ROWINF = 1
 
180
          ROWSUP = NY
 
181
 
 
182
      ELSE IF (COL.LT.STAMAX) THEN
 
183
 
 
184
          ROWINF = 0
 
185
          ROWSUP = 0
 
186
          DO 20  ROW = 1,NY
 
187
             IF (ORDSTA(ROW).LT.COL.AND.ROWINF.EQ.0)  ROWINF = ROW
 
188
20        CONTINUE
 
189
          DO 25  ROW = NY,1,-1
 
190
             IF (ORDSTA(ROW).LT.COL.AND.ROWSUP.EQ.0)  ROWSUP = ROW
 
191
25        CONTINUE
 
192
 
 
193
      ELSE IF (COL.GT.ENDMIN) THEN
 
194
 
 
195
          ROWINF = 0
 
196
          ROWSUP = 0
 
197
          DO 30  ROW = 1,NY
 
198
             IF (ORDEND(ROW).GT.COL.AND.ROWINF.EQ.0)  ROWINF = ROW
 
199
30        CONTINUE
 
200
          DO 35  ROW = NY,1,-1
 
201
             IF (ORDEND(ROW).GT.COL.AND.ROWSUP.EQ.0)  ROWSUP = ROW
 
202
35        CONTINUE
 
203
 
 
204
      ENDIF
 
205
 
 
206
C --- Fill in the buffer with defined values, prolonging the frame by a
 
207
C --- mirror effect.
 
208
 
 
209
      DO 40 ROW = ROWINF,ROWSUP
 
210
 
 
211
         DO 50 Y = (-1)*DELY , DELY
 
212
 
 
213
             ROWCUR = ROW+Y
 
214
 
 
215
             IF (ROWCUR.LT.ROWINF)
 
216
     +               BUFFER(Y) = INPFRAM(COL,2*ROWINF-ROWCUR)
 
217
 
 
218
             IF (ROWCUR.GT.ROWSUP)
 
219
     +           BUFFER(Y) = INPFRAM(COL,2*ROWSUP-ROWCUR)
 
220
 
 
221
             IF (ROWCUR.GE.ROWINF.AND.ROWCUR.LE.ROWSUP)
 
222
     +           BUFFER(Y) = INPFRAM(COL,ROWCUR)
 
223
 
 
224
50       CONTINUE
 
225
 
 
226
C ---    Sort the buffer values and get the median value.
 
227
C ---    Consider using HeapSort algorithm (Numerical Recipes)
 
228
 
 
229
60       CONTINUE
 
230
 
 
231
         FLAG = 0
 
232
 
 
233
         DO 70 Y = (-1)*DELY , DELY-1
 
234
 
 
235
            IF (BUFFER(Y).GT.BUFFER(Y+1)) THEN
 
236
                TRANS = BUFFER(Y)
 
237
                BUFFER(Y) = BUFFER(Y+1)
 
238
                BUFFER(Y+1) = TRANS
 
239
                FLAG = 1
 
240
            ENDIF
 
241
 
 
242
70       CONTINUE
 
243
 
 
244
         IF (FLAG.NE.0) GOTO 60
 
245
 
 
246
         RATIO = (INPFRAM(COL,ROW)-BUFFER(0))/(INPFRAM(COL,ROW)
 
247
     +                                                   +BUFFER(0))
 
248
 
 
249
         IF (ABS(RATIO).GT.THRES)   THEN
 
250
                 OUTFRAM(COL,ROW) = BUFFER(0)
 
251
         ELSE
 
252
                 OUTFRAM(COL,ROW) = INPFRAM(COL,ROW)
 
253
         ENDIF
 
254
 
 
255
40    CONTINUE
 
256
10    CONTINUE
 
257
 
 
258
      RETURN
 
259
      END
 
260
 
 
261
C ====================== Average ====================
 
262
 
 
263
      SUBROUTINE AVERAGE(INPFRAM,NX,NY,OUTFRAM,ORDSTA,ORDEND)
 
264
 
 
265
      IMPLICIT   NONE
 
266
 
 
267
      INTEGER    NX,NY,ROW,COL
 
268
      INTEGER    ORDSTA(NY),ORDEND(NY)
 
269
 
 
270
      REAL       INPFRAM(NX,NY),OUTFRAM(NX,NY)
 
271
 
 
272
      INTEGER    STAMIN,STAMAX,ENDMIN,ENDMAX
 
273
      INTEGER    ROWINF,ROWSUP,CNT
 
274
      REAL       SUM
 
275
 
 
276
      DO 5  ROW = 1,NY
 
277
 
 
278
             IF (ORDSTA(ROW).LT.STAMIN) STAMIN = ORDSTA(ROW)
 
279
             IF (ORDSTA(ROW).GT.STAMAX) STAMAX = ORDSTA(ROW)
 
280
             IF (ORDEND(ROW).LT.ENDMIN) ENDMIN = ORDEND(ROW)
 
281
             IF (ORDEND(ROW).GT.ENDMAX) ENDMAX = ORDEND(ROW)
 
282
 
 
283
5     CONTINUE
 
284
 
 
285
      DO 10  COL = 1,NX
 
286
 
 
287
C --- Find out what is ROWINF and ROWSUP for column COL
 
288
C --- This part is easier to program with a DO WHILE. Consider C language.
 
289
 
 
290
      IF (COL.GE.STAMAX.AND.COL.LE.ENDMIN) THEN
 
291
 
 
292
          ROWINF = 1
 
293
          ROWSUP = NY
 
294
 
 
295
      ELSE IF (COL.LT.STAMAX) THEN
 
296
 
 
297
          ROWINF = 0
 
298
          ROWSUP = 0
 
299
          DO 20  ROW = 1,NY
 
300
             IF (ORDSTA(ROW).LT.COL.AND.ROWINF.EQ.0)  ROWINF = ROW
 
301
20        CONTINUE
 
302
          DO 25  ROW = NY,1,-1
 
303
             IF (ORDSTA(ROW).LT.COL.AND.ROWSUP.EQ.0)  ROWSUP = ROW
 
304
25        CONTINUE
 
305
 
 
306
      ELSE IF (COL.GT.ENDMIN) THEN
 
307
 
 
308
          ROWINF = 0
 
309
          ROWSUP = 0
 
310
          DO 30  ROW = 1,NY
 
311
             IF (ORDEND(ROW).GT.COL.AND.ROWINF.EQ.0)  ROWINF = ROW
 
312
30        CONTINUE
 
313
          DO 35  ROW = NY,1,-1
 
314
             IF (ORDEND(ROW).GT.COL.AND.ROWSUP.EQ.0)  ROWSUP = ROW
 
315
35        CONTINUE
 
316
 
 
317
      ENDIF
 
318
 
 
319
C --- Fill in the buffer with defined values, prolonging the frame by a
 
320
C --- mirror effect.
 
321
 
 
322
      SUM = 0.
 
323
      DO 40 ROW = ROWINF,ROWSUP
 
324
 
 
325
          SUM = SUM + INPFRAM(COL,ROW)
 
326
 
 
327
40    CONTINUE
 
328
 
 
329
      CNT = ROWSUP - ROWINF + 1
 
330
      SUM = SUM/CNT
 
331
 
 
332
      DO 50 ROW = 1,NY
 
333
 
 
334
         OUTFRAM(COL,ROW) = SUM
 
335
 
 
336
50    CONTINUE
 
337
10    CONTINUE
 
338
 
 
339
      RETURN
 
340
      END