← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/surfphot/src/dtom.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/surfphot/src/dtom.for
 
1
C @(#)dtom.for  19.1 (ES0-DMD) 02/25/03 13:31:22
 
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
        PROGRAM DTOM
 
30
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
31
C.IDENT  program DTOM                           A. Lauberts     ESO  890710
 
32
C
 
33
C.KEYWORDS  characteristic curve, Schmidt plates, standard galaxies,
 
34
C           sensitometer spots
 
35
C
 
36
C.PURPOSE
 
37
C
 
38
C  Find least squares fit of
 
39
C  MPG = photographic magnitude/2.5 (in rings between circular apertures)
 
40
C                                    + innermost aperture
 
41
C  to
 
42
C  MEL = photoelectric magnitude/2.5 (in same rings as MPG)
 
43
C                                     curve fitted
 
44
C  and simultaneously
 
45
C
 
46
C   F  = A*log(D-DF) + B*log(DS-D) + CM (the D to log I function)
 
47
C  to  
 
48
C   E  = log I + CS (versus the spot densities D)
 
49
C
 
50
C  Minimize  S = Sum (MPG-MEL)**2 + Sum (F-E)**2
 
51
C      over   min0:max9 rings       min0:max21 spots
 
52
C  where 
 
53
C
 
54
C   D = Density measured by PDS (divide by 200 or 800!) 
 
55
C  DF = Density of the chemical fog  = X(1) = XX(1)  
 
56
C  DS = Density at saturation        = X(2) = XX(2)
 
57
C   G = Gamma = dD/d(log I)
 
58
C  RG = ln10/Gamma at infl           = X(3) = XX(7)       
 
59
C  DI = Density at inflexion point   = X(4) = XX(8)
 
60
C  DB = Density at sky background
 
61
C  CM = Zero point for mag/2.5       = X(5)
 
62
C       Zero point in sky units             = XX(5)
 
63
C  CS = Zero point for spots         = X(6) = XX(9)
 
64
C   A = RG*(DI-DF)**2/(DS-DF)               = XX(3)                     
 
65
C   B =-RG*(DS-DI)**2/(DS-DF)               = XX(4)
 
66
C  SM = Sky magnitude/20 micron aperture    = XX(6)
 
67
C
 
68
C.USAGE  MIDAS applic program
 
69
C        execute as @@ DTOM
 
70
C
 
71
C.PARAMETERS
 
72
C
 
73
C Input: promted by inquire statements 
 
74
C
 
75
C  TABLEA    C*8        standard galaxy photoel table file
 
76
C                       (format same as BVRLA.TBL used for the 
 
77
C                       automated ESOLV procedures)
 
78
C  TABLEC    C*8        table file with calibration parameters
 
79
C  FRAMEA    C*8        image frame with central galaxy
 
80
C  FRAMEB    C*8        table file with sensitometer spots:
 
81
C                       SE   R*4 array  log I   intensity
 
82
C                       SD   R*4 array  D       density
 
83
C  NES       R*4 array  1st elem, no of spots wanted
 
84
C  SPOTE     R*4        constant log I difference (option)
 
85
C  NA, NR    I*4        no of apertures (rings)
 
86
C  ESO       R*4        ESO ident = (field.object no) = ESO.OBJ
 
87
C  SCAL      R*4        plate scale arcseq/mm
 
88
C  WT        R*4 array  weights of innermost, outermost, internal rings
 
89
C  SM        R*4        sky mag/20 microns aperture (option)
 
90
C  XX        R*4 array  initial calibration values (4 first elems)
 
91
C  FIX       I*4 array  1/0 for fix/no-fix calibr parameters
 
92
C
 
93
C Descriptors:
 
94
C
 
95
C  BVR       R*4 array  ident ESO.OBJ + BVR vs log A coefficients read from 
 
96
C                       TABLEA (format same as BVRLA.TBL used for the
 
97
C                       automated ESOLV procedures).
 
98
C  GALCEN    R*4 array  Galaxy centre computed in program COMSKY
 
99
C                       or set by interactive cursor
 
100
C  SKYBGR    R*4 array  Sky background computed in COMSKY
 
101
C                       = command COMUTE/SKY
 
102
C
 
103
C Output:
 
104
C
 
105
C  CAL       R*4 array  ident ESO.OBJ + calibration parameters -> TABLEC
 
106
C
 
107
C.COMMENTS
 
108
C
 
109
C Note: In case of non-existing V-R (all coeffs = 0's), a substitute
 
110
C       is taken as
 
111
C
 
112
C       V-R = 0.134 + 0.453 * (B-V)  at BS = 22 mags/sqarcs
 
113
C       d(V-R)/dlogA = -0.025 = constant,
 
114
C
 
115
C       being a good approximation (R.M.S. err = 0.02)
 
116
C       for 0.8 < B-V < 1.2 (E and S0's)
 
117
C 
 
118
C       The photoelectric values are reduced to outside the atmosphere
 
119
C       (airmass zero). Therefore, the present 'sky parameter' becomes 
 
120
C       brighter than the true sky. The magnitude difference is simply
 
121
C       the actual extinction. Whenever needed we take Bext = 0.20 (ESO
 
122
C       manual 1978) and Rext = 0.10 (ESO manual gives 0.08, own observ
 
123
C       0.12 after 1982) per unit airmass.
 
124
C
 
125
C       This program calls former NAGLIB subroutine E04GAF, 
 
126
C       which - together with supplementary smaller routines -
 
127
C       presently resides in module MNAG
 
128
C
 
129
C-------------------------------------------------------------------------------
 
130
        COMMON /VMR/MADRID(1)
 
131
        COMMON /DENS/ID(4096,9),MAXD,NR,WT(30),NS,DB,SM,A,B,XX,STEP,
 
132
     +          SD(21),SE(21),EL(9),ORG,FIX
 
133
        DOUBLE PRECISION  RR,F,X(6),X02AAF,XTOL(6),SCALE(6),R(33),W(100),
 
134
     +          A,B,DB,ORG(6)
 
135
        REAL  XX(9),APER(9),MEL(9),RADA(9),WEIGHT(9),CEL(9),
 
136
     +          STEPA(2),STARTA(2),CEN(2),LOGA(9),BVR(31),CAL(9),SEL(9)
 
137
        INTEGER  NPIXA(2),STAT,PNTRA,NES(2),COLA(31),COLC(9),FIX(6)
 
138
        CHARACTER*1   CALC,CHR
 
139
        CHARACTER*60  FRAMEA,FRAMEB,SKYMAG,TABLEA,TABLEC,LABEL(2)
 
140
        CHARACTER*20  IDENTA
 
141
        CHARACTER*24  CUNITA
 
142
        CHARACTER*9   FSSQ
 
143
        CHARACTER*72  CALIB,RESIDS,SURFB
 
144
        CHARACTER*3   ERRNO,RIN,SPO
 
145
        LOGICAL       NULL(31),IFL
 
146
C
 
147
        PARAMETER     NCOLA=31, NCOLC=9
 
148
        DATA          COLC/1,2,3,4,5,6,7,8,9/,
 
149
     +                LABEL/'DENSITY','LOG_INT'/
 
150
C
 
151
        EXTERNAL  RESID,LSQ,MONIT
 
152
C
 
153
        G(D)=C3*ALOG10(D-C1)+C4*ALOG10(C2-D)+C5
 
154
C
 
155
        DO I=1,31
 
156
          COLA(I)=I
 
157
        ENDDO
 
158
C
 
159
C  set up MIDAS environment, access tables (already existing)
 
160
C
 
161
        CALL SXSPRO('DTOM',STAT)
 
162
        CALL SXKGET('TABLEA','C',1,60,IAV,TABLEA,STAT)
 
163
        CALL SXKGET('TABLEC','C',1,60,IAV,TABLEC,STAT)
 
164
        CALL TXTOPN(TABLEA,STAT)
 
165
        CALL TXTOPN(TABLEC,STAT)
 
166
C  get plate info
 
167
        CALL SXKGET('INPUTC','C',1,1,IAV,CHR,STAT)
 
168
        CALL SXKGET('SCALE','R',1,1,IAV,SCAL,STAT)
 
169
        SKY=0.
 
170
        DB=0.
 
171
        DO 40 J=1,30
 
172
        WT(J)=1.
 
173
40      CONTINUE
 
174
        CALL SXKGET('IN_B','C',1,60,IAV,FRAMEB,STAT)
 
175
        CALL TXTOPN(FRAMEB,STAT)
 
176
        CALL TXIGET(FRAMEB,NCOLB,NROWB,IAV,STAT)
 
177
        NS=NROWB
 
178
        DO I=1,NS
 
179
          CALL TXRGET(FRAMEB,'R',I,2,COLA,XX,NULL,STAT)
 
180
          SE(I)=XX(1)
 
181
          SD(I)=XX(2)
 
182
        ENDDO
 
183
        WRITE(6,111) (SD(I),I=1,NS)
 
184
111     FORMAT(' SPOTD',10F7.3)
 
185
48      WRITE(6,*) ' first element, no of spots: '
 
186
        READ(5,*) NES
 
187
        NE=NES(1)
 
188
        NS=NES(2)
 
189
        IF(NS.EQ.1.OR.NS.GT.21) THEN
 
190
          CALL SXTPUT('NO OF SPOTS MUST BE 0 OR 2-21',STAT)
 
191
          GOTO 48
 
192
        ENDIF
 
193
        IF(NS.EQ.0) GOTO 55
 
194
        WRITE(6,*) ' const log intensity difference 
 
195
     + (if =0 then 1st col in Cxxx): '
 
196
        READ(5,*) SPOTE
 
197
49      CONTINUE
 
198
        IF(NE.NE.1) THEN
 
199
          NE=NE-1
 
200
          DO I=1,NS
 
201
            SD(I)=SD(I+NE)
 
202
            SE(I)=SE(I+NE)
 
203
          ENDDO
 
204
        ENDIF
 
205
        IF(SPOTE.GT.0.) THEN
 
206
          DO I=1,NS
 
207
            SE(I)=SPOTE*FLOAT(I)
 
208
          ENDDO
 
209
          S=-0.5
 
210
          DO I=1,NROWB
 
211
            CALL TXEPUT(FRAMEB,'R',I,1,S,STAT)
 
212
            S=S+SPOTE
 
213
          ENDDO
 
214
        ENDIF
 
215
        WRITE(6,112) (SE(I),I=1,NS)
 
216
112     FORMAT(' SPOTE',10F7.3)
 
217
        CALL TXTCLO(FRAMEB,STAT)
 
218
        NR=0
 
219
55      CALL SXKGET('NAPE','I',1,1,IAV,NA,STAT)
 
220
        IF(NA.GT.9) THEN
 
221
          CALL SXTPUT(' >9 APERTURES/OBJECT',STAT)
 
222
          GOTO 55
 
223
        ENDIF
 
224
C  get image frame 
 
225
        CALL SXKGET('IN_A','C',1,60,IAV,FRAMEA,STAT)
 
226
        CALL SXIGET(FRAMEA,'I',2,NAXISA,NPIXA,STARTA,STEPA,
 
227
     +              IDENTA,CUNITA,PNTRA,STAT)
 
228
        IF(NA.NE.0) CALL SXDGET(FRAMEA,'GALCEN','R',1,2,IAV,CEN,STAT)
 
229
        CALL SXDGET(FRAMEA,'SKYBGR','R',1,1,IAV,SKY,STAT)
 
230
        CALL SXDGET(FRAMEA,'O_TIME','R*8',1,1,IAV,A,STAT)
 
231
        IF(A.GT.1983.5D0) THEN
 
232
          MAXD=4096
 
233
          DB=SKY/800.
 
234
         ELSE
 
235
          MAXD=1024
 
236
          DB=SKY/200.
 
237
        ENDIF
 
238
C  convert sky to mags/(20*20) square micron aperture
 
239
        STEP=STEPA(1)
 
240
        IF(NA.EQ.0) THEN
 
241
          CALL SXKGET('SKYMAG','R',1,1,IAV,SM,STAT)
 
242
          SM=SM-5.*ALOG10(SCAL*0.020)
 
243
        ENDIF
 
244
        NDIM=NPIXA(1)*NPIXA(2)
 
245
        IF(NA.EQ.0) GOTO 100
 
246
        NA1=9
 
247
        S=5000./SCAL
 
248
        A1=10.**0.1
 
249
        DO I=1,NA1
 
250
          RADA(I)=S
 
251
          S=S*A1
 
252
        ENDDO
 
253
        A1=0.9
 
254
        M=100*MAXD/1024
 
255
        CALL INTEGR(NDIM,MADRID(PNTRA),STARTA,STEPA,NPIXA,CEN,RADA,
 
256
     +              NA1)
 
257
        CALL SXTPUT(' D-UP    0.5  1.0  1.5  2.0  2.5  3.0  3.5  
 
258
     +4.0  4.5  5.0',STAT)
 
259
        DO J=1,NA1
 
260
          DO I=1,M*10,M
 
261
            DO K=1,M-1
 
262
              ID(I,J)=ID(I,J)+ID(I+K,J)
 
263
            ENDDO
 
264
          ENDDO
 
265
          A1=A1+0.1
 
266
          WRITE(SURFB,204) A1, (ID(I,J), I=1,M*10,M)
 
267
          CALL SXTPUT(SURFB,STAT)
 
268
        ENDDO
 
269
        S=RADA(NA1)
 
270
        NA1=5
 
271
        A1=10.**0.1
 
272
        DO I=1,NA1
 
273
          RADA(I)=S
 
274
          S=S*A1
 
275
        ENDDO
 
276
        A1=1.8
 
277
        CALL INTEGR(NDIM,MADRID(PNTRA),STARTA,STEPA,NPIXA,CEN,RADA,
 
278
     +              NA1)
 
279
        DO J=2,NA1
 
280
          DO I=1,M*10,M
 
281
            DO K=1,M-1
 
282
              ID(I,J)=ID(I,J)+ID(I+K,J)
 
283
            ENDDO
 
284
          ENDDO
 
285
          A1=A1+0.1
 
286
          WRITE(SURFB,204) A1, (ID(I,J), I=1,M*10,M)
 
287
          CALL SXTPUT(SURFB,STAT)
 
288
        ENDDO
 
289
        WRITE(6,*) ' first, last log aperture: '
 
290
        READ(5,*) LOGA(1),LOGA(2)
 
291
C  specify table entry
 
292
60      CALL SXKGET('ESOOBJ','R',1,1,IAV,ESO,STAT)
 
293
        APER(1)=10.**LOGA(1)
 
294
        RADA(1)=APER(1)*500./SCAL
 
295
       IF(NA.GT.1) THEN
 
296
         S=(LOGA(2)-LOGA(1))/FLOAT(NA-1)
 
297
         A1=10.**S
 
298
        DO I=2,NA
 
299
         APER(I)=APER(I-1)*A1
 
300
         LOGA(I)=LOGA(I-1)+S
 
301
         RADA(I)=RADA(I-1)*A1
 
302
        ENDDO
 
303
       ENDIF
 
304
C  store density distribution in common block /DENS/ID
 
305
C  new RADA out!
 
306
        CALL INTEGR(NDIM,MADRID(PNTRA),STARTA,STEPA,NPIXA,CEN,RADA,
 
307
     +              NA)
 
308
C  refine apertures considering actual no of pixels
 
309
        DO I=1,NA
 
310
          APER(I)=RADA(I)*SCAL*0.002
 
311
          LOGA(I)=ALOG10(APER(I))
 
312
        ENDDO
 
313
C  read BVR vs log A polynomial fit coefficients 
 
314
C  BVR(1) = ESO ident in the fld.obj no form 
 
315
C  Example: 271.012 means ESO field = 271, 12 = object no in ESO/Upp cat
 
316
C
 
317
        ERR=0.0005
 
318
        CALL TXESER(TABLEA,'R*4',1,ESO,ERR,1,NSEQ,STAT)
 
319
        CALL TXRGET(TABLEA,'R*4',NSEQ,NCOLA,COLA,BVR,NULL,STAT)
 
320
C  compute array of U,B,V or R values vs aperture A
 
321
C  as well as surface brightness in mags/sqarcsec
 
322
        CALL POLY('B',BVR,LOGA,MEL,SEL,NA)
 
323
        CALL SXTPUT('B SURFBR MAGS/SQARCSEC',STAT)
 
324
        WRITE(SURFB,201) (SEL(I), I=1,NA)
 
325
        CALL SXTPUT(SURFB,STAT)
 
326
        IF(CHR.NE.'B') THEN
 
327
          CALL POLY(CHR,BVR,LOGA,MEL,CEL,NA)
 
328
          CALL SXTPUT(' B-'//CHR,STAT)
 
329
          DO I=1,NA
 
330
            CEL(I)=SEL(I)-CEL(I)
 
331
          ENDDO
 
332
          WRITE(SURFB,201) (CEL(I), I=1,NA)
 
333
          CALL SXTPUT(SURFB,STAT)
 
334
        ENDIF
 
335
        CALL SXTPUT(' AT APERTURES ARCSEC',STAT)
 
336
        WRITE(SURFB,201) (APER(I), I=1,NA)
 
337
        CALL SXTPUT(SURFB,STAT)
 
338
C  take .4 times magnitude of inner aperture
 
339
        EL(1)=.4*MEL(1)
 
340
C  transform magnitudes to intensities
 
341
        DO 80 I=1,NA
 
342
80      MEL(I)=10.**(-.4*MEL(I))
 
343
        IF(NA.EQ.1) GOTO 100
 
344
C  take negative log of intensities in rings bound by apertures
 
345
        DO 90 I=2,NA
 
346
90      EL(I) = - ALOG10(MEL(I)-MEL(I-1))
 
347
100     NR=NA
 
348
C  free this area for next input
 
349
        CALL SXFCLO(FRAMEA,STAT)
 
350
120     N=6
 
351
        M=NR+NS
 
352
        IF(NR.EQ.0) GOTO 125
 
353
        CALL SXKGET('WEIGHT','R',1,9,IAV,WEIGHT,STAT)
 
354
        WT(1)=WEIGHT(1)
 
355
        WT(NA)=WEIGHT(2)
 
356
        DO I=2,NA-1
 
357
          WT(I)=WEIGHT(I+1)
 
358
        ENDDO
 
359
C  equal importance photoel vs photographic
 
360
        S=1.
 
361
        IF(NS.NE.0) S=NS/FLOAT(NR)
 
362
        DO I=1,NR
 
363
          WT(I)=WT(I)*S
 
364
        ENDDO
 
365
C  get initial values XX
 
366
125     CALL SXKGET('INPUTR','R',1,4,IAV,XX,STAT)
 
367
        CALL SXKGET('INPUTI','I',1,N,IAV,FIX,STAT)
 
368
      DO I=1,4
 
369
        X(I)=XX(I)
 
370
      ENDDO
 
371
        X(5)=0.
 
372
        X(6)=0.
 
373
        CALL RESID(M,N,X,R,IFL)
 
374
        IF(NR.EQ.0) GOTO 130
 
375
      A1=0.
 
376
      DO I=1,NR
 
377
        A1=A1+WT(I)
 
378
        X(5)=X(5)-R(I)
 
379
      ENDDO
 
380
        X(5)=X(5)/A1
 
381
130     CONTINUE
 
382
      A1=0.
 
383
      DO I=NR+1,M
 
384
        A1=A1+WT(I)
 
385
        X(6)=X(6)-R(I)
 
386
      ENDDO
 
387
        IF(NS.NE.0) X(6)=X(6)/A1
 
388
        XTOL(1)=DSQRT(X02AAF(RR))
 
389
      DO I=1,N
 
390
        ORG(I)=X(I)
 
391
      ENDDO
 
392
      DO I=2,N
 
393
        XTOL(I)=XTOL(1)
 
394
      ENDDO
 
395
        METHOD=1
 
396
        IW=N*(N+4)+M    
 
397
        IPRINT=1    
 
398
        MAXCAL=20      
 
399
        IFAIL=0
 
400
        CALL SXTPUT('E04GAF RESULTS:',STAT)
 
401
        CALL E04GAF(M,N,X,R,F,XTOL,METHOD,SCALE,W,IW,RESID,LSQ,   
 
402
     +              MONIT,IPRINT,MAXCAL,IFAIL)
 
403
160     WRITE(ERRNO,202) IFAIL
 
404
        CALL SXTPUT('THIS HAS ERROR NUMBER'//ERRNO,STAT)
 
405
C  update calibration parameters in TABLEC (one row only)
 
406
        CALL SXKGET('ESOFLD','I',1,1,IAV,NSEQ,STAT)
 
407
        IF(NA.NE.0.) THEN
 
408
          CAL(1)=ESO
 
409
          ELSE
 
410
          CAL(1)=NSEQ
 
411
        ENDIF
 
412
        I=ESO
 
413
        IF(I.NE.NSEQ) TYPE *,ESO,' not in field ',NSEQ
 
414
       DO I=1,8
 
415
        CAL(I+1)=XX(I)
 
416
       ENDDO
 
417
        CALL TXRPUT(TABLEC,'R',1,NCOLC,COLC,CAL,STAT)
 
418
999     CALL TXTCLO(TABLEA,STAT)
 
419
        CALL TXTCLO(TABLEC,STAT)
 
420
        CALL SXSEPI
 
421
200     FORMAT(E9.3)
 
422
201     FORMAT(9F8.3)
 
423
202     FORMAT(I3)
 
424
203     FORMAT(F8.2)
 
425
204     FORMAT(1X,F3.1,3X,10I5)
 
426
        END                                                     
 
427
C   
 
428
C    
 
429
      SUBROUTINE RESID(M,N,XC,RC,IFL)     
 
430
C  CALLED BY EO4GAF AND LSQ 
 
431
C  CALCULATES THE VALUES OF THE RESIDUALS RC AT XC
 
432
        COMMON /DENS/ID(4096,9),MAXD,NR,WT(30),NS,DB,SM,A,B,XX,STEP,
 
433
     +          SD(21),SE(21),EL(9),ORG,FIX
 
434
      LOGICAL IFL        
 
435
      DOUBLE PRECISION  XC,RC,Q,X41,X24,XC1,XC2,XC4,X,A,B,IP,DB,ORG(6)
 
436
      DIMENSION XC(N),RC(M),XX(9)
 
437
      INTEGER FIX(6)
 
438
C
 
439
        F(X) = A*DLOG10(X-XC1) + B*DLOG10(XC2-X)
 
440
        DO I=1,N
 
441
          IF(FIX(I)) XC(I)=ORG(I)
 
442
        ENDDO
 
443
        XC1=XC(1)
 
444
        XC2=XC(2)
 
445
        XC4=XC(4)
 
446
        Q=XC(3)/(XC2-XC1)
 
447
        X41=XC4-XC1
 
448
        X24=XC2-XC4
 
449
      IF(XC1.GT.DB*.1) GOTO 3
 
450
      IF(X41.GT.0.D0.AND.X24.GT.0.D0) GOTO 4
 
451
        GOTO 3
 
452
4       A= Q*X41**2
 
453
        B=-Q*X24**2
 
454
        IF(NR.EQ.0) THEN
 
455
          XC(5)=0.D0
 
456
          GOTO 11
 
457
        ENDIF
 
458
      DO 1 I=1,NR
 
459
        CALL INTRING(XC,I,IP,IFL)
 
460
        IF(IFL.EQ..TRUE.) RETURN
 
461
1     RC(I) = ( DLOG10(IP) + EL(I) + XC(5) ) * WT(I)
 
462
11      CONTINUE
 
463
      DO 2 I=1,NS
 
464
        X=SD(I)
 
465
        IF(X.GT.XC1.AND.X.LT.XC2) GOTO 2
 
466
        GOTO 3
 
467
2     RC(I+NR) = ( F(X) - SE(I) +XC(6) ) * WT(I+NR)
 
468
        IFL=.FALSE.
 
469
        RETURN
 
470
3       IFL=.TRUE.
 
471
        CALL SXTPUT('IFL=.TRUE. FROM RESID',STAT)
 
472
        RETURN              
 
473
        END                                                                  
 
474
C  
 
475
C   
 
476
      SUBROUTINE LSQ(M,N,XC,RC,AJTJC,GC)    
 
477
C  CALCULATES THE JACOBIAN RC AT XC   
 
478
C  CALLED FROM E04GAF    
 
479
      DOUBLE PRECISION  XC,RC,AJTJC,GC,AJC,SUM,DELTA,DELTA2,XJ,YFIT
 
480
      DIMENSION XC(N),RC(M),GC(N),AJTJC(N,N),AJC(30,6),YFIT(30)     
 
481
      LOGICAL IFL
 
482
      PARAMETER (DELTA=.001D0)
 
483
      PARAMETER (DELTA2=.002D0)
 
484
C  NON ANALYTICAL DERIVATIVE
 
485
      DO 3 J=1,N
 
486
      XJ=XC(J)
 
487
      XC(J)=XJ+DELTA
 
488
      CALL RESID(M,N,XC,RC,IFL)
 
489
      DO 1 I=1,M
 
490
1     YFIT(I)=RC(I)
 
491
      XC(J)=XJ-DELTA
 
492
      CALL RESID(M,N,XC,RC,IFL)
 
493
      DO 2 I=1,M
 
494
2     AJC(I,J)=(YFIT(I)-RC(I))/DELTA2
 
495
3     XC(J)=XJ
 
496
      DO 80 I=1,N    
 
497
         SUM=0.    
 
498
         DO 20 K=1,M   
 
499
            SUM=SUM+AJC(K,I)*RC(K)     
 
500
20       CONTINUE     
 
501
         GC(I)=SUM      
 
502
         II=I    
 
503
         DO 60 J=1,II   
 
504
            SUM=0.  
 
505
            DO 40 K=1,M  
 
506
               SUM=SUM+AJC(K,I)*AJC(K,J)  
 
507
40          CONTINUE  
 
508
            AJTJC(J,I)=SUM  
 
509
60       CONTINUE  
 
510
80    CONTINUE     
 
511
      RETURN   
 
512
      END    
 
513
C     
 
514
C   
 
515
      SUBROUTINE MONIT(M,N,XC,RC,FC,GC,NCALL)     
 
516
C  CALLED BY E04GAF   
 
517
C  PRINTS THE VALUES EVERY PRINT ITERATIONS
 
518
        COMMON /DENS/ID(4096,9),MAXD,NR,WT(30),NS,DB,SM,A,B,XX,STEP
 
519
        DOUBLE PRECISION  XC,GC,RC,FC,DB,A,B                         
 
520
        DIMENSION XC(N),GC(N),RC(M),R(30),XX(9)
 
521
        CHARACTER*3   CALL,RIN,SPO
 
522
        CHARACTER*9   SSQ
 
523
        CHARACTER*72  CALIB,RESIDS
 
524
        WRITE(CALL,200) NCALL
 
525
200     FORMAT(I3)
 
526
        WRITE(SSQ,201) FC
 
527
201     FORMAT(E9.3)
 
528
        CALL SXTPUT('AFTER'//CALL//' CALLS OF RESID, THE SUM OF
 
529
     + SQUARES IS '//SSQ//' AT',STAT)
 
530
        CALL SXTPUT('   DF      DS      A       B       CM
 
531
     +      SM      RG      DI',STAT)
 
532
      DO 1 I=1,5
 
533
1     XX(I)=XC(I)
 
534
        XX(3)=A
 
535
        XX(4)=B
 
536
        XX(9)=XC(6)
 
537
C  scale intensities to sky units through adjusting XX(5)
 
538
        S=DB
 
539
        IF(S.LE.XX(1)) THEN
 
540
          WRITE(6,*) ' SKYBGR < FOG!'
 
541
          RETURN
 
542
        ENDIF
 
543
        S=A*ALOG10(S-XX(1))+B*ALOG10(XX(2)-S)+XX(5)
 
544
        XX(5)=XX(5)-S
 
545
C  sky magnitude/20 micron aperture
 
546
        IF(NR.NE.0) THEN
 
547
          XX(6)=-2.5*(S - 2.*ALOG10(STEP/20.))
 
548
         ELSE
 
549
          XX(6)=SM
 
550
        ENDIF
 
551
        XX(7)=XC(3)
 
552
        XX(8)=XC(4)
 
553
        WRITE(CALIB,202) (XX(I),I=1,8)
 
554
202     FORMAT(9F8.3)
 
555
        CALL SXTPUT(CALIB,STAT)
 
556
        WRITE(RIN,200) NR
 
557
        WRITE(SPO,200) NS
 
558
        CALL SXTPUT('WITH RESIDUES OF'//RIN//' RINGS AND'
 
559
     +//SPO//' SPOTS',STAT)
 
560
      DO I=1,M
 
561
        R(I)=RC(I)/WT(I)
 
562
      ENDDO
 
563
        IF(NR.NE.0) THEN
 
564
          WRITE(RESIDS,202) (R(I), I=1,NR)
 
565
          CALL SXTPUT(RESIDS,STAT)
 
566
        ENDIF
 
567
        DO 20 I1=NR+1,37,9
 
568
        IF(M.LT.I1) GOTO 30
 
569
        I2=I1+8
 
570
        IF(M.LT.I2) I2=M
 
571
        WRITE(RESIDS,202) (R(I), I=I1,I2)
 
572
20      CALL SXTPUT(RESIDS,STAT)
 
573
30      CONTINUE
 
574
      RETURN
 
575
      END
 
576
C
 
577
C
 
578
        SUBROUTINE INTRING(X,IR,IP,IFL)
 
579
C  called from RESID
 
580
C  computes integrated intensities in rings
 
581
        COMMON /DENS/ID(4096,9),MAXD,NR,WT(30),NS,DB,SM,A,B
 
582
        DOUBLE PRECISION  X(6),A,B,DB,DF,DF1,DS,DS1,DK,DI,IP,IB,IT
 
583
        LOGICAL IFL
 
584
C
 
585
        IF(MAXD.EQ.4096) THEN
 
586
          DI=0.00125D0
 
587
        ELSE
 
588
          DI=0.005D0
 
589
        ENDIF
 
590
        DF=X(1)
 
591
        DS=X(2)
 
592
        DF1=DF+.01D0
 
593
        DS1=DS-.01D0
 
594
        IF(DB.GT.DF1.AND.DB.LT.DS1) GOTO 4
 
595
        GOTO 2
 
596
4       IB = (DB-DF)**A * (DS-DB)**B
 
597
        IP=0.
 
598
        DK=0.
 
599
      DO 1 I=1,MAXD
 
600
        II=ID(I,IR)
 
601
        IF(II.EQ.0) GOTO 1
 
602
        IF(DK.GT.DF1.AND.DK.LT.DS1) GOTO 3
 
603
        GOTO 2
 
604
3       IT = (DK-DF)**A * (DS-DK)**B
 
605
        IP = IP + (IT-IB)*FLOATJ(II)
 
606
1     DK=DK+DI
 
607
        IFL=.FALSE.
 
608
        IF(IP.GT.0.D0) RETURN
 
609
2       IFL=.TRUE.
 
610
        CALL SXTPUT('IFL=.TRUE. FROM INTRING',STAT)
 
611
        RETURN
 
612
        END
 
613
C
 
614
C
 
615
        SUBROUTINE INTEGR(NDIM,AA,STARTA,STEPA,NPIXA,CEN,RADA,NA)
 
616
C  called from DTOM
 
617
C  calculates distribution of densities in ring-shaped regions
 
618
        COMMON /DENS/ID(4096,9),MAXD
 
619
        REAL  AA(NDIM),STARTA(2),STEPA(2),CEN(2),RADA(9),RAD2(9)
 
620
        INTEGER   NPIXA(2),NUMR(9)
 
621
        DO J=1,NA
 
622
          DO I=1,MAXD
 
623
            ID(I,J)=0
 
624
          ENDDO
 
625
          NUMR(J)=0
 
626
        ENDDO
 
627
      DO I=1,NA
 
628
        RAD2(I)=RADA(I)**2
 
629
      ENDDO
 
630
C  put restrictions on max region to search
 
631
        X0=CEN(1)
 
632
        Y0=CEN(2)
 
633
        RAD=RADA(NA)
 
634
        X1=X0-RAD
 
635
        X2=X0+RAD
 
636
        Y1=Y0-RAD
 
637
        Y2=Y0+RAD
 
638
        S=1.
 
639
        J1=(Y1-STARTA(2))/STEPA(2) -S
 
640
        J2=(Y2-STARTA(2))/STEPA(2) +S +2.
 
641
        IF(J1.LT.1) J1=1
 
642
        IF(J2.GT.NPIXA(2)) J2=NPIXA(2)
 
643
        IF(J1.GE.J2) RETURN
 
644
        Y1=STARTA(2)+(J1-1)*STEPA(2)
 
645
        I1=(X1-STARTA(1))/STEPA(1) -S
 
646
        I2=(X2-STARTA(1))/STEPA(1) +S +2.
 
647
        IF(I1.LT.1) I1 = 1
 
648
        IF(I2.GT.NPIXA(1)) I2 = NPIXA(1)
 
649
        IF(I1.GE.I2) RETURN
 
650
        X1=STARTA(1)+(I1-1)*STEPA(1)
 
651
        Y=Y1
 
652
        DO 3 J = J1,J2
 
653
        LOC=(J-1)*NPIXA(1)+I1
 
654
        YY=(Y-Y0)**2
 
655
        X=X1
 
656
        DO 2 I = I1,I2
 
657
        D2=(X-X0)**2 + YY
 
658
        DO 4 IR=1,NA
 
659
        IF(D2.GT.RAD2(IR)) GOTO 4
 
660
        IL=JIFIX(AA(LOC)+1.)
 
661
        IL=JMIN0(IL,MAXD)
 
662
C  form distribution of densities in ringshaped regions
 
663
        ID(IL,IR)=ID(IL,IR)+1
 
664
        NUMR(IR)=NUMR(IR)+1
 
665
        GOTO 1
 
666
4       CONTINUE
 
667
1       LOC=LOC+1
 
668
2       X=X+STEPA(1)
 
669
3       Y=Y+STEPA(2)
 
670
C  refine radii considering actual no of pixels
 
671
        IF(NA.GT.1) THEN
 
672
          DO I=NA,2,-1
 
673
            DO J=1,I-1
 
674
                NUMR(I)=NUMR(I)+NUMR(J)
 
675
            ENDDO
 
676
          ENDDO
 
677
        ENDIF
 
678
        DO I=1,NA
 
679
          RADA(I)=SQRT(NUMR(I)/3.141592)*STEPA(1)
 
680
        ENDDO
 
681
        RETURN
 
682
        END
 
683
C
 
684
        SUBROUTINE POLY ( T,C,LOGA,MEL,SEL,NA )
 
685
C  compute magnitudes MEL and mean surface brightess SEL
 
686
C
 
687
        CHARACTER*1  T
 
688
        REAL         C(1),MEL(NA),SEL(NA),LOGA(9)
 
689
        DO I=1,NA
 
690
          X=LOGA(I)
 
691
          IF(T.EQ.'B') THEN
 
692
            IF(C(18).EQ.0.) THEN
 
693
                C(14)=-0.3456+C(20)*2.2502
 
694
                C(15)=-0.0344+C(21)*2.1358
 
695
                C(16)= 0.0451+C(22)*2.1070
 
696
                C(17)=-0.0143+C(23)*2.0661
 
697
            ENDIF
 
698
            MEL(I)=C(2)+C(14) + X*(C(3)+C(15) + X*(C(4)+C(16) +
 
699
     +             X*(C(5)+C(17))))
 
700
           ELSE IF(T.EQ.'R') THEN
 
701
            IF(C(24).EQ.0.) THEN
 
702
                C(20)= 0.1536+C(14)*0.4444
 
703
                C(21)= 0.0161+C(15)*0.4682
 
704
                C(22)=-0.0214+C(16)*0.4746
 
705
                C(23)= 0.0069+C(17)*0.4840
 
706
            ENDIF
 
707
            MEL(I)=C(2)-C(20) + X*(C(3)-C(21) + X*(C(4)-C(22) +
 
708
     +             X*(C(5)-C(23))))
 
709
           ELSE IF(T.EQ.'U') THEN
 
710
            MEL(I)=C(2)+C(8)+C(14)+X*(C(3)+C(9)+C(15)+
 
711
     +             X*(C(4)+C(10)+C(16)+X*(C(5)+C(11)+C(17))))
 
712
           ELSE IF(T.EQ.'V') THEN
 
713
            MEL(I)=C(2)+X*(C(3)+X*(C(4)+X*C(5)))
 
714
           ELSE
 
715
            CALL SXTPUT(' LEADING CHAR MUST BE U, B, V or R !',STAT)
 
716
            RETURN
 
717
          ENDIF
 
718
          IF(I.NE.1.AND.MEL(I).LT.MEL(I1)) THEN
 
719
            SEL(I)= MEL(I) +5.*X -0.262
 
720
     +             -2.5*ALOG10(1.-10.**(-.4*(MEL(I1)-MEL(I))))
 
721
     +             +2.5*ALOG10(1.-10.**(2.*(X1-X)))
 
722
           ELSE
 
723
            SEL(I)=MEL(I) +5.*X -0.262
 
724
          ENDIF
 
725
          I1=I
 
726
          X1=X
 
727
        ENDDO
 
728
        END