← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/daophot/libsrc/addstar.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:
contrib/daophot/libsrc/addstar.for
 
1
C===========================================================================
 
2
C Copyright (C) 1995,2004 European Southern Observatory (ESO)
 
3
C
 
4
C This program is free software; you can redistribute it and/or 
 
5
C modify it under the terms of the GNU General Public License as 
 
6
C published by the Free Software Foundation; either version 2 of 
 
7
C the License, or (at your option) any later version.
 
8
C
 
9
C This program is distributed in the hope that it will be useful,
 
10
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
C GNU General Public License for more details.
 
13
C
 
14
C You should have received a copy of the GNU General Public 
 
15
C License along with this program; if not, write to the Free 
 
16
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
C MA 02139, USA.
 
18
C
 
19
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
C       Internet e-mail: midas@eso.org
 
21
C       Postal address: European Southern Observatory
 
22
C                       Data Management Division 
 
23
C                       Karl-Schwarzschild-Strasse 2
 
24
C                       D 85748 Garching bei Muenchen 
 
25
C                       GERMANY
 
26
C===========================================================================
 
27
C
 
28
      SUBROUTINE  ADDSTR  (PAR, MAXPAR, PSF, MAXPSF, MAXEXP,
 
29
     .     F, NCOL, NROW, WATCH)
 
30
C
 
31
C=======================================================================
 
32
C
 
33
C This subroutine generates random x,y coordinates and magnitudes,
 
34
C appropriately scales the point-spread function, and adds these
 
35
C artificial stars into a copy of the original picture at the
 
36
C appropriate locations.  As an alternative, it will read positions
 
37
C and magnitudes in from a disk file, and add synthetic stars with
 
38
C simulated photon noise.
 
39
C
 
40
C             OFFICIAL DAO VERSION:  1991 April 18
 
41
C
 
42
C.VERSIONS
 
43
C 040721        last modif
 
44
C 
 
45
C=======================================================================
 
46
C
 
47
      IMPLICIT NONE
 
48
      INTEGER MAXPSF, MAXPAR, MAXEXP, NCOL, NROW
 
49
C
 
50
C Parameters
 
51
C
 
52
C MAXPSF is the largest permissible number of elements in the look-up
 
53
C        table for the point-spread function.
 
54
C
 
55
      REAL F(NCOL,NROW), PSF(MAXPSF,MAXPSF,MAXEXP)
 
56
      REAL PAR(MAXPAR), RMAG(2)
 
57
      INTEGER ISEED(3)
 
58
C
 
59
      REAL SQRT, AMAX1, DAORAN, USEPSF
 
60
      INTEGER RDPSF, MIN0, MAX0, NINT
 
61
C
 
62
      CHARACTER*80 LINE
 
63
      CHARACTER*30 ADDPIC, OUTSTM, ADDFIL, EXTEND
 
64
      CHARACTER*30 COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL, SWITCH
 
65
      CHARACTER CASE*5, NTOCHR*2
 
66
      DOUBLE PRECISION SUMPHOT, SUMERR
 
67
      REAL LOBAD, NRML, STRMAG, X, Y, DIFMAX, DIFSQ, DX, DY, DYSQ
 
68
      REAL DIFF, SCALE, DELTAX, DELTAY, DVDXC, DVDYC, ERR
 
69
      REAL HIBAD, THRESH, AP1, DUMMY, READNS, FRAD, XWIDE, YWIDE
 
70
      REAL PSFRSQ, WATCH, PSFMAG, BRIGHT, XPSF, YPSF, PSFRAD, PHPADU
 
71
      REAL DATUM, SKY
 
72
      INTEGER I, J, ID, ISTAR, IDUM, IFRAME, LX, LY, NX, NY, ISTAT
 
73
      INTEGER IPSTYP, NPSF, NPAR, NEXP, NFRAC, NSTAR
 
74
      INTEGER NFRAME, NL
 
75
      LOGICAL RANDOM
 
76
C
 
77
      COMMON /FILNAM/ COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL
 
78
C
 
79
C-----------------------------------------------------------------------
 
80
C
 
81
      CALL TBLANK                                   ! Type a blank line
 
82
  950 CALL GETNAM ('File with the PSF:', PSFFIL)
 
83
      IF ((PSFFIL .EQ. 'END OF FILE') .OR.
 
84
     .     (PSFFIL .EQ. 'GIVE UP')) THEN
 
85
         PSFFIL = ' '
 
86
         RETURN
 
87
      END IF
 
88
      PSFFIL=EXTEND(PSFFIL, CASE('psf'))
 
89
C
 
90
C Read in the PSF.
 
91
C
 
92
      ISTAT = RDPSF(PSFFIL, IPSTYP, PAR, MAXPAR, NPAR,
 
93
     .     PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC,
 
94
     .     PSFMAG, BRIGHT, XPSF, YPSF)
 
95
      IF (ISTAT .NE. 0) THEN
 
96
         PSFFIL = 'GIVE UP'
 
97
         GO TO 950
 
98
      END IF
 
99
C
 
100
      PSFRAD = (REAL(NPSF-1)/2. - 1.)/2.
 
101
      PSFRSQ = PSFRAD**2
 
102
      CALL SEED3 (ISEED)
 
103
      CALL GETDAT ('Photons per ADU:', PHPADU,1)
 
104
      IF (PHPADU .LE. 0) RETURN
 
105
      ADDFIL='RANDOM STARS'
 
106
 1015 CALL GETNAM ('Input data file:', ADDFIL)
 
107
      IF ((ADDFIL .EQ. 'END OF FILE') .OR. (ADDFIL .EQ. 'EXIT') .OR.
 
108
     .     (ADDFIL .EQ. 'GIVE UP')) THEN
 
109
         RETURN
 
110
      ELSE IF (ADDFIL .EQ. 'RANDOM STARS') THEN
 
111
         RANDOM=.TRUE.
 
112
         WRITE (6,610) PSFMAG
 
113
  610    FORMAT (/' Magnitude of PSF star is', F7.3/)
 
114
         CALL GETDAT ('Minimum, maximum magnitudes desired:', RMAG, 2)
 
115
         IF (RMAG(1) .LT. -1.E38) RETURN            ! CTRL-Z was entered
 
116
C
 
117
         CALL GETDAT ('Number of stars to add to each frame:',
 
118
     .        DATUM, 1)
 
119
         IF (DATUM .LT. -1.E38) RETURN             ! CTRL-Z was entered
 
120
         NSTAR=MAX0(1, NINT(DATUM))
 
121
C
 
122
         CALL GETDAT ('Number of new frames:', DATUM, 1)
 
123
         IF (DATUM .LT. -1.E38) RETURN              ! CTRL-Z was entered
 
124
         NFRAME=MAX0(1, MIN0(NINT(DATUM), 99))
 
125
         OUTSTM=' '
 
126
         CALL GETNAM ('File-name stem:', OUTSTM)
 
127
         IF (OUTSTM .EQ. 'END OF FILE') RETURN     ! CTRL-Z was entered
 
128
         OUTSTM=EXTEND(OUTSTM, CASE('add'))
 
129
C
 
130
      ELSE
 
131
         ADDFIL = EXTEND(ADDFIL, CASE('add'))
 
132
         CALL INFILE (2, ADDFIL, ISTAT)
 
133
         IF (ISTAT .NE. 0) THEN
 
134
            CALL STUPID ('Error opening input file '//ADDFIL)
 
135
            ADDFIL = 'GIVE UP'
 
136
            GO TO 1015
 
137
         END IF
 
138
         ADDPIC=SWITCH(ADDFIL, ' ')
 
139
         CALL GETNAM ('Output picture name:', ADDPIC)
 
140
         NL=-1
 
141
         CALL RDHEAD (2, NL, IDUM, IDUM, LOBAD, HIBAD, THRESH, AP1,
 
142
     .        DUMMY, READNS, FRAD)
 
143
         IF (NL .LE. 0) NL=1
 
144
      END IF
 
145
C
 
146
C-----------------------------------------------------------------------
 
147
C
 
148
C SECTION 2
 
149
C
 
150
C Do it.
 
151
C
 
152
      IF (RANDOM) THEN
 
153
         IF (WATCH .GT. 0.5) THEN
 
154
            CALL TBLANK
 
155
            CALL OVRWRT (
 
156
     .      '  Star       Picture                         Data file', 1)
 
157
         END IF
 
158
C
 
159
C Beginning of loop over output frames.
 
160
C
 
161
         DO 2900 IFRAME=1,NFRAME
 
162
C
 
163
C Build up output stellar-data filename, and open file.
 
164
C
 
165
         ADDFIL=NTOCHR(IFRAME)//CASE('.add')
 
166
         ADDFIL=SWITCH(OUTSTM,ADDFIL)
 
167
  952    CALL OUTFIL (3, ADDFIL, ISTAT)
 
168
         IF (ISTAT .NE. 0) THEN
 
169
            CALL STUPID ('Error opening output file '//ADDFIL)
 
170
            ADDFIL = 'GIVE UP'
 
171
            CALL GETNAM ('New output file name:', ADDFIL)
 
172
            IF ((ADDFIL .EQ. 'END OF FILE') .OR.
 
173
     .           (ADDFIL .EQ. 'GIVE UP')) RETURN
 
174
            ADDFIL = EXTEND(ADDFIL, CASE('add'))
 
175
            GO TO 952
 
176
         END IF
 
177
         CALL WRHEAD (3, 1, NCOL, NROW, 7, 0., 0., 0., 0., 0., 0., 0.)
 
178
C
 
179
C Create output picture filename.
 
180
C
 
181
         ADDPIC=SWITCH(ADDFIL, ' ')
 
182
C
 
183
C Copy the input picture verbatim into the output picture.
 
184
C
 
185
         CALL COPPIC (ADDPIC, F, NCOL, NROW, ISTAT)
 
186
         IF (ISTAT .NE. 0) THEN
 
187
            CALL STUPID ('Error creating output picture.')
 
188
            RETURN
 
189
         END IF
 
190
         LX = 1
 
191
         LY = 1
 
192
         NX = NCOL
 
193
         NY = NROW
 
194
         CALL RDARAY ('COPY', LX, LY, NX, NY, NCOL, F, ISTAT)
 
195
C
 
196
C Beginning of loop over artificial stars.
 
197
C
 
198
         XWIDE = REAL(NCOL)-1.
 
199
         YWIDE = REAL(NROW)-1.
 
200
         DO 2500 ISTAR=1,NSTAR
 
201
C
 
202
C Make up centroid and magnitude.
 
203
C
 
204
         X=1.+XWIDE*DAORAN(ISEED(1))
 
205
         Y=1.+YWIDE*DAORAN(ISEED(2))
 
206
         STRMAG=RMAG(1)+DAORAN(ISEED(3))*(RMAG(2)-RMAG(1))
 
207
C
 
208
C Write them to the data file.
 
209
C
 
210
         WRITE (3,320) ISTAR+8999, X, Y, STRMAG
 
211
  320    FORMAT (1X, I5, 14F9.3)
 
212
         IF (WATCH .GT. 0.5) THEN
 
213
            WRITE (LINE,622) ISTAR, ADDPIC, ADDFIL
 
214
  622       FORMAT (I6, 7X, A30, 2X, A30)
 
215
            CALL OVRWRT (LINE(1:75), 2)
 
216
         END IF
 
217
         DELTAX=(X-1.)/XPSF - 1.
 
218
         DELTAY=(Y-1.)/YPSF - 1.
 
219
         LX = MAX0( 1, INT(X-PSFRAD)+1 )
 
220
         LY = MAX0( 1, INT(Y-PSFRAD)+1 )
 
221
         NX = MIN0( NCOL, INT(X+PSFRAD) )
 
222
         NY = MIN0( NROW, INT(Y+PSFRAD) )
 
223
         SCALE=10.**(0.4*(PSFMAG-STRMAG))
 
224
         DIFMAX=0.
 
225
         DIFSQ=0.
 
226
C
 
227
C Beginning of double loop over pixels in the subarray.
 
228
C
 
229
         SUMPHOT = 0.0D0
 
230
         SUMERR = 0.0D0
 
231
         DO 2210 J=LY,NY
 
232
            DY=FLOAT(J)-Y
 
233
            DYSQ=DY**2
 
234
C
 
235
            DO 2200 I=LX,NX
 
236
               DX=FLOAT(I)-X
 
237
               IF (DX**2+DYSQ .GE. PSFRSQ) THEN
 
238
                  IF (DX .GT. 0.) GO TO 2210
 
239
               ELSE
 
240
                  DIFF=SCALE*USEPSF(IPSTYP, DX, DY, BRIGHT, PAR, PSF,
 
241
     .                 NPSF, NPAR, NEXP, NFRAC, DELTAX, DELTAY, DVDXC,
 
242
     .                 DVDYC)
 
243
C
 
244
C DIFF represents the value of the stellar profile at this pixel.
 
245
C Compute a Poisson random error using a normal approximation,
 
246
C sigma(DIFF)=sqrt(DIFF/PHPADU).
 
247
C
 
248
                  ERR=SQRT(AMAX1(0.,DIFF/PHPADU))*
 
249
     .                 NRML(DAORAN(ISEED(MOD(I+J,3)+1)))
 
250
                  SUMPHOT = SUMPHOT + DBLE(DIFF)
 
251
                  SUMERR = SUMERR + DBLE(ERR)
 
252
                  F(I,J)=F(I,J)+DIFF+ERR
 
253
               END IF
 
254
 2200       CONTINUE
 
255
 2210    CONTINUE
 
256
 2500    CONTINUE                              ! End of loop over stars
 
257
         LX = 1
 
258
         LY = 1
 
259
         NX = NCOL
 
260
         NY = NROW
 
261
         CALL WRARAY ('COPY', LX, LY, NX, NY, NCOL, F, ISTAT)
 
262
         CALL CLPIC ('COPY')                   ! Close copy
 
263
         CALL CLFILE (3)
 
264
         CALL OVRWRT (' ', 1)
 
265
C
 
266
 2900    CONTINUE                              ! End of loop over frames
 
267
      ELSE
 
268
         CALL COPPIC (ADDPIC, F, NCOL, NROW, ISTAT)
 
269
         IF (ISTAT .NE. 0) THEN
 
270
            CALL STUPID ('Error creating output picture.')
 
271
            RETURN
 
272
         END IF
 
273
         IF (WATCH .GT. 0.5) THEN
 
274
            CALL TBLANK
 
275
            CALL OVRWRT ('  Star', 1)
 
276
         END IF
 
277
         LX = 1
 
278
         LY = 1
 
279
         NX = NCOL
 
280
         NY = NROW
 
281
         CALL RDARAY ('COPY', LX, LY, NX, NY, NCOL, F, ISTAT)
 
282
C
 
283
C Loop over stars.
 
284
C
 
285
         ISTAR=0
 
286
 3000    ISTAR=ISTAR+1
 
287
 3010    CALL RDSTAR (2, NL, ID, X, Y, STRMAG, SKY)
 
288
         IF (ID .LT. 0) GO TO 3900             ! End-of-file encountered
 
289
         IF (ID .EQ. 0) GO TO 3010             ! Ignore a blank line
 
290
         IF (STRMAG .GE. 99.) GO TO 3000       ! Ignore a bad star
 
291
         IF (WATCH .GT. 0.5) THEN
 
292
            WRITE (LINE,622) ISTAR
 
293
            CALL OVRWRT (LINE(1:6), 2)
 
294
         END IF
 
295
         DELTAX=(X-1.)/XPSF-1.
 
296
         DELTAY=(Y-1.)/YPSF-1.
 
297
         SCALE=10.**(0.4*(PSFMAG-STRMAG))
 
298
C
 
299
C Add the shifted scaled PSF
 
300
C
 
301
         LX = MAX0(1, INT(X-PSFRAD)+1 )
 
302
         LY = MAX0(1, INT(Y-PSFRAD)+1 )
 
303
         NX = MIN0(NCOL, INT(X+PSFRAD) )
 
304
         NY = MIN0(NROW, INT(Y+PSFRAD) )
 
305
         SUMPHOT = 0.0D0
 
306
         SUMERR = 0.0D0
 
307
         DO 3030 J=LY,NY
 
308
            DY=FLOAT(J)-Y
 
309
            DYSQ=DY**2
 
310
            DO 3020 I=LX,NX
 
311
               DX=FLOAT(I)-X
 
312
               IF (DX**2+DYSQ .GE. PSFRSQ) THEN
 
313
                  IF (DX .GT. 0.) GO TO 3030
 
314
               ELSE
 
315
                  DIFF=SCALE*USEPSF(IPSTYP, DX, DY, BRIGHT, PAR, PSF,
 
316
     .                 NPSF, NPAR, NEXP, NFRAC, DELTAX, DELTAY, DVDXC,
 
317
     .                 DVDYC)
 
318
C
 
319
C DIFF represents the value of the stellar profile at this pixel.
 
320
C Compute a Poisson random error using a normal approximation,
 
321
C sigma(DIFF)=sqrt(DIFF/PHPADU).
 
322
C
 
323
                  ERR=DAORAN(ISEED(MOD(I+J,3)+1))
 
324
                  ERR=SQRT(AMAX1(0.,DIFF/PHPADU))*NRML(ERR)
 
325
                  SUMPHOT = SUMPHOT + DBLE(DIFF)
 
326
                  SUMERR = SUMERR + DBLE(ERR)
 
327
                  F(I,J)=F(I,J)+DIFF+ERR
 
328
               END IF
 
329
 3020       CONTINUE
 
330
 3030    CONTINUE
 
331
         GO TO 3000                           ! Go to next star
 
332
C
 
333
 3900    LX = 1
 
334
         LY = 1
 
335
         NX = NCOL
 
336
         NY = NROW
 
337
         CALL WRARAY ('COPY', LX, LY, NX, NY, NCOL, F, ISTAT)
 
338
         CALL CLPIC ('COPY')
 
339
         CALL CLFILE (2)
 
340
         ADDFIL='EXIT'
 
341
         CALL OVRWRT (' ', 1)
 
342
         GO TO 1015
 
343
      END IF
 
344
C
 
345
      RETURN
 
346
C
 
347
      END!
 
348
C
 
349
C#######################################################################
 
350
C
 
351
      CHARACTER*2  FUNCTION  NTOCHR (N)
 
352
C
 
353
C Converts an integer in the range 1-99 to two characters representing
 
354
C the number.
 
355
C
 
356
      IMPLICIT NONE
 
357
      INTEGER N, ITENS
 
358
      IF ((N .GT. 0) .AND. (N .LT. 100)) GO TO 1010
 
359
      NTOCHR='00'
 
360
      RETURN
 
361
C
 
362
 1010 ITENS=N/10
 
363
      NTOCHR=CHAR(48+ITENS)//CHAR(48+N-10*ITENS)
 
364
      RETURN
 
365
C
 
366
      END!
 
367
C
 
368
C###################################################################
 
369
C
 
370
      SUBROUTINE  SEED3 (ISEED)
 
371
C
 
372
C Seed the random number generator, based on a floating-point number
 
373
C entered by the user.
 
374
C
 
375
      IMPLICIT NONE
 
376
      INTEGER ISEED(3)
 
377
C
 
378
      REAL DAORAN
 
379
      INTEGER INT
 
380
C
 
381
      REAL SEED
 
382
      INTEGER I
 
383
C
 
384
      CALL GETDAT ('Seed (any integer):', SEED, 1)
 
385
      I = NINT(SEED)
 
386
      ISEED(1) = INT(524288.*DAORAN(I)) + 1
 
387
      ISEED(2) = INT(524288.*DAORAN(I)) + 1
 
388
      ISEED(3) = INT(524288.*DAORAN(I)) + 1
 
389
      RETURN
 
390
C
 
391
      END!
 
392
C
 
393
C######################################################################
 
394
C
 
395
      REAL  FUNCTION  NRML  (RANNUM)
 
396
C
 
397
C Convert a uniform probability distribution to a Gaussian distribution
 
398
C with mean zero and standard deviation unity.
 
399
C
 
400
      IMPLICIT NONE
 
401
C
 
402
      REAL SQRT, ALOG
 
403
C
 
404
      REAL P, RANNUM, SIGN, T
 
405
 1000 P=RANNUM
 
406
      SIGN=-1.
 
407
      IF (P .GT. 0.5) THEN
 
408
         P=P-0.5
 
409
         SIGN=1.
 
410
      ELSE IF (P .LE. 0.) THEN
 
411
         NRML = -1.E20
 
412
         GO TO 1000
 
413
      END IF
 
414
      T=SQRT(ALOG(1/P**2))
 
415
      NRML=SIGN*(T- (2.30753+.27061*T) / (1.+T*(.99229+T*.04481)) )
 
416
      RETURN
 
417
      END!
 
418
C
 
419
C#######################################################################
 
420
C
 
421
      REAL FUNCTION DAORAN (IDUM)
 
422
      IMPLICIT NONE
 
423
C
 
424
C RAN2 from Numerical Recipes
 
425
C
 
426
      INTEGER IDUM
 
427
      INTEGER IFF
 
428
      INTEGER IY
 
429
      INTEGER M
 
430
      INTEGER IA, IC
 
431
      INTEGER J
 
432
      REAL    RM
 
433
      INTEGER IR(97)
 
434
C 
 
435
      SAVE    IY
 
436
      SAVE    IFF
 
437
      SAVE    IR
 
438
C
 
439
      DATA    IFF /0/
 
440
      DATA    M/714025/
 
441
      DATA    IA/1366/
 
442
      DATA    IC/150889/
 
443
      DATA    RM/1.400511E-6/
 
444
C
 
445
1000  CONTINUE
 
446
      IF ((IDUM .LT. 0) .OR. (IFF .EQ. 0)) THEN
 
447
         IFF = 1
 
448
         IDUM = MOD(IABS(IC-IDUM), M)
 
449
         DO J=1,97
 
450
            IDUM = MOD(IA*IDUM+IC,M)
 
451
            IR(J) = IDUM
 
452
         END DO
 
453
         IDUM = MOD(IA*IDUM+IC,M)
 
454
         IY = IDUM
 
455
      END IF
 
456
C 
 
457
      J = 1+(97*IY)/M
 
458
      IF ((J .GT. 97) .OR. (J .LT. 1)) 
 
459
     +   CALL STETER(33,'Problems in random number generator DAORAN')
 
460
 
 
461
CC    was:   IF ((J .GT. 97) .OR. (J .LT. 1)) PAUSE
 
462
 
 
463
      IY = IR(J)
 
464
      DAORAN = IY*RM
 
465
      IDUM = MOD(IA*IDUM+IC,M)
 
466
      IR(J) = IDUM
 
467
      IF (DAORAN .LE. 0.) GO TO 1000             ! Stetson's modification
 
468
C
 
469
      RETURN
 
470
      END