← Back to branch summary

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

« back to all changes in this revision

Viewing changes to prim/general/src/rcosmic.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/rcosmic.for
 
1
C @(#)rcosmic.for       13.1.1.1 (ES0-DMD) 06/02/98 18:17:32
 
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 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
31
C                                         all rights reserved
 
32
C.IDENTIFICATION: RCOSMIC
 
33
C.PURPOSE:        Remove cosmic ray events on single ccd exposure and replace
 
34
C                 them by interpolation on neighbourhood pixels.
 
35
C.LANGUAGE:       F77+ESOext
 
36
C.AUTHOR:         P.Magain, M.Remy, Institut d'Astrophysique de Liege 
 
37
C.INPUT/OUTPUT:   IN_A/C/1/60      = Input frame.
 
38
C                 IN_B/C/1/60      = median filter of the input frame
 
39
C                 OUT_B/C/1/12     = Output frame.
 
40
C                 PARAMS/R/1/1     = Mean value of the background.
 
41
C                 PARAMS/R/2/1     = Readout noise in ADU units.
 
42
C                 PARAMS/R/3/1     = Inverse gain factor (e-/ADU).
 
43
C VERSION  1.2    New ST interfaces M.Peron 221191
 
44
C VERSION  1.3    remove limitation of size, suggestion from Hans Schwengeler
 
45
C ----------------------------------------------------------------------------
 
46
      PROGRAM RCOSMIC
 
47
      INTEGER*4 IAV,I
 
48
      INTEGER*4 STATUS,MADRID,SIZEX,IOMODE
 
49
      INTEGER*4 NAXIS,NPIX(2),IMNI,IMNO,IMNF,IMNC
 
50
      INTEGER*8 PNTRI,PNTRF,PNTRO,PNTRC
 
51
      INTEGER*4 KUN,KNUL
 
52
      CHARACTER*60 IMAGE,OBJET,COSMIC
 
53
      CHARACTER*72 IDENT1,IDENT2,IDENT3
 
54
      CHARACTER*48 CUNIT
 
55
      DOUBLE PRECISION START(2),STEP(2)
 
56
      REAL*4 SKY,GAIN,RON,NS,RC,PARAM(5),CUTS(2)
 
57
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
58
      COMMON/VMR/MADRID(1)
 
59
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
60
      DATA IDENT1 /' '/
 
61
      DATA IDENT2 /' '/
 
62
      DATA IDENT3 /'cosmic ray mask '/
 
63
      DATA CUNIT /' '/
 
64
      CALL STSPRO('RCOSMIC')
 
65
      CALL STKRDC('IN_A',1,1,60,IAV,IMAGE,KUN,KNUL,STATUS)
 
66
      CALL STIGET(IMAGE,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
67
     1                2,NAXIS,NPIX,START,STEP
 
68
     1                ,IDENT1,CUNIT,PNTRI,IMNI,STATUS)
 
69
 
 
70
      CALL STKRDR('PARAMS',1,5,IAV,PARAM,KUN,KNUL,STATUS)
 
71
      CALL STIGET('middumma',D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
72
     1                2,NAXIS,NPIX,START,STEP
 
73
     1                ,IDENT2,CUNIT,PNTRF,IMNF,STATUS)
 
74
      SKY = PARAM(1)
 
75
      GAIN = PARAM(2)
 
76
      RON = PARAM(3)
 
77
      NS = PARAM(4)
 
78
      RC = PARAM(5)
 
79
      CALL STKRDC('OUTIMA',1,1,60,IAV,OBJET,KUN,KNUL,STATUS)
 
80
      CALL STIPUT(OBJET,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
81
     1                 NAXIS,NPIX,START,STEP
 
82
     1                ,IDENT1,CUNIT,PNTRO,IMNO,STATUS)
 
83
 
 
84
      SIZEX = 1
 
85
      DO I=1,NAXIS
 
86
         SIZEX = SIZEX*NPIX(I)
 
87
      ENDDO
 
88
      CALL STKRDC('COSMIC',1,1,60,IAV,COSMIC,KUN,KNUL,STATUS)
 
89
      IF (COSMIC(1:1).EQ.'+') THEN
 
90
            COSMIC = 'dummy_frame'
 
91
            IOMODE = F_X_MODE
 
92
      ELSE
 
93
            IOMODE = F_O_MODE
 
94
      ENDIF    
 
95
      CALL STIPUT(COSMIC,D_I2_FORMAT,IOMODE,F_IMA_TYPE
 
96
     1                 ,NAXIS,NPIX,START,STEP
 
97
     1                ,IDENT3,CUNIT,PNTRC,IMNC,STATUS)
 
98
      CALL COSROUT(MADRID(PNTRI),MADRID(PNTRC),NPIX(1),NPIX(2),
 
99
     1             RON,GAIN,NS,SKY,RC
 
100
     1            ,MADRID(PNTRF),MADRID(PNTRO))
 
101
 
 
102
      CUTS(1) = 0
 
103
      CUTS(2) = 1
 
104
      IF (IOMODE.EQ.F_O_MODE) 
 
105
     + CALL STDWRR(IMNC,'LHCUTS',CUTS,1,2,KUN,STATUS)
 
106
      CALL DSCUPT(IMNI,IMNO,' ',STATUS) 
 
107
      CALL STSEPI
 
108
      END
 
109
      SUBROUTINE COSROUT(AI,COSMIC,I_IMA,J_IMA,RON,GAIN,
 
110
     1                   NS,SKY,RC,AM,AO)
 
111
      INTEGER I_IMA,J_IMA,NUM
 
112
      INTEGER ORD(10000)
 
113
      INTEGER K,L
 
114
      INTEGER IDUMAX,JDUMAX,I1,I2,J1,II,JJ
 
115
      INTEGER I,J,IMAX,JMAX,IMIN,JMIN
 
116
      INTEGER FIRST(2),NEXT(2)
 
117
      INTEGER*2 COSMIC(I_IMA,J_IMA)
 
118
      REAL*4 VECTEUR(10000),FMAX,ASUM,RC
 
119
      REAL*4 AI(I_IMA,J_IMA),AO(I_IMA,J_IMA),AM(I_IMA,J_IMA)
 
120
      REAL*4 SIGMA,SKY,S1,S2
 
121
      REAL*4 RON,GAIN,NS,AMEDIAN
 
122
      DO 10 J=1,J_IMA
 
123
      DO 5 I=1,I_IMA
 
124
      AO(I,J)=AI(I,J)
 
125
      COSMIC(I,J)= 0
 
126
    5 CONTINUE
 
127
   10 CONTINUE
 
128
C
 
129
C     La boucle suivante selectionne les pixels qui sont
 
130
C     significativement au dessus de l'image filtree medianement.
 
131
C
 
132
C
 
133
C     COSMIC =-1 ----> candidate for cosmic
 
134
C            = 0 ----> not a cosmic
 
135
C            = 1 -----> a cosmic (at the end)
 
136
C            = 2 ----> member of the group
 
137
C            = 3 ----> member of a group which has been examined
 
138
C            = 4 ----> neighbourhood  of the group
 
139
      K=1
 
140
      DO 80 J=2,J_IMA-1
 
141
      DO 70 I=2,I_IMA-1
 
142
      SIGMA=SQRT(RON**2+AM(I,J)/GAIN)
 
143
      IF ((AI(I,J)-AM(I,J)).GE.(NS*SIGMA)) THEN
 
144
            COSMIC(I,J) = -1
 
145
            K = K+1
 
146
      ENDIF
 
147
   70 CONTINUE
 
148
   80 CONTINUE
 
149
C
 
150
C ------------------------
 
151
C
 
152
C
 
153
C
 
154
C
 
155
C     Ces pixels sont regroupes par ensembles connexes dans la boucle
 
156
C
 
157
C
 
158
C  -----------------------------
 
159
      FIRST(1) = 2
 
160
      FIRST(2) = 2
 
161
  100 CALL FINDNEXT(COSMIC,I_IMA,J_IMA,FIRST,NEXT)
 
162
      IF (NEXT(1).EQ.-1) RETURN
 
163
      I = NEXT(1)
 
164
      J = NEXT(2) 
 
165
      COSMIC(I,J) = 2
 
166
      IMIN = I
 
167
      IMAX = I 
 
168
      JMIN = J
 
169
      JMAX = J
 
170
      IDUMAX = I
 
171
      JDUMAX = J
 
172
      FMAX = AI(I,J)
 
173
  110 I1 = 0
 
174
      I2 = 0
 
175
      CONTINUE
 
176
      DO 125 L=1,2
 
177
          DO 115 K=1,2
 
178
               II = I+K-L
 
179
               JJ = J+K+L-3
 
180
               IF (COSMIC(II,JJ).EQ.-1) THEN
 
181
                   I1 = II
 
182
                   J1 = JJ  
 
183
                   IMIN = MIN(IMIN,II) 
 
184
                   IMAX = MAX(IMAX,II)
 
185
                   JMIN = MIN(JMIN,JJ)
 
186
                   JMAX = MAX(JMAX,JJ)
 
187
                   COSMIC(II,JJ) = 2
 
188
                   IF (AI(II,JJ).GT.FMAX) THEN
 
189
                         FMAX = AI(II,JJ)
 
190
                         IDUMAX = II
 
191
                         JDUMAX = JJ
 
192
                   ENDIF
 
193
                ELSE IF (COSMIC(II,JJ).EQ.0) THEN
 
194
                   COSMIC(II,JJ) = 4
 
195
                ENDIF
 
196
  115     CONTINUE 
 
197
  125 CONTINUE 
 
198
      COSMIC(I,J) = 3
 
199
      IF (I1.NE.0) THEN
 
200
      I = I1
 
201
      J = J1
 
202
      GOTO 110
 
203
      ENDIF        
 
204
      DO 140 L = JMIN,JMAX  
 
205
         DO 130 K = IMIN,IMAX
 
206
              IF (COSMIC(K,L).EQ.2) THEN
 
207
                 I = K
 
208
                 J = L
 
209
                 GOTO 110
 
210
              ENDIF
 
211
  130 CONTINUE
 
212
  140 CONTINUE   
 
213
      FIRST(1) = NEXT(1)+1
 
214
      FIRST(2) = NEXT(2) 
 
215
C We start here the real work....
 
216
C 1- decide if the pixel's group is a cosmic
 
217
C 2-replace these values by another one
 
218
      
 
219
      S1 = AI(IDUMAX-1,JDUMAX-1) + AI(IDUMAX+1,JDUMAX-1)
 
220
     +     + AI(IDUMAX,JDUMAX-1) +AI(IDUMAX,JDUMAX+1)
 
221
      S2 = AI(IDUMAX-1,JDUMAX+1) + AI(IDUMAX+1,JDUMAX+1)
 
222
     +     + AI(IDUMAX-1,JDUMAX) + AI(IDUMAX+1,JDUMAX)
 
223
      ASUM = (S1+S2)/8.-SKY
 
224
      IF ((FMAX-SKY).GT.RC*ASUM) THEN
 
225
         NUM = 1
 
226
         DO L = JMIN-1,JMAX+1
 
227
            DO K = IMIN-1,IMAX+1
 
228
               IF (COSMIC(K,L).EQ.4) THEN
 
229
                   VECTEUR(NUM) = AI(K,L)
 
230
                   NUM = NUM+1
 
231
               ENDIF    
 
232
            ENDDO
 
233
         ENDDO
 
234
         CALL SORT(NUM-1,VECTEUR,ORD)
 
235
         AMEDIAN = VECTEUR(ORD((NUM-1)/2))
 
236
         DO L = JMIN-1,JMAX+1
 
237
            DO K = IMIN-1,IMAX+1
 
238
               IF (COSMIC(K,L).EQ.3) THEN
 
239
                   COSMIC(K,L) = 1
 
240
                   AO(K,L) = AMEDIAN
 
241
               ELSE IF (COSMIC(K,L).EQ.4) THEN
 
242
                   COSMIC(K,L) = 0
 
243
               ENDIF
 
244
            ENDDO
 
245
         ENDDO
 
246
      ELSE 
 
247
         DO L = JMIN-1,JMAX+1
 
248
            DO K = IMIN-1,IMAX+1
 
249
               IF (COSMIC(K,L).NE.-1) COSMIC(K,L) = 0
 
250
            ENDDO
 
251
          ENDDO
 
252
      ENDIF
 
253
        
 
254
      
 
255
 
 
256
      IF (NEXT(1).GT.0) GOTO 100
 
257
C
 
258
C
 
259
C
 
260
      RETURN
 
261
      END
 
262
 
 
263
      SUBROUTINE FINDNEXT(COSMIC,I_IMA,J_IMA,FIRST,NEXT)
 
264
      INTEGER I_IMA,J_IMA,FIRST(2),NEXT(2)
 
265
      INTEGER I,J
 
266
      INTEGER*2 COSMIC(I_IMA,J_IMA)
 
267
      DO J = FIRST(2), J_IMA
 
268
          DO I = 2, I_IMA
 
269
             IF (COSMIC(I,J).EQ.-1) THEN
 
270
                 NEXT(1) = I
 
271
                 NEXT(2) = J
 
272
                 RETURN
 
273
             ENDIF
 
274
          ENDDO
 
275
      ENDDO 
 
276
      NEXT(1) = -1
 
277
      NEXT(2) = -1
 
278
      RETURN
 
279
      END
 
280
 
 
281
      SUBROUTINE SORT(KMAX,INP,ORD)
 
282
      INTEGER KMAX,IMIN,IMAX,I,J,K,L
 
283
      INTEGER ORD(10000)
 
284
      REAL*4 INP(10000),F
 
285
      DO 4100 J=1,KMAX
 
286
      ORD(J)=J
 
287
 4100 CONTINUE
 
288
      IF (INP(1).GT.INP(2)) THEN 
 
289
             ORD(1)=2
 
290
             ORD(2)=1
 
291
      END IF
 
292
      DO 4400 J=3,KMAX
 
293
      F=INP(J)
 
294
      L=ORD(J-1)
 
295
      IF (INP(L).LE.F) GO TO 4400
 
296
      L=ORD(1)
 
297
      IMIN=1
 
298
      IF (F.LE.INP(L)) GO TO 4250
 
299
      IMAX=J-1
 
300
 4200 I=(IMIN+IMAX)/2
 
301
      L=ORD(I)
 
302
      IF (INP(L).LT.F) THEN
 
303
              IMIN=I
 
304
              ELSE
 
305
              IMAX=I
 
306
      END IF
 
307
      IF (IMAX.GT.(IMIN+1)) GO TO 4200
 
308
      IMIN=IMAX
 
309
 4250 DO 4300 K=J-1,IMIN,-1
 
310
      ORD(K+1)=ORD(K)
 
311
 4300 CONTINUE
 
312
      ORD(IMIN)=J
 
313
 4400 CONTINUE
 
314
      RETURN
 
315
      END