← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/daophot/libsrc/substar.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/substar.for
 
1
C @(#)substar.for       19.1 (ES0-DMD) 02/25/03 13:23:51
 
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
      SUBROUTINE  SUBSTR  (PAR, MAXPAR, PSF, MAXPSF, MAXEXP, 
 
30
     .     F, NCOL, NROW, WATCH)
 
31
C
 
32
C=======================================================================
 
33
C
 
34
C This subroutine scales and shifts the point-spread function according
 
35
C to each star's magnitude and centroid, and subtracts the resulting
 
36
C profile from a copy of the original picture.
 
37
C
 
38
C             OFFICIAL DAO VERSION:  1991 April 18
 
39
C
 
40
C Arguments
 
41
C
 
42
C  WATCH (INPUT) governs whether information relating to the progress 
 
43
C        of the reductions is to be typed on the terminal screen
 
44
C        during execution.
 
45
C
 
46
C=======================================================================
 
47
C
 
48
      IMPLICIT NONE
 
49
      INTEGER MAXBOX, MAXPSF, MAXPAR, MAXEXC, MAXEXP
 
50
      PARAMETER  (MAXBOX=69, MAXEXC=200)
 
51
C
 
52
C Parameters
 
53
C
 
54
C MAXBOX is the side of the square subarray containing the largest
 
55
C        (circular) PSF that can be subtracted from the picture.
 
56
C
 
57
C MAXPSF is the largest permissible number of elements on a side of the
 
58
C        (square) look-up table for the point-spread function.
 
59
C
 
60
C        MAXBOX = (MAXPSF-7)/2.
 
61
C
 
62
      INTEGER NCOL, NROW
 
63
      REAL PAR(MAXPAR)
 
64
      REAL F(NCOL,NROW), PSF(MAXPSF,MAXPSF,MAXEXP)
 
65
      INTEGER IEXC(MAXEXC)
 
66
C
 
67
      REAL USEPSF
 
68
      INTEGER RDPSF, MAX0, MIN0, INT
 
69
C
 
70
      CHARACTER*80 LINE
 
71
      CHARACTER*30 SUBPIC
 
72
      CHARACTER*30 COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL, SWITCH
 
73
      CHARACTER CASE*5, ANSWER*1
 
74
      REAL LOBAD, X, Y, STRMAG, SKY, DX, DY, DYSQ, SCALE, COL, ROW
 
75
      REAL DELTAX, DELTAY, DVDXC, DVDYC, XPSF, YPSF, PSFRAD, DUM
 
76
      REAL HIBAD, THRESH, AP1, PHPADU, READNS, FRAD, WATCH, PSFMAG
 
77
      REAL BRIGHT, PSFRSQ
 
78
      INTEGER I, J, L, NL, IDUM, NEXC, LX, LY, NX, NY, ISTAR, ID
 
79
      INTEGER ISTAT, IPSTYP, NPSF, NPAR, NEXP, NFRAC
 
80
C
 
81
      COMMON /FILNAM/ COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL
 
82
C
 
83
C-----------------------------------------------------------------------
 
84
C
 
85
C SECTION 1
 
86
C
 
87
C Get file names and set up the needed numerical constants.
 
88
C
 
89
      CALL TBLANK                                   ! Type a blank line
 
90
  950 CALL GETNAM ('File with the PSF:', PSFFIL)
 
91
      IF ((PSFFIL .EQ. 'END OF FILE') .OR.
 
92
     .     (PSFFIL .EQ. 'GIVE UP')) THEN
 
93
         PSFFIL = ' '
 
94
         RETURN
 
95
      END IF
 
96
      ISTAT = RDPSF (PSFFIL, IPSTYP, PAR, MAXPAR, NPAR,
 
97
     .     PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC, 
 
98
     .     PSFMAG, BRIGHT, XPSF, YPSF)
 
99
      IF (ISTAT .LT. 0) THEN
 
100
         PSFFIL = 'GIVE UP'
 
101
         GO TO 950
 
102
      END IF
 
103
C
 
104
      PSFRAD=(REAL(NPSF-1)/2. - 1.)/2.
 
105
      PSFRSQ=PSFRAD**2
 
106
      COL=FLOAT(NCOL)
 
107
      ROW=FLOAT(NROW)
 
108
C
 
109
      CALL CLFILE (2)
 
110
C
 
111
 1015 CALL GETNAM ('File with photometry:', PROFIL)
 
112
      IF ((PROFIL .EQ. 'END OF FILE') .OR.
 
113
     .     (PROFIL .EQ. 'GIVE UP')) THEN
 
114
         PROFIL = ' '
 
115
         RETURN
 
116
      END IF
 
117
      CALL INFILE (2, PROFIL, ISTAT)
 
118
      IF (ISTAT .NE. 0) THEN
 
119
         CALL STUPID ('Error opening input file '//PROFIL)
 
120
         PROFIL = 'GIVE UP'
 
121
         GO TO 1015
 
122
      END IF
 
123
C
 
124
      CALL GETYN ('Do you have stars to leave in?', ANSWER)
 
125
      IF (ANSWER .EQ. 'Y') THEN
 
126
         SUBPIC = SWITCH(PROFIL, CASE('.lst'))
 
127
 1025    CALL GETNAM ('File with star list:', SUBPIC)
 
128
         IF ((SUBPIC .EQ. 'END OF FILE') .OR. (SUBPIC .EQ.
 
129
     .        'GIVE UP')) THEN
 
130
            CALL CLFILE (2)
 
131
            RETURN
 
132
         END IF
 
133
         CALL INFILE (3, SUBPIC, ISTAT)
 
134
         IF (ISTAT .NE. 0) THEN
 
135
            CALL STUPID ('Error opening file '//SUBPIC)
 
136
            SUBPIC = 'GIVE UP'
 
137
            GO TO 1025
 
138
         END IF
 
139
         NL = -1
 
140
         CALL RDHEAD (3, NL, IDUM, IDUM, DUM, DUM, DUM, DUM,
 
141
     .        DUM, DUM, DUM)
 
142
         NEXC = 0
 
143
 1035    NEXC = NEXC+1
 
144
 1036    CALL RDSTAR (3, NL, IEXC(NEXC), DUM, DUM, DUM, DUM)
 
145
         IF (IEXC(NEXC) .EQ. 0) GO TO 1036
 
146
         IF (IEXC(NEXC) .LT. 0) THEN
 
147
            NEXC = NEXC-1
 
148
         ELSE
 
149
            IF (NEXC .LT. MAXEXC) GO TO 1035
 
150
         END IF
 
151
         IF (NEXC .LE. 0) ANSWER = 'N'
 
152
         CALL CLFILE (3)
 
153
      END IF
 
154
C
 
155
      CALL RDHEAD (2, NL, IDUM, IDUM, LOBAD, HIBAD, THRESH, AP1, 
 
156
     .     PHPADU, READNS, FRAD)
 
157
C
 
158
      SUBPIC=SWITCH(PROFIL, CASE('s'))
 
159
      CALL GETNAM ('Name for subtracted image:', SUBPIC)
 
160
      IF (SUBPIC .EQ. 'END OF FILE') GO TO 9010    ! CTRL-Z was entered
 
161
C
 
162
C Copy the input picture verbatim into the output picture.
 
163
C
 
164
      CALL COPPIC (SUBPIC, F, NCOL, NROW, ISTAT)
 
165
      IF (ISTAT .NE. 0) THEN
 
166
         CALL STUPID ('Error creating output image.')
 
167
         CALL CLFILE (2)
 
168
         RETURN
 
169
      END IF
 
170
C
 
171
      IF (WATCH .GT. 0.5) THEN
 
172
         CALL TBLANK
 
173
         CALL OVRWRT (' Star', 1)
 
174
      END IF
 
175
C
 
176
C Read the entire image in from the copy.
 
177
C
 
178
      LX = 1
 
179
      LY = 1
 
180
      NX = NCOL
 
181
      NY = NROW
 
182
      CALL RDARAY ('COPY', LX, LY, NX, NY, NCOL, F, ISTAT)
 
183
C
 
184
C-----------------------------------------------------------------------
 
185
C
 
186
C SECTION 2
 
187
C
 
188
C Loop over stars.
 
189
C
 
190
      ISTAR=0
 
191
 2000 ISTAR=ISTAR+1
 
192
 2010 CALL RDSTAR (2, NL, ID, X, Y, STRMAG, SKY)
 
193
      IF (ID .LT. 0) GO TO 9000                ! End-of-file encountered
 
194
      IF (ID .EQ. 0) GO TO 2010                ! Ignore a blank line
 
195
      IF ((1.-X .GT. PSFRAD) .OR. (X-COL .GT. PSFRAD) .OR.
 
196
     .    (1.-Y .GT. PSFRAD) .OR. (Y-ROW .GT. PSFRAD)) GO TO 2000
 
197
      IF (ANSWER .EQ. 'Y') THEN
 
198
         L = NEXC
 
199
         DO I=1,L
 
200
            IF (ID .EQ. IEXC(I)) THEN
 
201
               IEXC(I) = IEXC(NEXC)
 
202
               NEXC = NEXC-1
 
203
               IF (NEXC .EQ. 0) ANSWER = 'N'
 
204
               GO TO 2010
 
205
            END IF
 
206
         END DO
 
207
      END IF
 
208
      IF (STRMAG .GE. 99.) GO TO 2000          ! Ignore a bad star
 
209
      IF (WATCH .GT. 0.5) THEN
 
210
         WRITE (LINE,620) ISTAR
 
211
  620    FORMAT (I5)
 
212
         CALL OVRWRT (LINE(1:5), 2)
 
213
      END IF
 
214
      DELTAX=(X-1.)/XPSF-1.
 
215
      DELTAY=(Y-1.)/YPSF-1.
 
216
      LX=MAX0(1, INT(X-PSFRAD)+1)
 
217
      LY=MAX0(1, INT(Y-PSFRAD)+1)
 
218
      NX=MIN0(NCOL, INT(X+PSFRAD))
 
219
      NY=MIN0(NROW, INT(Y+PSFRAD))
 
220
      SCALE=10.**(0.4*(PSFMAG-STRMAG))
 
221
C
 
222
C Subtract the shifted scaled PSF
 
223
C
 
224
      DO 2030 J=LY,NY
 
225
         DY=FLOAT(J)-Y
 
226
         DYSQ=DY**2
 
227
         DO 2020 I=LX,NX
 
228
            IF (F(I,J) .GT. HIBAD) GO TO 2020
 
229
            DX=FLOAT(I)-X
 
230
            IF (DX**2+DYSQ .GE. PSFRSQ) THEN
 
231
               IF (DX .GT. 0.) GO TO 2030
 
232
            ELSE
 
233
               F(I,J)=F(I,J)-SCALE*USEPSF(IPSTYP, DX, DY, BRIGHT, PAR, 
 
234
     .              PSF, NPSF, NPAR, NEXP, NFRAC, DELTAX, DELTAY, 
 
235
     .              DVDXC, DVDYC)
 
236
            END IF
 
237
 2020    CONTINUE
 
238
 2030 CONTINUE
 
239
      GO TO 2000                              ! Go to next star
 
240
C
 
241
C-----------------------------------------------------------------------
 
242
C
 
243
C Normal return.
 
244
C
 
245
 9000 CONTINUE
 
246
      IF (WATCH .GT. 0.5) CALL OVRWRT (' ', 4)
 
247
C
 
248
C Write the modified image back into the copy.
 
249
C
 
250
      LX = 1
 
251
      LY = 1
 
252
      NX = NCOL
 
253
      NY = NROW
 
254
      CALL WRARAY ('COPY', LX, LY, NX, NY, NCOL, F, ISTAT)
 
255
      CALL CLPIC ('COPY')
 
256
      CALL STUPID ('    Done.  ')
 
257
 9010 CALL CLFILE (2)
 
258
      RETURN
 
259
C
 
260
      END!