← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/echelle/src/necford.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/necford.for
 
1
C @(#)necford.for       19.1 (ESO-DMD) 02/25/03 14:20:23
 
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 ECHFOR
 
30
C
 
31
C+
 
32
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
33
C                                         all rights reserved
 
34
C
 
35
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  21:22 - 3 DEC 1987
 
36
C
 
37
C.LANGUAGE: F77+ESOext
 
38
C
 
39
C.AUTHOR: D.PONZ
 
40
C
 
41
C.IDENTIFICATION
 
42
C
 
43
C  program  ECHSORD  version 1.6 830310
 
44
C
 
45
C.MODIFICATIONS
 
46
C
 
47
C   910128   P. Ballester   Define and set variables for mmake -i
 
48
C
 
49
C.KEYWORDS
 
50
C
 
51
C  ECHELLE, CASPEC
 
52
C
 
53
C.PURPOSE
 
54
C
 
55
C SEARCH FOR ORDERS IN AN IMAGE. TYPE ECHELLE - FLAT FIELD
 
56
C OUTPUT IS A TABLE WITH THE FOLLOWING COLUMNS :
 
57
C COLUMN NO.  LABEL   UNIT      DESCRIPTION
 
58
C    1        :ORDER  -         ORDER NUMBER
 
59
C    2        :X      PIXEL     SAMPLE NUMBER
 
60
C    3        :Y      PIXEL     1. MOMENT AS SUM(F*Y)/SUM(F)
 
61
C    4        :YBKG   PIXEL     POSITION OF THE BACKGROUND
 
62
C    5        :BKG    DN        BACKGROUND LEVEL
 
63
C
 
64
C.ALGORITHM
 
65
C
 
66
C  MAXIMA FOLLOWING ALGORITHM IS USED
 
67
C
 
68
C  CONSTRAINS : MAXIMUM NUMBER OF ORDERS IN THE IMAGE IS 100
 
69
C               FOR EACH ORDER 19 POINTS ARE FOUND
 
70
C
 
71
C.INPUT/OUTPUT
 
72
C
 
73
C  COMMAND :
 
74
C  SEARCH/ORDER INPUT W,T,S OUTPUT METH
 
75
C  input keywords are
 
76
C  ECHC(1:8)   INPUT   - IMAGE NAME
 
77
C  ECHC(21:28)  OUTPUT  - TABLE NAME
 
78
C         W,T,S   - PARAMETERS WHERE :
 
79
C  ECHR(2)              T THRESHOLD IN THE DETECTION
 
80
C
 
81
C 010226                last modif
 
82
C 
 
83
C----------------------------------------------------------
 
84
 
 
85
      IMPLICIT  NONE
 
86
      INTEGER   MADRID
 
87
      INTEGER   ACTVAL,I
 
88
      INTEGER   IXSTR
 
89
      INTEGER   NACOL,NAXIS,NDIM,NORDER,NP,NPOINT
 
90
      INTEGER   NTOT
 
91
      INTEGER   STATUS,WINDOW,IACOL,ID
 
92
      INTEGER   NPIX(3),ICOL(7),KUN,KNUL,IMNI
 
93
      INTEGER   UPPER(100),LOWER(100),ORDPOS(100)
 
94
      INTEGER*8 IADD
 
95
      INTEGER*8 PNTRA
 
96
C 
 
97
      REAL    RPAR(4),LEVEL
 
98
C 
 
99
      DOUBLE PRECISION    START(3),STEP(3)
 
100
C 
 
101
      CHARACTER*8 FRMIN,TABLE
 
102
      CHARACTER*72 IDENA
 
103
      CHARACTER*80 LINE
 
104
      CHARACTER*64 CUNITA
 
105
      CHARACTER*2 METH
 
106
      CHARACTER*16 LABCOL(5),UNIT(5)
 
107
      CHARACTER*6 FORM(5)
 
108
C
 
109
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
110
      COMMON /VMR/MADRID(1)
 
111
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
112
C 
 
113
      DATA LABCOL(1)/'ORDER'/,FORM(1)/'I6'/
 
114
      DATA LABCOL(2)/'X'/,FORM(2)/'F7.1'/
 
115
      DATA LABCOL(3)/'Y'/,FORM(3)/'F7.1'/
 
116
      DATA LABCOL(4)/'YBKG'/,FORM(4)/'F7.1'/
 
117
      DATA LABCOL(5)/'BKG'/,FORM(5)/'E12.3'/
 
118
      DATA UNIT(1)/'UNITLESS'/
 
119
      DATA UNIT(2)/'PIXEL'/
 
120
      DATA UNIT(3)/'PIXEL'/
 
121
      DATA UNIT(4)/'PIXEL'/
 
122
      DATA UNIT(5)/'DN'/
 
123
      DATA NACOL/5/
 
124
C
 
125
      NDIM   = 3
 
126
C
 
127
C ... initialize system and read parameters
 
128
C
 
129
      CALL STSPRO('ECHFOR')
 
130
      CALL STKRDC('ECHC',1,1,8,ACTVAL,FRMIN,KUN,KNUL,STATUS)
 
131
      CALL STKRDC('ECHC',1,17,4,ACTVAL,METH,KUN,KNUL,STATUS)
 
132
      TABLE  = 'ORDER'
 
133
      CALL STKRDR('ECHR',1,4,ACTVAL,RPAR,KUN,KNUL,STATUS)
 
134
      LEVEL  = RPAR(2)
 
135
      IF (METH.EQ.'CO') THEN
 
136
          LABCOL(3) = 'YBKG'
 
137
          LABCOL(4) = 'Y'
 
138
      END IF
 
139
      CALL STTPUT(' ECHELLE DEFINITION',STATUS)
 
140
      CALL STTPUT(' ------------------',STATUS)
 
141
      CALL STTPUT(' INPUT IMAGE  : '//FRMIN,STATUS)
 
142
      CALL STTPUT(' OUTPUT TABLE : '//TABLE,STATUS)
 
143
      CALL STTPUT(' PARAMETERS  ',STATUS)
 
144
      WRITE (LINE,9010) RPAR(2)
 
145
      CALL STTPUT(LINE,STATUS)
 
146
C
 
147
C ... read frame
 
148
C
 
149
      CALL STIGET(FRMIN,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
150
     .            NDIM,NAXIS,NPIX,START,STEP,IDENA,
 
151
     .            CUNITA,PNTRA,IMNI,STATUS)
 
152
C
 
153
C ... find order positions in the middle of the image
 
154
C
 
155
      NORDER = 0
 
156
      IXSTR  = NPIX(1)/2
 
157
      WINDOW = RPAR(1)
 
158
      NP     = NPIX(2) - 2*WINDOW
 
159
      CALL FINDM1(MADRID(PNTRA),NPIX(1),NPIX(2),IXSTR,WINDOW,NP,LEVEL,
 
160
     +            NORDER,ORDPOS,UPPER,LOWER)
 
161
      IF (NORDER.EQ.0) THEN
 
162
          CALL STTPUT('Order detection failed',STATUS)
 
163
          CALL STTPUT('Use a lower threshold',STATUS)
 
164
          GO TO 20
 
165
      END IF
 
166
      WRITE (LINE,9030) NORDER
 
167
      CALL STTPUT(LINE,STATUS)
 
168
C
 
169
C ... allocate table
 
170
C
 
171
      NTOT   = (NORDER+1)*NPIX(1)/10
 
172
      IACOL  = 12
 
173
      CALL TBTINI(TABLE,F_TRANS,F_O_MODE,IACOL,NTOT,ID,STATUS)
 
174
      DO 10 I = 1,NACOL
 
175
          CALL TBCINI(ID,D_R4_FORMAT,1,FORM(I),
 
176
     .                UNIT(I),LABCOL(I),ICOL(I),
 
177
     +                STATUS)
 
178
   10 CONTINUE
 
179
      CALL TBCMAP(ID,0,IADD,STATUS)
 
180
C
 
181
C ... follow orders
 
182
C
 
183
      CALL FOLLOW(MADRID(PNTRA),NPIX(1),NPIX(2),NORDER,ORDPOS,UPPER,
 
184
     +            LOWER,WINDOW,LEVEL,MADRID(IADD),NTOT,NACOL,NPOINT)
 
185
C
 
186
C ... update number of elements
 
187
C
 
188
      NTOT   = NPOINT
 
189
      CALL STKWRI('ECHI',NORDER,4,1,KUN,STATUS)
 
190
      CALL STKWRI('ECHI',NPIX,7,2,KUN,STATUS)
 
191
      CALL TBIPUT(ID,NACOL,NTOT,STATUS)
 
192
      IDENA = 'ORDER POSITION'
 
193
C      CALL SXDPUT(TABLEFULL,'IDENT','C',IDENA,1,72,STATUS)
 
194
C      CALL DSCUPD(TABLEFULL,TABLEFULL,' ',STATUS)
 
195
      CALL TBSINI(ID,STATUS)
 
196
      CALL TBTCLO(ID,STATUS)
 
197
C
 
198
C ... free data
 
199
C
 
200
   20 CALL STSEPI
 
201
C 9000 FORMAT ('  ORDER WIDTH : ',F6.1)
 
202
 9010 FORMAT ('    THRESHOLD : ',F6.1)
 
203
C 9020 FORMAT ('        SLOPE : ',F6.3)
 
204
 9030 FORMAT (' NUMBER OF DETECTED ORDERS : ',I4)
 
205
      END
 
206
 
 
207
      SUBROUTINE FOLLOW(X,N1,N2,NORDER,ORDPOS,UPPER,LOWER,WINDOW,THRES,
 
208
     +                  TAB,NTOT,NACOL,NP)
 
209
C
 
210
C  DEFINE ORDER POSTIONS BY FOLLOWING STRUCTURES ABOVE A GIVEN THRESHOLD
 
211
C      IMPLICIT NONE
 
212
C  INPUT ARGUMENTS
 
213
C
 
214
C  X(N1, N2) REAL*4 INPUT ARRAY
 
215
C  N1, N2       INTG*4 DIMENSIONS OF X
 
216
C  NORDER INTG*4 NUMBER OF ORDERS
 
217
C  ORDPOS INTG*4 Y POSITION OF THE ORDERS IN PIXELS
 
218
C  LOWER,UPPER  INTG*4 LOWER AND UPPER ORDER POSITIONS
 
219
C  WINDOW INTG*4 NUMBER OF PIXELS IN THE EXCLUDING WINDOW
 
220
C  THRES REAL THRESHOLD
 
221
C  NTOT  INTG*4 ALLOCATED SPACE
 
222
C
 
223
C  OUTPUT ARGUMENTS
 
224
C
 
225
C  TAB   REAL ARRAY WITH THE LOCATED POSITIONS IN FORMAT
 
226
C                   (NTOT,0:NACOL) 1 - ORDER NUMBER
 
227
C               2 - SAMPLE
 
228
C       3 - MOMENT ORDER 1
 
229
C                            4 - position of back.
 
230
C       5 - BACKGROUND
 
231
C  NP  INTG NUMBER OF LOCATED POSITIONS
 
232
C
 
233
C-------------------------------------------------------------
 
234
C
 
235
C
 
236
      IMPLICIT NONE
 
237
 
 
238
      INTEGER       N1,N2,NORDER,NTOT,NP,NACOL
 
239
      INTEGER       NSAMP,ISTAT,IORD,IRL,IL
 
240
      REAL          THRES
 
241
      REAL          X(N1,N2)
 
242
      INTEGER       ORDPOS(NORDER),LOWER(NORDER),UPPER(NORDER)
 
243
      INTEGER       ICPOS(20),ILPOS(20),IUPOS(20),WINDOW
 
244
      INTEGER       TINULL
 
245
      DOUBLE PRECISION TDNULL
 
246
      REAL          TAB(NTOT,0:NACOL),TRNULL
 
247
      CHARACTER*80  LINE
 
248
      LOGICAL       NEXT
 
249
      INTEGER       NF,IHW,IXSTR,IBACK,ISTEP,IX,IY,I1,IU
 
250
      INTEGER       KZ6902,NPTS,NFOUND
 
251
C
 
252
      CALL TBMNUL(TINULL,TRNULL,TDNULL)
 
253
      NSAMP  = 10
 
254
      IXSTR = N1/2
 
255
C
 
256
C ... DISPLAY RESULTS
 
257
C
 
258
      CALL STTPUT(' ',ISTAT)
 
259
      WRITE (LINE,9020)
 
260
      CALL STTPUT(LINE,ISTAT)
 
261
      WRITE (LINE,9030)
 
262
      CALL STTPUT(LINE,ISTAT)
 
263
C
 
264
C ... ITERATION ON ORDERS
 
265
C
 
266
      NP     = 0
 
267
      DO 40 IORD = 1,NORDER
 
268
C
 
269
C ... INCLUDE START + END POINTS
 
270
C
 
271
          NF     = 0
 
272
          IF (IORD.NE.NORDER) IHW    = (ORDPOS(IORD+1)-ORDPOS(IORD))/2
 
273
          NP     = NP + 1
 
274
          TAB(NP,1) = IORD
 
275
          TAB(NP,2) = 1.
 
276
          TAB(NP,3) = TRNULL
 
277
          TAB(NP,4) = TRNULL
 
278
          TAB(NP,5) = TRNULL
 
279
          NP     = NP + 1
 
280
          TAB(NP,1) = IORD
 
281
          TAB(NP,2) = N1
 
282
          TAB(NP,3) = TRNULL
 
283
          TAB(NP,4) = TRNULL
 
284
          TAB(NP,5) = TRNULL
 
285
          NP     = NP + 1
 
286
          TAB(NP,1) = IORD
 
287
          TAB(NP,2) = IXSTR
 
288
          TAB(NP,3) = ORDPOS(IORD)
 
289
          IBACK  = ORDPOS(IORD) + IHW
 
290
          IF (IBACK.LT.N2-WINDOW) THEN
 
291
              TAB(NP,4) = IBACK
 
292
              TAB(NP,5) = X(IXSTR,IBACK)
 
293
 
 
294
          ELSE
 
295
              TAB(NP,4) = TRNULL
 
296
              TAB(NP,5) = TRNULL
 
297
          END IF
 
298
C
 
299
C ... LEFT/RIGHT  PART OF THE IMAGE
 
300
C
 
301
          DO 30 IRL = 1,2
 
302
              IF (IRL.EQ.1) THEN
 
303
                  ISTEP  = 1
 
304
 
 
305
              ELSE
 
306
                  ISTEP  = -1
 
307
              END IF
 
308
 
 
309
              IX     = IXSTR
 
310
              IL     = LOWER(IORD)
 
311
              IU     = UPPER(IORD)
 
312
              NEXT   = .TRUE.
 
313
99            KZ6902 = 1
 
314
                  IF ( .NOT. (NEXT)) GO TO 20
 
315
                  IX     = IX + ISTEP
 
316
                  IY     = MAX(IL-3,WINDOW)
 
317
                  I1     = MIN(N2-WINDOW,IU+3)
 
318
                  NPTS   = I1 - IY + 1
 
319
                  CALL FINDM1(X,N1,N2,IX,IY,NPTS,THRES,NFOUND,ICPOS,
 
320
     +                        IUPOS,ILPOS)
 
321
C
 
322
C ... CHECK IF ONE ORDER FOUND
 
323
C
 
324
                  IF (NFOUND.EQ.1) THEN
 
325
C
 
326
C ... STORE INTO TABLE
 
327
C
 
328
                      IF (MOD(IX,NSAMP).EQ.0) THEN
 
329
                          NF     = NF + 1
 
330
                          NP     = NP + 1
 
331
                          TAB(NP,1) = IORD
 
332
                          TAB(NP,2) = IX
 
333
                          TAB(NP,3) = ICPOS(1)
 
334
                          IBACK  = ICPOS(1) + IHW
 
335
                          IF (IBACK.LT.N2-WINDOW) THEN
 
336
                              TAB(NP,4) = IBACK
 
337
                              TAB(NP,5) = X(IXSTR,IBACK)
 
338
 
 
339
                          ELSE
 
340
                              TAB(NP,4) = TRNULL
 
341
                              TAB(NP,5) = TRNULL
 
342
                          END IF
 
343
 
 
344
                      END IF
 
345
 
 
346
                      IL     = ILPOS(1)
 
347
                      IU     = IUPOS(1)
 
348
                      NEXT   = IL .GT. WINDOW .AND. IU .LT.
 
349
     +                         (N2-WINDOW) .AND. IX .GT. WINDOW .AND.
 
350
     +                         IX .LT. (N1-WINDOW)
 
351
 
 
352
                  ELSE
 
353
                      NEXT   = .FALSE.
 
354
                  END IF
 
355
 
 
356
              GOTO 99
 
357
   20         CONTINUE
 
358
   30     CONTINUE
 
359
          WRITE (LINE,9010) IORD,IXSTR,ORDPOS(IORD),IHW,NF
 
360
          CALL STTPUT(LINE,ISTAT)
 
361
   40 CONTINUE
 
362
      WRITE (LINE,9040)
 
363
      CALL STTPUT(LINE,ISTAT)
 
364
      RETURN
 
365
 
 
366
C 9000 FORMAT (I4)
 
367
 9010 FORMAT (1X,I7,2X,I8,2X,I8,2X,I10,2X,I8)
 
368
 9020 FORMAT (' SEQ.NO.  X CENTER  Y CENTER  INTERORDER  TEMPLA')
 
369
 9030 FORMAT (' -------  --------  --------  ----------  --------')
 
370
 9040 FORMAT (' -------------------------------------------------')
 
371
      END
 
372
 
 
373
      SUBROUTINE FINDM1(X,N1,N2,IX,IY,NP,THRES,NORDER,ORDPOS,UPPER,
 
374
     +                  LOWER)
 
375
C
 
376
      IMPLICIT NONE
 
377
C 
 
378
C  FIND CENTER OF THE ORDERS IN THE MIDDLE OF THE IMAGE
 
379
C
 
380
C  INPUT ARGUMENTS
 
381
C  X(N1,N2) REAL INPUT ARRAY
 
382
C  N1,N2 INTG ARRAY DIMENSION
 
383
C  IX, IY INTG    STARTING PIXEL POSITIONS
 
384
C  NP  INTG    NUMBER OF PIXELS TO SEARCH
 
385
C  THRES REAL THRESHOLD
 
386
C
 
387
C  OUTPUT ARGUMENTS
 
388
C  NORDER INTG NUMBER OD ORDERS
 
389
C  ORDPOS INTG    POSITION OF THE ORDER (UPPER-LOWER)/2
 
390
C  UPPER INTG UPPER ORDER LIMIT
 
391
C  LOWER INTG LOWER ORDER LIMIT
 
392
C
 
393
      INTEGER    N1,N2,IX,IY,NP,NORDER,ISTART,NO,I,NFIRST,NPIX
 
394
      INTEGER ORDPOS(1),UPPER(1),LOWER(1),WIDTH,WIDTH1,STATUS
 
395
C 
 
396
      REAL       THRES
 
397
      REAL X(N1,N2)
 
398
C 
 
399
      LOGICAL FIRST
 
400
C
 
401
C 
 
402
C 
 
403
      ISTART = IX
 
404
      WIDTH  = 0
 
405
C
 
406
C ... FIND UPPER LEVEL
 
407
C
 
408
      FIRST  = .TRUE.
 
409
      NO     = 0
 
410
      DO 10, I = IY,IY + NP - 1
 
411
          IF (X(ISTART,I).GT.THRES) THEN
 
412
C
 
413
C ... ABOVE THE THRESHOLD
 
414
C
 
415
              IF (FIRST) THEN
 
416
C
 
417
C ... FIRST VALUE ABOVE THE THRESHOLD
 
418
C
 
419
                  FIRST  = .FALSE.
 
420
                  NFIRST = I
 
421
                  NPIX   = 1
 
422
              END IF
 
423
 
 
424
          ELSE
 
425
C
 
426
C ... BELOW THE THRESHOLD
 
427
C
 
428
              IF ( .NOT. FIRST) THEN
 
429
C
 
430
C ... FIRST VALUE BELOW THE THRESHOLD
 
431
C
 
432
                  FIRST  = .TRUE.
 
433
                  NO     = NO + 1
 
434
                  UPPER(NO) = I - 1
 
435
                  LOWER(NO) = NFIRST
 
436
                  ORDPOS(NO) = LOWER(NO) + (UPPER(NO)-LOWER(NO))/2
 
437
                  IF (WIDTH.EQ.0) THEN
 
438
                      WIDTH  = UPPER(NO) - LOWER(NO)
 
439
 
 
440
                  ELSE
 
441
                      WIDTH1 = UPPER(NO) - LOWER(NO)
 
442
                      IF (ABS(WIDTH1-WIDTH).GT.0.1*WIDTH) 
 
443
     +                   CALL STTPUT(
 
444
     +                  'Warning: Order width changes',STATUS)
 
445
                      WIDTH  = WIDTH1
 
446
                  END IF
 
447
 
 
448
              END IF
 
449
 
 
450
          END IF
 
451
 
 
452
   10 CONTINUE
 
453
      NORDER = NO
 
454
      RETURN
 
455
 
 
456
      END