← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/long/src/spoext.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/long/src/spoext.for
 
1
C @(#)spoext.for        19.1 (ES0-DMD) 02/25/03 14:25:24
 
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 @(#)spoext.for        19.1  (ESO)  02/25/03  14:25:24 
 
30
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
31
C .COPYRIGHT   (C) 1993 European Southern Observatory    
 
32
C .IDENT       .prg                                      
 
33
C .AUTHORS     Pascal Ballester (ESO/Garching)           
 
34
C              Cristian Levin   (ESO/La Silla)           
 
35
C .KEYWORDS    Spectroscopy, Long-Slit                   
 
36
C .PURPOSE                                               
 
37
C .VERSION     1.0  Package Creation  17-MAR-1993        
 
38
C -------------------------------------------------------
 
39
 
 
40
C @(#)spoext.for        4.2 (ESO-IPG) 11/20/92 14:47:36
 
41
          PROGRAM SPOEXT
 
42
C+
 
43
C
 
44
C.IDENTIFICATION
 
45
C
 
46
C  PROGRAM  SPOEXT                VERSION 1.0        880420
 
47
C                                        
 
48
C  Otmar Stahl                        Landessternwarte Heidelberg
 
49
C     M. Peron  , modif suggested by T. Oosterloo                                      900830
 
50
C.KEYWORDS
 
51
C
 
52
C  SPECTRA, CCD, DATA EXTRACTION
 
53
C
 
54
C.PURPOSE
 
55
C
 
56
C  EXTRACT STELLAR SPECTRA FROM A CCD IMAGE. THE INDIVIDUAL ROWS OF 
 
57
C  THE CCD IMAGE WHICH CONTAIN THE FLUX ARE ADDED. THE WEIGHTS ARE
 
58
C  CHOSEN SUCH THAT THE SUMMED SPECTRUM HAS THE OPTIMAL S/N-RATIO.
 
59
C  COSMIC RAY HITS ARE LOCATED BY ANALYSING THE PROFILE PERPENDICULAR 
 
60
C  TO THE DISPERSION AND MASKED OUT. 
 
61
C  
 
62
C.ALGORITHM
 
63
C
 
64
C  THE METHOD IS DESCRIBED IN DETAIL IN THE PAPER BY 
 
65
C  Keith Horne in PASP 98, pg. 609. 
 
66
C  THE MASK USED TO REMOVE THE COSMIC RAY HITS IS KEPT IN THE
 
67
C  MASK FILE cosmics.bdf.
 
68
C
 
69
C  EXTRACT/SPECTRUM  output input sky order,iter read-out-noise,g,thresh
 
70
C
 
71
C  OUT_A                output spectrum
 
72
C  IN_A                 input image (must be bias and dark subtracted 
 
73
C                       (if relevant) and divided by a normalized flat 
 
74
C                       field, not sky subtracted)
 
75
C  IN_B                 sky image (normally derived from IN_A in advance)
 
76
C  INPUTI               fit order for the polynomials, number of 
 
77
C                       iterations for cosmic ray removal
 
78
C  INPUTR                read-out-noise, conversion factor (e/ADU),
 
79
C                       threshold for cosmic ray removal
 
80
C
 
81
      IMPLICIT NONE
 
82
      CHARACTER*60    INFR,SKYFR,OUTFR
 
83
      CHARACTER       CUNIT*64,IDENT*72
 
84
      CHARACTER       CUNITS*64,IDENTS*72
 
85
      CHARACTER*80    OUTPUT
 
86
      INTEGER       NPIX(2),MPIX(2),IFIT(2),NUM
 
87
      INTEGER       IAV,ISTAT,NAXIS
 
88
      INTEGER*8     IPNTRA,IPNTRB,IPNTRC
 
89
      INTEGER       NBYTES
 
90
      INTEGER*8     IPNTR1,IPNTR2,IPNTR3,IPNTR4,IPNTR5,IPNTR6
 
91
      INTEGER       KUN,KNUL,IDI,IDS,IDO,MADRID,IDM,ISS
 
92
      REAL          CCD(3)
 
93
      DOUBLE PRECISION  START(2),STEP(2),STARTS(2),STEPS(2)
 
94
      INCLUDE       'MID_INCLUDE:ST_DEF.INC'
 
95
      COMMON/VMR/MADRID(1)
 
96
      INCLUDE       'MID_INCLUDE:ST_DAT.INC'
 
97
C
 
98
      CALL STSPRO('SPOEXT')
 
99
C
 
100
C get input and map images
 
101
C
 
102
      CALL STKRDC('IN_A',  1,1,60,IAV,INFR,KUN,KNUL,ISTAT)
 
103
      CALL STKRDC('IN_B',  1,1,60,IAV,SKYFR,KUN,KNUL,ISTAT)
 
104
      CALL STKRDC('OUT_A', 1,1,60,IAV,OUTFR,KUN,KNUL,ISTAT)
 
105
      CALL STKRDI('PUTI',1,2,IAV,IFIT,KUN,KNUL,ISTAT)
 
106
      CALL STKRDR('PUTR',1,3,IAV,CCD,KUN,KNUL,ISTAT)
 
107
C
 
108
      CALL STIGET(INFR,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
109
     .2,NAXIS,NPIX,START,STEP,IDENT,CUNIT,IPNTR1,IDI,ISTAT)
 
110
C
 
111
      CALL STIGET(SKYFR,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
112
     .2,NAXIS,MPIX,STARTS,STEPS,IDENTS,CUNITS,IPNTR6,IDS,ISTAT)
 
113
C
 
114
      IF ((MPIX(1).NE.NPIX(1)).OR.(MPIX(2).NE.NPIX(2))) THEN
 
115
          WRITE(OUTPUT,100) 
 
116
100       FORMAT('OBJECT AND SKY FRS TO NOT MATCH')
 
117
          CALL STTPUT(OUTPUT,ISTAT)
 
118
          GO TO 999
 
119
      ENDIF
 
120
C
 
121
      CALL STIPUT(OUTFR,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
122
     .1,NPIX(1),START(1),STEP(1),IDENT,CUNIT,IPNTR2,IDO,ISTAT)
 
123
C
 
124
      CALL STIPUT('cosmics',D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
125
     .NAXIS,NPIX,START,STEP,IDENT,CUNIT,IPNTR3,IDM,ISTAT)
 
126
C
 
127
C reserve  space for arrays
 
128
C
 
129
      NUM = NPIX(1)*NPIX(2)*4
 
130
      NBYTES = NPIX(1)*4
 
131
      CALL TDMGET(NUM,IPNTR4,ISS)     
 
132
      CALL TDMGET(NUM,IPNTR5,ISS)     
 
133
      CALL TDMGET(NBYTES,IPNTRA,ISS)
 
134
      CALL TDMGET(NBYTES,IPNTRB,ISS)
 
135
      CALL TDMGET(NBYTES,IPNTRC,ISS)
 
136
C
 
137
C DO THE REAL WORK
 
138
C
 
139
      CALL XTRACT (MADRID(IPNTR1),MADRID(IPNTR2),MADRID(IPNTR3), 
 
140
     &MADRID(IPNTR4),MADRID(IPNTR5),MADRID(IPNTR6),MADRID(IPNTRA),
 
141
     &MADRID(IPNTRB),MADRID(IPNTRC),NPIX,IFIT,CCD)
 
142
C
 
143
C COPY HISTORY
 
144
C
 
145
C      CALL COPYDSC_ST(INFR,OUTFR,4,'HISTORY',ISTAT)
 
146
C
 
147
      CALL TDMFRE(NUM,IPNTR4,ISS)     
 
148
      CALL TDMFRE(NUM,IPNTR5,ISS)     
 
149
      CALL TDMFRE(NBYTES,IPNTRA,ISS)
 
150
      CALL TDMFRE(NBYTES,IPNTRB,ISS)
 
151
      CALL TDMFRE(NBYTES,IPNTRC,ISS)
 
152
999   CALL STSEPI
 
153
      END
 
154
C
 
155
C**************************************************************************
 
156
      SUBROUTINE XTRACT(IN,OUT,MASK,PROF,VARI,SKY,X,Y,Z,NPIX,IFIT,CCD)
 
157
C 
 
158
C THIS ROUTINE DOES THE EXTRACTION OF THE SPECTRUM
 
159
C
 
160
      IMPLICIT NONE
 
161
C      INTEGER        MASK(1) ! DECLARED AS INTEGER*2 !!!
 
162
      REAL            MASK(1) ! DECLARED AS INTEGER*2 !!!
 
163
      INTEGER        NPIX(2),IFIT(2)
 
164
      INTEGER        J,K,KK,L,II,M,N,IORD,ITER,IWIND
 
165
      REAL           IN(1),OUT(1),PROF(1),VARI(1),SKY(1)
 
166
      REAL           A(10),B(10),S(10),W(10)     
 
167
      REAL           RON,G,SIGMA,V0,SUM,CHISQ,YY,DIFF,DIFF2,QUOT
 
168
      REAL           SUM1,SUM2,RMAX
 
169
      REAL           X(1),Y(1),Z(1),CCD(3)
 
170
C                                                                       
 
171
      M = NPIX(1)
 
172
      N = NPIX(2)
 
173
C
 
174
      IORD = IFIT(1)
 
175
      ITER = IFIT(2)
 
176
      IWIND = 3
 
177
C
 
178
      RON   = CCD(1)
 
179
      G     = CCD(2)
 
180
      SIGMA = CCD(3)
 
181
C
 
182
      SIGMA = SIGMA*SIGMA
 
183
      V0 = RON**2/G**2
 
184
C
 
185
C  EXTRACT INITIAL SPECTRUM (JUST THE SUM) AND SET mask TO 1
 
186
C
 
187
      DO 20 J = 1,M
 
188
        SUM = 0.
 
189
        DO 10 K = 1,N
 
190
          II = (K-1)*M+J
 
191
          MASK(II) = 1
 
192
          SUM = SUM+IN(II)-SKY(II)
 
193
10      CONTINUE
 
194
        X(J) = J
 
195
        OUT(J) = SUM
 
196
20    CONTINUE
 
197
C
 
198
C COMPUTE SPATIAL PROFILE, I.E. THE 2-DIM SPECTRUM 
 
199
C IS DIVIDED (ROW BY ROW) BY THE SUM
 
200
C                                                 
 
201
      DO 40 J = 1,M
 
202
        DO 30 K = 1,N
 
203
          II = (K-1)*M+J
 
204
          IF (OUT(J).EQ.0) OUT(J)=1
 
205
          PROF(II) = (IN(II)-SKY(II))/OUT(J)
 
206
30      CONTINUE
 
207
40    CONTINUE
 
208
C
 
209
C  FILTER SPATIAL PROFILE AND FIT WITH POLYNOMIAL
 
210
C
 
211
      DO 70 K = 1,N
 
212
        DO 50 J = 1,M
 
213
          II = (K-1)*M+J
 
214
          Y(J) = PROF(II)
 
215
50      CONTINUE
 
216
C
 
217
C  MEDIAN FILTER OF ROWS
 
218
C
 
219
        CALL FILTER(Y,Z,M,IWIND)
 
220
C 
 
221
C  FIT WITH POLYNOMIAL 
 
222
C
 
223
        CALL LSORTH(X,Y,A,B,S,W,M,CHISQ,IORD)
 
224
        DO 60 J = 1,M
 
225
          II = (K-1)*M+J
 
226
          CALL POLY(X(J),YY,A,B,S,W,IORD)
 
227
          PROF(II) = YY
 
228
60      CONTINUE
 
229
70    CONTINUE
 
230
C
 
231
C NORMALIZE SPATIAL PROFILE AND ENFORCE POSITIVITY
 
232
C
 
233
      DO 100 J = 1,M
 
234
        SUM=0.
 
235
        DO 80 K = 1,N
 
236
          II = (K-1)*M+J
 
237
          IF(PROF(II).LT.0.) PROF(II)=0.
 
238
          SUM = SUM+PROF(II)
 
239
80      CONTINUE
 
240
        DO 90 K = 1,N
 
241
          II = (K-1)*M+J
 
242
          PROF(II) = PROF(II)/SUM
 
243
90      CONTINUE
 
244
100   CONTINUE
 
245
C
 
246
C ITERATE (in the present version, the spatial profile is not 
 
247
C iterated, but only the mask and the variance)
 
248
C                                                                       
 
249
      DO 160 L = 1,ITER                                                  
 
250
C                                                                       
 
251
C  compute variance estimate                 
 
252
C
 
253
        DO 120 J = 1,M                                               
 
254
          DO 110 K = 1,N                  
 
255
           II = (K-1)*M+J
 
256
           VARI(II) = V0 + ABS(OUT(J)*PROF(II)+SKY(II))/G
 
257
110       CONTINUE
 
258
120     CONTINUE
 
259
C                                                                       
 
260
C MASK COSMIC RAY HITS (ONLY ONE PER COLUMN AND ITERATION)
 
261
C
 
262
        DO 150 J = 1,M
 
263
          RMAX = 1.
 
264
          KK = 0
 
265
          DO 130 K = 1,N
 
266
            II = (K-1)*M+J                                                     
 
267
            DIFF = IN(II)-SKY(II)-OUT(J)*PROF(II)
 
268
            DIFF2 = DIFF*DIFF
 
269
            QUOT = (DIFF2/(SIGMA*VARI(II)))*MASK(II)
 
270
            IF(QUOT.GT.RMAX) THEN
 
271
              RMAX = QUOT
 
272
              KK = K
 
273
            ENDIF
 
274
130       CONTINUE
 
275
          IF(KK.NE.0) THEN
 
276
            II = (KK-1)*M+J
 
277
            MASK(II) = 0
 
278
          ENDIF
 
279
C
 
280
C COMPUTE WEIGHTED MEAN WITH COSMICS MASKED
 
281
C
 
282
          SUM1 = 0.
 
283
          SUM2 = 0.            
 
284
          DO 140 K = 1,N
 
285
            II = (K-1)*M+J
 
286
            SUM1 = SUM1 + (MASK(II)*PROF(II)*(IN(II)-SKY(II)))/VARI(II)
 
287
            SUM2 = SUM2 + (MASK(II)*PROF(II)*PROF(II))        /VARI(II)
 
288
140       CONTINUE
 
289
          OUT(J) = SUM1/SUM2
 
290
150     CONTINUE
 
291
160   CONTINUE
 
292
C                                                                       
 
293
      RETURN
 
294
      END                                                           
 
295
C
 
296
C****************************************************************************
 
297
      SUBROUTINE FILTER(ARR,Z,M,IWIND)
 
298
C
 
299
      IMPLICIT NONE
 
300
      INTEGER     II,IWIND,J,M,K
 
301
      REAL        ARR(1),Z(1),Y(50),XMED
 
302
C
 
303
      II = 2*IWIND+1
 
304
C
 
305
C  SORT DATA IN ARRAY X AND CALL MEDIAN FILTER
 
306
C
 
307
      DO 20 J = IWIND+1,M-IWIND
 
308
        DO 10 K = 1,II
 
309
          Y(K) = ARR(J+K-IWIND-1)
 
310
10      CONTINUE
 
311
C
 
312
        CALL MDIAN1(Y,II,XMED)
 
313
C
 
314
        Z(J) = XMED
 
315
20    CONTINUE
 
316
C
 
317
C  COPY DATA
 
318
C
 
319
      DO 30 J = IWIND+1,M-IWIND
 
320
         ARR(J) = Z(J) 
 
321
30    CONTINUE
 
322
C
 
323
      RETURN
 
324
      END
 
325
C
 
326
C**************************************************************************
 
327
      SUBROUTINE MDIAN1(ARR,N,XMED)
 
328
C
 
329
C COMPUTE MEDIAN VALUE OF AN ARRAY
 
330
C
 
331
      IMPLICIT NONE
 
332
      INTEGER        N,N2
 
333
      REAL           ARR(N),XMED
 
334
C
 
335
C  SORT DATA
 
336
C
 
337
      CALL SORT(N,ARR)
 
338
C
 
339
C  COMPUTE MEDIAN
 
340
C
 
341
      N2 = N/2
 
342
      IF(2*N2.EQ.N) THEN
 
343
        XMED = 0.5*(ARR(N2)+ARR(N2+1))
 
344
      ELSE
 
345
        XMED = ARR(N2+1)
 
346
      ENDIF
 
347
C
 
348
      RETURN
 
349
      END
 
350
C
 
351
C************************************************************************
 
352
      SUBROUTINE SORT(N,ARR)
 
353
C
 
354
C SORTS AN ARRAY
 
355
C 
 
356
      IMPLICIT NONE
 
357
      INTEGER     N,I,J
 
358
      REAL        A,ARR(N)
 
359
C
 
360
      DO 20 J = 2,N
 
361
        A = ARR(J)      
 
362
        DO 5 I = J-1,1,-1
 
363
          IF(ARR(I).LE.A) GOTO 10
 
364
          ARR(I+1) = ARR(I)
 
365
5       CONTINUE
 
366
      I = 0
 
367
10    ARR(I+1) = A
 
368
20    CONTINUE
 
369
C
 
370
      RETURN
 
371
      END
 
372
C
 
373
C*****************************************************************************
 
374
C        POLY : INPUT: XFIT AND POLYNOMIAL COEFFICIENTS COMPUTED
 
375
C              WITH LSORTH
 
376
C               OUTPUT: YFIT
 
377
C                                                                      
 
378
        SUBROUTINE POLY(XFIT,YFIT,A,B,S,W,IORD)  
 
379
C
 
380
        IMPLICIT NONE
 
381
        INTEGER           K,IORD            
 
382
        REAL              A(10),B(10),S(10),W(10),P(10) 
 
383
        REAL              XFIT,YFIT
 
384
C   
 
385
        P(1)=1.  
 
386
        YFIT=S(1) 
 
387
C     
 
388
        P(2)=XFIT-A(1)   
 
389
        YFIT=YFIT+S(2)*P(2) 
 
390
C  
 
391
        DO 10 K=3,IORD+1   
 
392
          P(K)=(XFIT-A(K-1))*(P(K-1))-B(K-1)*P(K-2)
 
393
          YFIT=YFIT+S(K)*P(K)   
 
394
10      CONTINUE   
 
395
        RETURN     
 
396
        END 
 
397
C
 
398
C***************************************************************************   
 
399
C        LSORTH : FITS ORTHOGONAL POLYNOMIAL OF MAXIMUM DEGREE IORD
 
400
C                 THROUGH X-Y TABLE. AN F-TEST DECIDES WHETHER THE FIT
 
401
C                 STOPS ALREADY AT A LOWER ORDER THAN IORD.
 
402
C                 THE ARRAY Y IS RETURNED AS A VECTOR
 
403
C                 CONTAINING THE RESIDUALS.
 
404
C                 THE ROUTINE IS DESCRIBED IN PASP 91 (1979), PG. 546.    
 
405
C 
 
406
C                 THE POLYNOMIAL COEFFICIENTS ARE RETURNED
 
407
C
 
408
        SUBROUTINE LSORTH (X,Y,A,B,S,W,M,CHISQ,IORD)
 
409
C
 
410
        IMPLICIT  NONE
 
411
        INTEGER      IMAX,IMAXST,IFAUTO,IORD,I,N,M,I1,IFF,K,K1
 
412
        INTEGER      IFTEST
 
413
        REAL         X(1),Y(1),A(10),B(10),S(10),W(10),P(10)
 
414
        REAL         CHISQ,DEL,XNU,F,F95,DCHI,P2
 
415
C
 
416
        DATA IMAXST/10/,IFTEST/2/ 
 
417
        IMAX=IMAXST     
 
418
        IFAUTO=0   
 
419
        IF(IORD .LT. 10) GOTO 110
 
420
        IORD=9
 
421
        IFAUTO=1
 
422
 110    CONTINUE 
 
423
        IF(IORD .NE. 0) IMAX=MAX0(IABS(IORD)+1,2) 
 
424
        DO 10 I=1,IMAXST
 
425
        W(I)=0. 
 
426
        S(I)=0. 
 
427
        A(I)=0. 
 
428
10      B(I)=0.
 
429
        P(1)=1.
 
430
        DO 20 N=1,M
 
431
        W(1)=W(1)+1.   
 
432
        S(1)=S(1)+Y(N)  
 
433
20      A(1)=A(1)+X(N)
 
434
        S(1)=S(1)/W(1)
 
435
        A(1)=A(1)/W(1)    
 
436
        I=1   
 
437
        XNU=M-1
 
438
30      IFF=1  
 
439
40      I1=I   
 
440
        IF(IMAX .GT. I) I = I+1    
 
441
        CHISQ=0.  
 
442
        DO 70 N=1,M    
 
443
        P(2)=X(N)-A(1)    
 
444
        IF(I .EQ. 2) GOTO 60 
 
445
        DO 50 K=3,I
 
446
        K1=K-1    
 
447
50      P(K)=(X(N)-A(K1))*P(K1)-B(K1)*P(K-2)  
 
448
60      DEL=Y(N)-S(I1)*P(I1)   
 
449
        Y(N)=DEL  
 
450
        CHISQ=CHISQ+DEL*DEL 
 
451
        IF(IMAX .LE. I1) GOTO 70  
 
452
        P2=P(I)**2
 
453
        S(I)=S(I)+DEL*P(I)                                                     
 
454
        A(I)=A(I)+X(N)*P2  
 
455
        W(I)=W(I)+P2    
 
456
70      CONTINUE       
 
457
        IF(IMAX .LE. I1) GOTO 80  
 
458
        A(I)=A(I)/W(I)  
 
459
        B(I)=W(I)/W(I1)  
 
460
        S(I)=S(I)/W(I)      
 
461
        XNU=XNU-1.     
 
462
        DCHI=S(I)**2*W(I)  
 
463
        IF(DCHI .GE. CHISQ) GOTO 30   
 
464
        F=XNU*DCHI/(CHISQ-DCHI)   
 
465
        F95=3.84+(10.+(12.+(30.+105./XNU/XNU)/XNU)/XNU)/XNU    
 
466
        IF(F .GT. F95) GOTO 30  
 
467
        IF(IFAUTO.EQ.0) GOTO 30  ! F-TEST OFF
 
468
        XNU=XNU+1.   
 
469
        IFF=IFF+1      
 
470
        S(I)=0.   
 
471
        IF(IFF .LE. IFTEST) GOTO 40  
 
472
80      IORD=MIN0(IMAX-1,I1)-IFF+1
 
473
        RETURN 
 
474
        END