← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/invent/libsrc/twodim.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/invent/libsrc/twodim.for
 
1
C @(#)twodim.for        19.1 (ES0-DMD) 02/25/03 13:25:40
 
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
31
C
 
32
C
 
33
C
 
34
C-----------------------------------------------------------------------
 
35
      SUBROUTINE TWODIM(AS, MAS, JAPYS, IBUFS, LZ,
 
36
     &                  APSF, FPSF, IPSF, NCAT, PMTR,
 
37
     &                  PRCT, NCT, PMT, PRC, IARR,
 
38
     &                  RARR, SIGMA)
 
39
C
 
40
      IMPLICIT    NONE
 
41
      INCLUDE    'MID_REL_INCL:INVENT.INC/NOLIST'
 
42
C
 
43
      REAL        AS(1)
 
44
      INTEGER     MAS((-MAXSUB):MAXSUB,(-MAXSUB):MAXSUB)
 
45
      INTEGER     JAPYS(1)
 
46
      INTEGER     IBUFS(4)
 
47
      INTEGER     LZ
 
48
      REAL        APSF(0:MAXSUB)
 
49
      REAL        FPSF(1)
 
50
      INTEGER     IPSF(1)
 
51
      INTEGER     NCAT(NIPAR,MAXCNT)
 
52
      REAL        PMTR(NRPAR,MAXCNT)
 
53
      REAL        PRCT(0:MAXSUB,MAXCNT)
 
54
      INTEGER     NCT(NIPAR)
 
55
      REAL        PMT(NRPAR)
 
56
      REAL        PRC(0:MAXSUB)
 
57
      INTEGER     IARR(32)
 
58
      REAL        RARR(64)
 
59
      REAL        SIGMA
 
60
C
 
61
      INTEGER   NP
 
62
      INTEGER   LPXL
 
63
      INTEGER   LSBP
 
64
 
 
65
      INTEGER   IC
 
66
      INTEGER   ICM
 
67
      INTEGER   IDF
 
68
      INTEGER   ISC
 
69
      INTEGER   IOFF, IOFF1, IOFF0, IOFFI, IOFF2, IOFFC
 
70
      INTEGER   IS, ISOLD
 
71
      INTEGER   ISTART, IACT
 
72
      INTEGER   JCM
 
73
      INTEGER   JDF
 
74
      INTEGER   JSTART, JACT
 
75
      INTEGER   I0, I0OLD
 
76
      INTEGER   J0, J0OLD
 
77
      INTEGER   JSC
 
78
      INTEGER   JS, JSOLD
 
79
      INTEGER   KCLN
 
80
      INTEGER   KSAT
 
81
      INTEGER   LPLIM
 
82
      INTEGER   MSBP2, MSBP
 
83
      INTEGER   MPXL, MPXL2
 
84
      INTEGER   MPRF
 
85
      INTEGER   NNOR
 
86
      INTEGER   NPAS
 
87
      INTEGER   NPX, NPXOLD
 
88
C     REAL      TMP(5)
 
89
      REAL      BGRD
 
90
      REAL      CORX, CORY
 
91
      REAL      SCA, SCC, SCAA
 
92
      REAL      SIGM, SHLM, SIGM1
 
93
      REAL      XDF
 
94
      REAL      YDF
 
95
      DOUBLE PRECISION   AR(25) , X(5) , S(5)
 
96
      LOGICAL   CENTER , COMP , LAST , STAR
 
97
C     LOGICAL   JCNTR
 
98
 
 
99
      NP = MAXCNT
 
100
      LPXL = IARR(20)
 
101
      LSBP = IARR(21)
 
102
      BGRD = PMTR(1,LZ)
 
103
      KCLN = NCAT(4,LZ)
 
104
      KSAT = NCAT(6,LZ)
 
105
      IF ( NCAT(10,LZ) .GT. 0 ) THEN
 
106
          COMP = .TRUE.
 
107
      ELSE
 
108
          COMP = .FALSE.
 
109
      ENDIF
 
110
      IF ( KSAT .EQ. -1 .AND. ABS(PMTR(3,LZ)) .LT. RARR(5) ) THEN
 
111
          STAR = .TRUE.
 
112
      ELSE
 
113
          STAR = .FALSE.
 
114
      ENDIF
 
115
      MSBP = 2 * LSBP + 1
 
116
      MPXL = 2 * LPXL + 1
 
117
      MSBP2 = MSBP * MSBP
 
118
      MPXL2 = MPXL * MPXL
 
119
      IOFF = MSBP2 * MPXL2
 
120
      MPRF = NOSP + MSBP2
 
121
C
 
122
C *** Determine scaling parameter SCA for object and SCC for component.
 
123
C
 
124
      CALL SCLDET( LZ, NCAT , PMTR , PRCT , NCT ,
 
125
     &             PMT , PRC , APSF , NPAS , SCA ,
 
126
     &             SCC )
 
127
      IF ( SCA .LE. 5.0*RARR(3) ) THEN
 
128
          STAR = .FALSE.
 
129
      ENDIF
 
130
C
 
131
C *** Save starting object coordinates.
 
132
C
 
133
      ISTART = NCAT(1,LZ)
 
134
      JSTART = NCAT(2,LZ)
 
135
C
 
136
C *** Start loop.
 
137
C
 
138
      I0 = 0
 
139
      J0 = 0
 
140
      I0OLD = 0
 
141
      J0OLD = 0
 
142
      IS    = 0                            ! addition for initialization 
 
143
      JS    = 0                            ! addition for initialization 
 
144
      ISOLD = IS
 
145
      JSOLD = JS
 
146
      NPX = 0
 
147
      NPXOLD = -1
 
148
      IC = 1
 
149
      LAST = .FALSE.
 
150
      SIGM = BGRD
 
151
      SIGM1 = 0.0
 
152
C
 
153
   10 CONTINUE
 
154
          IF ( IC .GT. 6 ) GOTO 11
 
155
C
 
156
          IF ( IC .EQ. 1 ) THEN
 
157
              SHLM = 0.6
 
158
          ELSE
 
159
              SHLM = MIN( 0.6 , 1.0 / FLOAT(MSBP) )
 
160
          ENDIF
 
161
C
 
162
C ***     Identify object coordinates and subpixel.
 
163
C
 
164
          CALL SUBPXL( NCAT(1,LZ) , PMTR(1,LZ) , LSBP , IACT , JACT ,
 
165
     &                 IS , JS )
 
166
          I0 = IACT - ISTART
 
167
          J0 = JACT - JSTART
 
168
          IF ( IC .EQ. 6 .OR. ( IC .GT. 1 .AND. I0 .EQ. I0OLD .AND.
 
169
     &             J0 .EQ. J0OLD .AND. IS .EQ. ISOLD .AND.
 
170
     &             JS .EQ. JSOLD .AND. NPX .EQ. NPXOLD ) ) THEN
 
171
              LAST = .TRUE.
 
172
              IC = 6
 
173
          ENDIF
 
174
          NPXOLD = NPX
 
175
C
 
176
C ***     Check component's distance.
 
177
C
 
178
          IF ( COMP ) THEN
 
179
              XDF = PMT(10) - FLOAT(NCAT(1,LZ))
 
180
              YDF = PMT(11) - FLOAT(NCAT(2,LZ))
 
181
              IDF = NCT(1) - NCAT(1,LZ)
 
182
              JDF = NCT(2) - NCAT(2,LZ)
 
183
              IF ( ABS(IDF) .LE. LPXL .AND. ABS(JDF) .LE. LPXL ) THEN
 
184
                  CENTER = .TRUE.
 
185
C
 
186
C ***             Identify component's subpixel.
 
187
C
 
188
                  CALL SUBPXL( NCT , PMT , LSBP , ICM , JCM ,
 
189
     &                         ISC , JSC )
 
190
              ELSE
 
191
                  CENTER = .FALSE.
 
192
              ENDIF
 
193
          ELSE
 
194
              CENTER = .FALSE.
 
195
          ENDIF
 
196
C
 
197
C ***     IOFF0 - adress of psf for central pixel at given IS and JS.
 
198
C
 
199
          IOFF0 = ( (LSBP+JS) * MSBP + LSBP + IS ) * MPXL2 +
 
200
     &               LPXL * ( MPXL + 1 ) + 1
 
201
C
 
202
C ***     IOFFC - difference of adress of central
 
203
C ***     pixel at ISC and JSC minus IOFF0.
 
204
C
 
205
          IF ( CENTER ) THEN
 
206
              IOFFC = ( (JSC-JS) * MSBP + ISC - IS ) * MPXL2
 
207
          ENDIF
 
208
C
 
209
C ***     IOFF1 and IOFF2 - corresponding adresses of x and y gradients.
 
210
C
 
211
          IOFF1 = IOFF0 + IOFF
 
212
          IOFF2 = IOFF1 + IOFF
 
213
C
 
214
C ***     Check if this object can serve as standard star.
 
215
C
 
216
          LPLIM = INT( FLOAT(LPXL) * SQRT(2.0) )
 
217
          IF ( STAR .AND. LAST .AND. KCLN .GT. LPLIM
 
218
     &                         .AND. NNOR .EQ. 4 ) THEN
 
219
              IF ( SIGM1 .GT. 0.0 ) THEN
 
220
                  SCAA = SCA / SIGM1
 
221
              ELSE
 
222
                  SCAA = 0.0
 
223
              ENDIF
 
224
              CALL IFSTAR( IBUFS , SCAA , LPXL , LSBP , IS ,
 
225
     &                     JS , MPRF , IPSF , FPSF , IOFF0 ,
 
226
     &                     IOFF , STAR , IOFFI )
 
227
              IC = 7
 
228
          ENDIF
 
229
C
 
230
C ***     Calculate normal equations.
 
231
C
 
232
          CALL NRMLEQ( AS , MAS , JAPYS , IBUFS , LPXL ,
 
233
     &                 APSF , FPSF , IOFF0 , IOFF1 , IOFF2 ,
 
234
     &                 IOFFI , BGRD , I0 , J0 , SCA ,
 
235
     &                 COMP , SCC , CENTER , IDF , JDF ,
 
236
     &                 XDF , YDF , IOFFC , IC , SIGM ,
 
237
     &                 NPX , AR )
 
238
C
 
239
C ***     Solve normal equations if there is enough of good points.
 
240
C
 
241
          IF ( COMP ) THEN
 
242
              NNOR = 5
 
243
          ELSE
 
244
              NNOR = 4
 
245
          ENDIF
 
246
          IF ( NPX .GT. NNOR .AND. ( .NOT. LAST ) ) THEN
 
247
              CALL LSQSOL ( NNOR , AR , NPX , X , S )
 
248
              SIGM = SNGL(S(NNOR))
 
249
              SCA = SCA + SNGL(X(1))
 
250
              CORX = SNGL(X(2))
 
251
C
 
252
C ***         Save sigma from the first, non trimmed
 
253
C ***         determination of object intensity.
 
254
C
 
255
              IF ( IC .EQ. 1 ) THEN
 
256
                  SIGM1 = SNGL(S(1))
 
257
              ENDIF
 
258
              IF ( ABS(CORX) .GT. SHLM ) THEN
 
259
                  CORX = SHLM * CORX / ABS(CORX)
 
260
              ENDIF
 
261
              PMTR(10,LZ) = FLOAT(NCAT(1,LZ)) +
 
262
     &                      FLOAT(IS) / FLOAT(MSBP) - CORX
 
263
              PMTR(22,LZ) = SNGL(S(2))
 
264
              NCAT(1,LZ) = NINT(PMTR(10,LZ))
 
265
              CORY = SNGL(X(3))
 
266
              IF ( ABS(CORY) .GT. SHLM ) THEN
 
267
                  CORY = SHLM * CORY / ABS(CORY)
 
268
              ENDIF
 
269
              PMTR(11,LZ) = FLOAT(NCAT(2,LZ)) +
 
270
     &                      FLOAT(JS) / FLOAT(MSBP) - CORY
 
271
              PMTR(23,LZ) = SNGL(S(3))
 
272
              NCAT(2,LZ) = NINT(PMTR(11,LZ))
 
273
              IF ( COMP ) THEN
 
274
                  SCC = SCC + SNGL(X(4))
 
275
              ENDIF
 
276
          ENDIF
 
277
C
 
278
          ISOLD = IS
 
279
          JSOLD = JS
 
280
          I0OLD = I0
 
281
          J0OLD = J0
 
282
C
 
283
          IC = IC + 1
 
284
          GOTO 10
 
285
   11 CONTINUE
 
286
C
 
287
C *** Apply corrections.
 
288
C
 
289
      IF ( NPXOLD .GT. NNOR ) THEN
 
290
          PMTR(12,LZ) = SCA
 
291
          IF ( SCA .GT. 0.0 ) THEN
 
292
              PMTR(24,LZ) = 1.086 * SNGL(S(1)) / SCA
 
293
          ELSE
 
294
              PMTR(24,LZ) = 1.0
 
295
          ENDIF
 
296
          NCAT(9,LZ) = 1
 
297
          IF ( NCAT(10,LZ) .GT. 0 ) THEN
 
298
              PMT(2) = SCC / ( 9.0 / ( 1.0 + 8.0 * APSF(1) ) )
 
299
              PMT(12) = SCC
 
300
          ENDIF
 
301
      ELSE
 
302
          NCAT(9,LZ) = 0
 
303
      ENDIF
 
304
      SIGMA = SNGL(S(NNOR))
 
305
C
 
306
      RETURN
 
307
C
 
308
      END