← Back to branch summary

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

« back to all changes in this revision

Viewing changes to contrib/romafot/src/rfotrfit.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/romafot/src/rfotrfit.for
 
1
C @(#)rfotrfit.for      19.1 (ES0-DMD) 02/25/03 13:30:15
 
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 FIT
 
30
C---
 
31
C.IDENTIFICATION:  RFOTFIT
 
32
C.PURPOSE:         Determine characteristics of object by non-linear fitting
 
33
C.                 Data and results are taken from and written the int. table
 
34
C.AUTHOR:          R. Buonanno, G. Buscema, C. Corsi, I. Ferraro, G. Iannicola
 
35
C.                 Osservatorio Astronomico di Roma
 
36
C.VERSION:         05/14/87 date of creation
 
37
C.VERSION:         15/12/87 RHW implementation in MIDAS standard interfaces
 
38
C.VERSION:         10/03/88 RB  new version 
 
39
C.VERSION:         03/05/89 RHW change to ST interfaces and streamlining code
 
40
C.VERSION:         22/01/91 RHW all variables declared; implicit none added
 
41
C+++
 
42
       IMPLICIT    NONE
 
43
       INCLUDE     'MID_REL_INCL:RFOTDECL.INC'
 
44
C
 
45
       INTEGER     NDXY
 
46
       INTEGER     NCA
 
47
       INTEGER     MBB
 
48
       INTEGER     MPA
 
49
       INTEGER     NCD
 
50
       INTEGER     NTM
 
51
       INTEGER     NOB
 
52
       INTEGER     M12
 
53
       INTEGER     IDP
 
54
       INTEGER     IDB
 
55
       PARAMETER   (NDXY=55)
 
56
       PARAMETER   (NCA=200)
 
57
       PARAMETER   (MBB=256)
 
58
       PARAMETER   (MPA=10000)
 
59
       PARAMETER   (NCD=MBB-18)
 
60
       PARAMETER   (NTM=NCD/6.5)
 
61
       PARAMETER   (NOB=NCD/3)
 
62
       PARAMETER   (M12=NDXY*NDXY)
 
63
       PARAMETER   (IDP=(NTM*4)+3)
 
64
       PARAMETER   (IDB=NOB*3)
 
65
C
 
66
       INTEGER      ODISPL,NDISPL
 
67
       INTEGER      EC, ED, EL
 
68
       INTEGER      FL, GR
 
69
       INTEGER      I, IS, IH, IQ
 
70
       INTEGER      IPX, IPY
 
71
       INTEGER      IOBJ, IROW, IDUM
 
72
       INTEGER      IKKL, IFLG
 
73
       INTEGER*8    IPNTR
 
74
       INTEGER      IMF
 
75
       INTEGER      IP, IPA
 
76
       INTEGER      IND, INK, INK1
 
77
       INTEGER      ISTZE(100),ISTUN(100),ISTMA(100)
 
78
       INTEGER      IVX(M12),IVY(M12)
 
79
       INTEGER      ISC(7), ISX, ISY
 
80
       INTEGER      ISTAT, IAC
 
81
       INTEGER      IDI, IND1, IND2
 
82
       INTEGER      IRINT
 
83
       INTEGER      IDE, ICB
 
84
       INTEGER      JL0, JL1, JL2, JIL4
 
85
       INTEGER      J5,  JIL, JOLD, JQ4
 
86
       INTEGER      JQ, JS
 
87
       INTEGER      KFR(NTM)
 
88
       INTEGER      KUN, KNUL
 
89
       INTEGER      K3, KP, KG, KKL, KL, KTP, KTS, KRO, K13
 
90
       INTEGER      KAL, KSF, KPR1, KAS
 
91
       INTEGER      K6, K7, K9, K74
 
92
       INTEGER      MADRID(1)
 
93
       INTEGER      MAXX
 
94
       INTEGER      LI, LLI
 
95
       INTEGER      L2, L3, L4, L5, L6, L7
 
96
       INTEGER      LJ1, LJ3
 
97
       INTEGER      LG, LFL2
 
98
       INTEGER      LFLAG
 
99
       INTEGER      MASK(NDXY,NDXY)
 
100
       INTEGER      NCINT,NRINT,NSINT,NACINT,NARINT
 
101
       INTEGER      NOBJ, NGRP, NINT, NAXIS
 
102
       INTEGER      NPIX(3)
 
103
       INTEGER      NONF(NTM)
 
104
       INTEGER      NPL, NL, NPS
 
105
       INTEGER      NNS, NNF
 
106
       INTEGER      NITER, NC, NP, NCOM, NCO
 
107
       INTEGER      NCM, NUX, NRFR
 
108
       INTEGER      NHL, NCP, NCOMM, NCAL
 
109
       INTEGER      NPSG
 
110
       INTEGER      RM(NCA)
 
111
       INTEGER      TNIT, TINULL
 
112
       INTEGER      TIDINT
 
113
       INTEGER      WFLAG
 
114
C      
 
115
       INTEGER      PREVET,PREMAS,FRASUZ,RIPVET,CALINT,AZZMAS,CONPAR
 
116
C
 
117
       DOUBLE PRECISION START(3),STEP(3)
 
118
       DOUBLE PRECISION TDNULL,TDTRUE,TDFALS,DDUM
 
119
C
 
120
       REAL         A, ALAL, ABC, ALT1
 
121
       REAL         AIN, ANC
 
122
       REAL         B, BET,BETA, BEBE
 
123
       REAL         D, D1, D2, D3, D4, D6
 
124
       REAL         DELT, DPS, DXY2
 
125
       REAL         EEE
 
126
       REAL         FOG, FON, FON1, FOG1, FAT
 
127
       REAL         GIG
 
128
       REAL         HH, HGU
 
129
       REAL         P(IDP),PES(7),SIGP(1000),RIA(NDXY)
 
130
       REAL         PSA, PSAS, PAS
 
131
       REAL         PA(IDP),PB(IDB)
 
132
       REAL         RAG(NTM)
 
133
       REAL         RVAL(5)
 
134
       REAL         RK4, RMIY, RKN, RR, RINTQ, RMAY
 
135
       REAL         SAT, SIGMA, SOFOT, SIGG, SCAP, SAT1
 
136
       REAL         SLU, SQM , SIGQ, SQMS, SMT, SIGMA1
 
137
       REAL         TRNULL, TBLSEL
 
138
       REAL         U, US
 
139
       REAL         USC(NTM),USE(NTM)
 
140
       REAL         V, VAL, VG
 
141
       REAL         VALME(M12)
 
142
       REAL         VEZE(MPA),VEUN(MPA),VEMA(MPA)
 
143
       REAL         VZ(M12)
 
144
       REAL         WEI(M12)
 
145
C       
 
146
       CHARACTER    FITMET*20,CBETA*80,FITOPT*4,MOF*1,WE*1,AVEOPT*1
 
147
       CHARACTER    DIS*1
 
148
       CHARACTER    FRAME*60,IDENT*72,CUNIT*72
 
149
       CHARACTER    INTFIL*60
 
150
       CHARACTER*80 STRING
 
151
       CHARACTER*60 XXX
 
152
       CHARACTER    INME*1
 
153
       CHARACTER    SIGFI*1,SKYFI*1,POSFI*1,IUT*1,TILT*1
 
154
C
 
155
       LOGICAL      RTELOG
 
156
       LOGICAL      SFLAG
 
157
C
 
158
       INCLUDE      'MID_INCLUDE:ST_DEF.INC'
 
159
       COMMON       /VMR/MADRID
 
160
       INCLUDE      'MID_INCLUDE:ST_DAT.INC'
 
161
 
 
162
       DATA         PES/1.,1.,2.,2.,1.,1.,2./
 
163
       DATA         D/0.7/
 
164
       DATA         ISC/7*1/
 
165
       DATA         ISTZE/100*0/
 
166
       DATA         ISTUN/100*0/
 
167
       DATA         ISTMA/100*0/
 
168
 
 
169
 9001  FORMAT('Frame:              ',A40)
 
170
 9002  FORMAT('Intermediate table: ',A40)
 
171
 9003  FORMAT('Sigma:     ', G13.6,'; Saturation:   ', G13.6) 
 
172
 9004  FORMAT('Tolerance: ', G13.6,'; # Iterations: ', I5)
 
173
 9005  FORMAT('Fit method: MOFFET, UNIFORM weighting; BETA: ',G13.6)
 
174
 9006  FORMAT('Fit method: MOFFET, 1/N weighting; BETA: ', G13.6)
 
175
 9007  FORMAT('Fix sigma, fix position, fix sky, tilt sky plane: ',A)
 
176
C
 
177
 9011  FORMAT('Too few pixel in this window')
 
178
 9012  FORMAT('Window identification: ',I6,2X,
 
179
     2        ' Components: ',I2,2X,' Iterations: ',I2)
 
180
 9013  FORMAT('Fit errors:      ',4F9.4)
 
181
 9014  FORMAT(3F9.4)
 
182
C
 
183
 9021  FORMAT(80('-'))
 
184
 9022  FORMAT('Components not examined: ',I6)
 
185
 9023  FORMAT('Total  Mag. index      Histogram')
 
186
 9024  FORMAT(1X,I5,2X,F5.1,1X,F11.0,1X,'I',A50)
 
187
 9025  FORMAT(1X,I5,2X,F5.1,1X,F11.0,1X,'I')
 
188
 9026  FORMAT('Components succesfully fitted:  ',I6)
 
189
 9027  FORMAT('Components with no convergency: ',I6)
 
190
 
 
191
 9031  FORMAT('Iteration: ',I6,'; Added star: ',I6)
 
192
 9032  FORMAT('Component below the photometric threshold')
 
193
 9035  FORMAT('N O  C O N V E R G E N C E  on ',I3,' components')
 
194
 9036  FORMAT(1X,I6,' Component Parameters:        ',F15.1,2F9.1,F7.2)
 
195
 9037  FORMAT('N O  C O N V E R G E N C E')
 
196
 
 
197
C
 
198
C *** Start the program *****************************************************
 
199
       CALL STSPRO('FIT')
 
200
       CALL TBMNUL(TINULL,TRNULL,TDNULL)
 
201
       CALL TBMCON(TBLSEL,TDTRUE,TDFALS)
 
202
 
 
203
       JL0  = 0
 
204
       JL1  = 0
 
205
       JL2  = 0
 
206
       RMAY = -10.**15
 
207
       RMIY = -RMAY
 
208
C
 
209
C *** read the frame
 
210
       CALL STKRDC('IN_A',1,1,60,IAC,FRAME,KUN,KNUL,ISTAT)
 
211
       CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,3,NAXIS,
 
212
     2             NPIX,START,STEP,IDENT,CUNIT,IPNTR,IMF,ISTAT)
 
213
       NPL   = NPIX(1)
 
214
       NL    = NPIX(2)
 
215
C
 
216
C *** read the intermediate file
 
217
       CALL STKRDC('IN_B',1,1,60,IAC,STRING,KUN,KNUL,ISTAT)  ! intermediate file
 
218
       IFLG = INDEX(INTFIL,'/')                              ! look for the flag
 
219
       IF (IFLG.EQ.0) THEN
 
220
          INTFIL = STRING
 
221
       ELSE
 
222
          IUT = 'F'
 
223
          INTFIL = STRING(1:IFLG-1)
 
224
       ENDIF
 
225
       CALL STECNT('GET',EC,ED,EL)
 
226
       CALL STECNT('PUT',1,0,0)
 
227
       CALL TBTOPN(INTFIL,F_IO_MODE,TIDINT,ISTAT)
 
228
       IF (ISTAT.NE.0) THEN
 
229
          STRING = '*** FATAL: Problems with opening intermediate'//
 
230
     2             ' table ... '
 
231
          CALL STTPUT(STRING,ISTAT)
 
232
          CALL STSEPI
 
233
       ENDIF
 
234
C
 
235
C *** get table info
 
236
       CALL TBIGET(TIDINT,NCINT,NRINT,NSINT,NACINT,NARINT,ISTAT)
 
237
       IF (ISTAT.NE.0) THEN
 
238
          STRING = '*** FATAL: Problems with getting info for '//
 
239
     2             ' intermediate table; Try again ... '
 
240
          CALL STTPUT(STRING,ISTAT)
 
241
          CALL STSEPI
 
242
       ENDIF
 
243
       IF (NRINT.EQ.0) THEN
 
244
          STRING = '*** FATAL: No data points in intermediate table'
 
245
          CALL STTPUT(STRING,ISTAT)
 
246
          CALL STSEPI
 
247
       ENDIF
 
248
       CALL STECNT('PUT',EC,ED,EL)
 
249
C
 
250
C *** read table descriptor
 
251
       CALL INTDRD(TIDINT,NGRP,NOBJ,NINT,
 
252
     2             SAT,FAT,SIGMA,BETA,SOFOT,AIN,FOG)
 
253
C
 
254
C *** get threshold, sky background
 
255
       CALL STKRDR('INPUTR',1,2,IAC,RVAL,KUN,KNUL,ISTAT)    !thresh., sky backgr
 
256
       IF (RVAL(1) .GT. 0.0) THEN
 
257
          SOFOT = RVAL(1)
 
258
       ENDIF
 
259
C
 
260
       FON1    = RVAL(2)
 
261
       FOG1    = 0.0
 
262
       IF (FON1.NE.0.0) THEN
 
263
          FON  = FON1
 
264
          FOG  = FOG1
 
265
       ENDIF
 
266
C
 
267
C      FON1    = RVAL(2)
 
268
C      ALFA     = 200.0
 
269
C      FOG1     = 0.0
 
270
C      IF (FOG1.NE.0.0 .OR. FON1.NE.0.0 .OR. A1.NE.0.0) THEN
 
271
C         A1   = 200.
 
272
C         FON  = FON1
 
273
C         ALFA = A1
 
274
C         FOG  = FOG1
 
275
C      ENDIF
 
276
C
 
277
C *** get sigma, sateration, tolerance, max. number of iterations
 
278
       CALL STKRDR('INPUTR',3,4,IAC,RVAL,KUN,KNUL,ISTAT)    
 
279
       SIGMA1 = RVAL(1)
 
280
       SAT1   = RVAL(2)
 
281
       IF (SIGMA1.NE.0. .OR. SAT1.NE.0.) THEN
 
282
          SIGMA = SIGMA1
 
283
          SAT   = SAT1
 
284
       ENDIF
 
285
       PSA    = RVAL(3)
 
286
       PSAS   = RVAL(3)**2
 
287
       TNIT   = INT(RVAL(4))
 
288
       SIGG   = SIGMA
 
289
       ISX    = 0
 
290
       ISY    = 0
 
291
C
 
292
C *** get fitting method
 
293
       CALL STKRDC('INPUTC',1,1,20,IAC,FITMET,KUN,KNUL,ISTAT)       !fit method
 
294
       CALL UPCAS(FITMET,FITMET)
 
295
       IF (FITMET(1:2) .EQ.'MU') THEN
 
296
          MOF   = 'M'
 
297
          WE    = 'C'
 
298
          NCOMM = INDEX(FITMET,',')
 
299
          CBETA = FITMET(NCOMM+1:20)
 
300
          IF (NCOMM.GT.0) THEN
 
301
             CALL GENCNV(CBETA,2,1,IDUM,BETA,DDUM,ISTAT)
 
302
C            CALL USRINP(BETA,1,'R',CBETA)
 
303
             IF (ABS(BETA).LT.1.0E-30) THEN
 
304
                BETA = 4.0
 
305
             ENDIF
 
306
          ENDIF
 
307
 
 
308
       ELSE IF (FITMET(1:2).EQ.'MI') THEN
 
309
          MOF   = 'M'
 
310
          WE    = 'N'
 
311
          NCOMM = INDEX(FITMET,',')
 
312
          CBETA = FITMET(NCOMM+1:20)
 
313
          IF (NCOMM.GT.0) THEN
 
314
             CALL GENCNV(CBETA,2,1,IDUM,BETA,DDUM,ISTAT)
 
315
C            CALL USRINP(BETA,1,'R',CBETA)
 
316
             IF (ABS(BETA).LT.1.0E-30) THEN
 
317
                BETA = 4.0
 
318
             ENDIF
 
319
          ENDIF
 
320
 
 
321
       ELSE IF (FITMET(1:2).EQ.'GU') THEN
 
322
          MOF  = 'G'
 
323
          WE   = 'C'
 
324
          BETA = 0.0
 
325
 
 
326
       ELSE IF (FITMET(1:2).EQ.'GI') THEN
 
327
          MOF  = 'G'
 
328
          WE   = 'N'
 
329
          BETA = 0.0
 
330
 
 
331
       ELSE
 
332
          MOF  = 'M'
 
333
          WE   = 'C'
 
334
          BETA = 4.0
 
335
       ENDIF
 
336
C
 
337
C *** get fit option
 
338
       CALL STKRDC('INPUTC',1,21,4,IAC,FITOPT,KUN,KNUL,ISTAT)       !fit option
 
339
       CALL UPCAS(FITOPT,FITOPT)
 
340
       SIGFI = FITOPT(1:1)
 
341
       SKYFI = FITOPT(2:2)
 
342
       POSFI = FITOPT(3:3)
 
343
       TILT  = FITOPT(4:4)
 
344
C
 
345
C *** get the fix option, sigma, position, sky background, sky tilt
 
346
 1000  CONTINUE
 
347
       IF (SKYFI.EQ.'Y' .AND. POSFI.EQ.'Y') THEN
 
348
          CALL STTPUT('*** WARNING: Illegal combination, '//
 
349
     2                '... try again',ISTAT)
 
350
          FITOPT = 'YNNN'
 
351
          CALL STKWRC('INPUTC',1,FITOPT,21,4,KUN,ISTAT)
 
352
          CALL STKPRC(STRING,'INPUTC',1,21,4,IAC,FITOPT,KUN,KNUL,ISTAT) 
 
353
          CALL UPCAS(FITOPT,FITOPT)
 
354
          SIGFI = FITOPT(1:1)
 
355
          SKYFI = FITOPT(2:2)
 
356
          POSFI = FITOPT(3:3)
 
357
          TILT  = FITOPT(4:4)
 
358
          GO TO 1000
 
359
       END IF
 
360
C
 
361
       IF (SIGFI.NE.'N') SIGFI = 'Y'
 
362
       IF (POSFI.NE.'Y') POSFI = 'N'
 
363
       IF (SKYFI.NE.'Y') THEN
 
364
          SKYFI = 'N'
 
365
          IF (POSFI.NE.'Y') THEN
 
366
             TILT = FITOPT(4:4)
 
367
          END IF
 
368
          IF (TILT.NE.'Y') THEN
 
369
             TILT = 'N'
 
370
          ENDIF
 
371
       END IF
 
372
C
 
373
       CALL STKRDC('INPUTC',1,41,1,IAC,AVEOPT,KUN,KNUL,ISTAT)    ! aver. option
 
374
       CALL UPCAS(AVEOPT,AVEOPT)
 
375
       INME = AVEOPT
 
376
C
 
377
       CALL STKRDC('INPUTC',1,61,1,IAC,DIS,KUN,KNUL,ISTAT)        ! display
 
378
       CALL UPCAS(DIS,DIS)
 
379
       IF (DIS.EQ.'Y') THEN
 
380
          NDISPL = 0
 
381
       ELSE
 
382
          NDISPL = 1
 
383
       ENDIF
 
384
C
 
385
C ***  write the info to the screen
 
386
       WRITE(STRING,9001) FRAME
 
387
       CALL STTPUT(STRING,ISTAT)
 
388
       WRITE(STRING,9002) INTFIL
 
389
       CALL STTPUT(STRING,ISTAT)
 
390
       WRITE(STRING,9003) SIGMA, SAT       
 
391
       CALL STTPUT(STRING,ISTAT)
 
392
       WRITE(STRING,9004) PSA,TNIT
 
393
       CALL STTPUT(STRING,ISTAT)
 
394
C
 
395
       IF (FITMET(1:2).EQ.'MU') THEN
 
396
          WRITE(STRING,9005) BETA
 
397
       ELSE IF (FITMET(1:2).EQ.'MI') THEN
 
398
          WRITE(STRING,9006) BETA
 
399
       ELSE IF (FITMET(1:2).EQ.'GU') THEN
 
400
          STRING = 'Fit method: GAUSS, UNIFORM weighting'
 
401
       ELSE IF (FITMET(1:2).EQ.'GI') THEN
 
402
          STRING = 'Fit method: GAUSS, 1/N weighting'
 
403
       ENDIF
 
404
       CALL STTPUT(STRING,ISTAT)
 
405
C
 
406
       WRITE(STRING,9007) FITOPT
 
407
       CALL STTPUT(STRING,ISTAT)
 
408
C
 
409
       IF (AVEOPT.EQ.'S') THEN
 
410
          STRING = 'Trial values: mean sigma'
 
411
       ELSE IF (AVEOPT.EQ.'B') THEN
 
412
          STRING = 'Trial values: mean sky background'
 
413
       ELSE IF (AVEOPT.EQ.'A') THEN
 
414
          STRING = 'Trial values: mean sigma and mean sky background'
 
415
       ELSE
 
416
          STRING = 'Trial values: none'
 
417
       ENDIF
 
418
       CALL STTPUT(STRING,ISTAT)
 
419
C
 
420
C *** set the display flag
 
421
       CALL STKRDI('LOG',4,1,IAC,ODISPL,KUN,KNUL,ISTAT)
 
422
       CALL STKWRI('LOG',NDISPL,4,1,KUN,ISTAT)
 
423
       WRITE(STRING,9021)
 
424
       CALL STTPUT(STRING,ISTAT)
 
425
C *** do the work *****************************************************
 
426
C
 
427
       NNS  = 1
 
428
       NNF  = 1
 
429
       KAL  = 1
 
430
       ALT1 = 1.
 
431
       KSF  = 0.
 
432
       IOBJ = 0
 
433
       IROW = 1
 
434
C
 
435
 1001  CONTINUE
 
436
          CALL INTWRD(TIDINT,IROW,NCP,NHL)                     ! read the window
 
437
          CALL TBSGET(TIDINT,IROW,SFLAG,ISTAT)
 
438
          IF (SFLAG) THEN
 
439
             D1     = PARINT(1)
 
440
             D2     = PARINT(2)
 
441
             V      = PARINT(3)
 
442
             B      = PARINT(4)
 
443
             U      = PARINT(5)
 
444
             D3     = PARINT(6)
 
445
             D4     = PARINT(7)
 
446
             P(1)   = PARINT(8)
 
447
             P(2)   = PARINT(9)
 
448
             P(3)   = PARINT(10)
 
449
             BET    = PARINT(11)
 
450
             SCAP   = PARINT(12)
 
451
             D6     = PARINT(13)
 
452
             FL     = INT(PARINT(14))
 
453
C
 
454
             DO 1011 IS = 1,NCP
 
455
                P((IS-1)*4+4)  = FITCMP((IS-1)*6+1)
 
456
                P((IS-1)*4+5)  = FITCMP((IS-1)*6+2)
 
457
                P((IS-1)*4+6)  = FITCMP((IS-1)*6+3)
 
458
                P((IS-1)*4+7)  = FITCMP((IS-1)*6+4)
 
459
                USC(IS)        = FITCMP((IS-1)*6+5)
 
460
                USE(IS)        = FITCMP((IS-1)*6+6)
 
461
 1011        CONTINUE
 
462
C
 
463
             DO 1012 IH = 1,NHL
 
464
                PB((IH-1)*3+1) =  FITHOL((IH-1)*3+1)
 
465
                PB((IH-1)*3+2) =  FITHOL((IH-1)*3+2)
 
466
                PB((IH-1)*3+3) =  FITHOL((IH-1)*3+3)
 
467
 1012        CONTINUE
 
468
C
 
469
             DO 1013 K3 = 1,NTM                      ! initialize the flag array
 
470
                IF (FLGCMP(K3).GE.2) THEN
 
471
                   FLGCMP(K3) = 1
 
472
                ENDIF
 
473
                NONF(K3) = FLGCMP(K3)
 
474
                KFR(K3)  = 0
 
475
 1013        CONTINUE
 
476
C
 
477
             FL = 0
 
478
 
 
479
             IF (TILT.EQ.'N') THEN
 
480
                P(1) = 0.
 
481
                P(2) = 0.
 
482
             ENDIF   
 
483
             
 
484
             IF (FON1.GT.0 .OR. (INME.EQ.'B' .OR. INME.EQ.'A')) THEN
 
485
                IF (FON.NE.0) THEN
 
486
                   P(3) = FON
 
487
                ELSE
 
488
                   FON = P(3)
 
489
                END IF
 
490
                P(1) = 0.
 
491
                P(2) = 0.
 
492
             END IF
 
493
C
 
494
             D1  = D1+ISX
 
495
             D2  = D2+ISY
 
496
             IPX = D1
 
497
             IPY = D2
 
498
             NP  = D3
 
499
             NC  = D4
 
500
C
 
501
             ASSIGN 30011 TO PREVET                            ! PREPARA VETTORE
 
502
             GO TO 30001
 
503
30011        CONTINUE
 
504
C
 
505
             IF (NCOM.GT.0 .AND. FL.EQ.0) THEN
 
506
                IF (NHL.GT.0) THEN
 
507
                   ASSIGN 30012 TO PREMAS                     ! PREPARA MASCHERA
 
508
                   GO TO 30002
 
509
30012              CONTINUE
 
510
                END IF
 
511
C
 
512
                KP  = 0
 
513
                LJ1 = IPX-START(1)+1
 
514
                LJ3 = IPY-START(2)+1
 
515
                DO 1014 J5 = LJ3,LJ3+NC-1
 
516
                   JOLD = J5-LJ3+1
 
517
                   CALL REALIN(NPL,NL,J5,LJ1,NP,MADRID(IPNTR),RIA)
 
518
                   DO 1015 I = 1,NP
 
519
                      IF (MASK(JOLD,I).EQ.0) THEN
 
520
                         VAL = RIA(I) - FOG
 
521
                         IF (VAL.LT.SAT.AND.VAL.GT.0.) THEN
 
522
                            KP      = KP+1
 
523
                            IVX(KP) = I
 
524
                            IVY(KP) = JOLD
 
525
                            VZ(KP)  = VAL
 
526
                            IF (WE.EQ.'N') THEN
 
527
                               IF (VAL.GE.1) THEN
 
528
                                  WEI(KP) = 1./VAL
 
529
                               ENDIF
 
530
                            ELSE
 
531
                               WEI(KP) = 1
 
532
                            ENDIF
 
533
                         END IF
 
534
                      END IF
 
535
 1015              CONTINUE
 
536
 1014           CONTINUE
 
537
 
 
538
                IF (KP.GT.NCOM*4+3) THEN
 
539
                   ASSIGN 30013 TO FRASUZ                    !FRASER SUZUKI
 
540
                   GO TO 30003
 
541
30013              CONTINUE
 
542
 
 
543
                   IF (WFLAG.EQ.1) THEN
 
544
                      GR = 1
 
545
                   ELSE
 
546
                      GR = 2
 
547
                   END IF
 
548
                   WRITE (STRING,9012) IDNGRP,NCP,NITER
 
549
                   CALL STTPUT(STRING,ISTAT)
 
550
C
 
551
                   ASSIGN 30014 TO RIPVET                  !ripristina vettore
 
552
                   GO TO 30004
 
553
30014              CONTINUE
 
554
C
 
555
C                  IF (GR.EQ.1) THEN
 
556
                      DO 1016 K3 = 1,NTM
 
557
                         FLGCMP(K3) = NONF(K3)
 
558
 1016                 CONTINUE
 
559
C                  END IF
 
560
C
 
561
                   IF (AIN.GT.0) THEN
 
562
                      ASSIGN 30015 TO CALINT              !calcola interquartile
 
563
                      GO TO 30005
 
564
30015                 CONTINUE
 
565
                   END IF
 
566
 
 
567
                   DO 1017 KKL = 1,NCOM
 
568
                      IKKL = KTP+(KKL-1)*KTS
 
569
                      WRITE(STRING,9013) 
 
570
     2                   (SIGP(KL),KL=IKKL+1,IKKL+KTS)
 
571
                      IF (KKL.EQ.1) THEN
 
572
                         IF (KTP.GT.0) THEN
 
573
                            WRITE(STRING(53:),9014) 
 
574
     2                         (SIGP(KL),KL=1,KTP)
 
575
                         ENDIF
 
576
                      ENDIF
 
577
                      CALL STTPUT(STRING,ISTAT)
 
578
 1017              CONTINUE
 
579
                   CALL STTPUT(' ',ISTAT)
 
580
                      
 
581
                ELSE
 
582
                   GR   = 2
 
583
                   WRITE(STRING,9011)
 
584
                   CALL STTPUT(STRING,ISTAT)
 
585
                END IF
 
586
C
 
587
C *** write the output row for this window
 
588
                PARINT(1)  = D1
 
589
                PARINT(2)  = D2
 
590
                PARINT(3)  = V
 
591
                PARINT(4)  = B
 
592
                PARINT(5)  = U
 
593
                PARINT(6)  = FLOAT(NP)
 
594
                PARINT(7)  = FLOAT(NC)
 
595
                PARINT(8)  = P(1)
 
596
                PARINT(9)  = P(2)
 
597
                PARINT(10) = P(3)
 
598
                PARINT(11) = BETA
 
599
                PARINT(12) = SCAP
 
600
                PARINT(13) = FLOAT(NITER)
 
601
                PARINT(14) = FLOAT(GR)
 
602
C
 
603
                DO 1018 IS = 1,NCP
 
604
                   FITCMP((IS-1)*6+1) = P((IS-1)*4+4)
 
605
                   FITCMP((IS-1)*6+2) = P((IS-1)*4+5)
 
606
                   FITCMP((IS-1)*6+3) = P((IS-1)*4+6)
 
607
                   FITCMP((IS-1)*6+4) = P((IS-1)*4+7)
 
608
                   FITCMP((IS-1)*6+5) = USC(IS)
 
609
                   FITCMP((IS-1)*6+6) = USE(IS)
 
610
 1018           CONTINUE
 
611
C
 
612
                DO 1019 IH = 1,NHL
 
613
                   FITHOL((IH-1)*3+1) = PB((IH-1)*3+1)
 
614
                   FITHOL((IH-1)*3+3) = PB((IH-1)*3+2)
 
615
                   FITHOL((IH-1)*3+3) = PB((IH-1)*3+3)
 
616
 1019           CONTINUE
 
617
C
 
618
                CALL INTWWR(TIDINT,IROW,NCP,NHL)
 
619
                DO 1020 KRO = 1,NCP
 
620
                   RK4 = P(4*KRO)
 
621
                   IF (FLGCMP(KRO).EQ.0) THEN
 
622
                      JL0 = JL0+1
 
623
                      IF (RK4.GT.1) THEN
 
624
                         VEZE(JL0) = -2.5*ALOG10(RK4)
 
625
                      ELSE
 
626
                         VEZE(JL0) = 0
 
627
                      END IF
 
628
                      RMAY = AMAX1(RMAY,VEZE(JL0))
 
629
                      RMIY = AMIN1(RMIY,VEZE(JL0))
 
630
C
 
631
                   ELSE IF (FLGCMP(KRO).EQ.1) THEN
 
632
                      JL1 = JL1+1
 
633
                      IF (RK4.GT.1) THEN
 
634
                         VEUN(JL1) = -2.5*ALOG10(RK4)
 
635
                      ELSE
 
636
                         VEUN(JL1) = 0
 
637
                      END IF
 
638
                      RMAY = AMAX1(RMAY,VEUN(JL1))
 
639
                      RMIY = AMIN1(RMIY,VEUN(JL1))
 
640
C
 
641
                   ELSE IF (FLGCMP(KRO).GT.1)THEN
 
642
                      JL2 = JL2+1
 
643
                      IF (RK4.GT.1) THEN
 
644
                         VEMA(JL2) = -2.5*ALOG10(RK4)
 
645
                      ELSE
 
646
                         VEMA(JL2) = 0
 
647
                      END IF
 
648
                      RMAY = AMAX1(RMAY,VEMA(JL2))
 
649
                      RMIY = AMIN1(RMIY,VEMA(JL2))
 
650
                   END IF
 
651
 1020           CONTINUE
 
652
C
 
653
                CALL INTDWR(TIDINT,NGRP,NOBJ,NINT,SAT,FAT,SIGMA,
 
654
     2                                 BETA,SOFOT,AIN,FOG)
 
655
                IF (NHL.GT.0) THEN
 
656
                   ASSIGN 30016 TO AZZMAS                 !AZZERA MASCHERA
 
657
                   GO TO 30006
 
658
30016              CONTINUE
 
659
                END IF
 
660
             END IF
 
661
          END IF
 
662
          IROW = IROW + NCP + NHL
 
663
          IF (IROW.LE.NRINT) THEN
 
664
             GO TO 1001
 
665
          ENDIF
 
666
C
 
667
C *** display the histogram
 
668
       PAS  = .5
 
669
       MAXX = 0
 
670
       DO 2001 I = 1,JL0
 
671
          NCAL        = (RMAY-VEZE(I))/PAS+1
 
672
          IF (NCAL.GT.100) NCAL=100
 
673
          ISTZE(NCAL) = ISTZE(NCAL)+1
 
674
          MAXX        = MAX0(MAXX,ISTZE(NCAL))
 
675
 2001  CONTINUE
 
676
 
 
677
       DO 2002 I = 1,JL1
 
678
          NCAL        = (RMAY-VEUN(I))/PAS+1
 
679
          IF (NCAL.GT.100) NCAL=100
 
680
          ISTUN(NCAL) = ISTUN(NCAL)+1
 
681
          MAXX        = MAX0(MAXX,ISTUN(NCAL))
 
682
 2002  CONTINUE
 
683
 
 
684
       DO 2003 I = 1,JL2
 
685
          NCAL        = (RMAY-VEMA(I))/PAS+1
 
686
          IF (NCAL.GT.100) NCAL=100
 
687
          ISTMA(NCAL) = ISTMA(NCAL)+1
 
688
          MAXX        = MAX0(MAXX,ISTMA(NCAL))
 
689
 2003  CONTINUE
 
690
 
 
691
       NCM = (RMAY-RMIY)/PAS+1
 
692
       IF (JL0.GT.0) THEN
 
693
          WRITE(STRING,9021)
 
694
          CALL STTPUT(STRING,ISTAT)
 
695
          WRITE(STRING,9022) JL0
 
696
          CALL STTPUT(STRING,ISTAT)
 
697
          CALL STTPUT(' ',ISTAT)
 
698
          WRITE (STRING,9023)
 
699
          CALL STTPUT(STRING,ISTAT)
 
700
          CALL STTPUT(' ',ISTAT)
 
701
 
 
702
          DO 2004 I = 1,NCM
 
703
             NUX = 50.*ISTZE(I)/MAXX
 
704
             ANC = RMAY-PAS*(I-1)
 
705
             HH  = 10.**(-.4*ANC)
 
706
             IF (NUX.GT.0) THEN
 
707
                XXX = ' '
 
708
                DO 2005 K13 = 1,NUX
 
709
                   XXX(K13:K13) = 'X'
 
710
 2005           CONTINUE
 
711
                WRITE(STRING,9024) ISTZE(I),ANC,HH,XXX
 
712
                CALL STTPUT(STRING,ISTAT)
 
713
             ELSE
 
714
                WRITE(STRING,9025) ISTZE(I),ANC,HH
 
715
                CALL STTPUT(STRING,ISTAT)
 
716
             END IF
 
717
 2004     CONTINUE
 
718
       END IF
 
719
 
 
720
       IF (JL1.GT.0) THEN
 
721
          WRITE(STRING,9021)
 
722
          CALL STTPUT(STRING,ISTAT)
 
723
          WRITE(STRING,9026) JL1
 
724
          CALL STTPUT(STRING,ISTAT)
 
725
          CALL STTPUT(' ',ISTAT)
 
726
          DO 2006 I = 1,NCM
 
727
             NUX = 50.*ISTUN(I)/MAXX
 
728
             ANC = RMAY-PAS*(I-1)
 
729
             HH  = 10.**(-.4*ANC)
 
730
             IF (NUX.GT.0) THEN
 
731
                XXX = ' '
 
732
                DO 2007 K13 = 1,NUX
 
733
                   XXX(K13:K13) = 'X'
 
734
 2007           CONTINUE
 
735
                WRITE(STRING,9024) ISTUN(I),ANC,HH,XXX
 
736
                CALL STTPUT(STRING,ISTAT)
 
737
             ELSE
 
738
                WRITE(STRING,9025) ISTUN(I),ANC,HH
 
739
                CALL STTPUT(STRING,ISTAT)
 
740
             END IF
 
741
 2006     CONTINUE
 
742
       END IF
 
743
C
 
744
       IF (JL2.GT.0) THEN
 
745
          WRITE(STRING,9021)
 
746
          CALL STTPUT(STRING,ISTAT)
 
747
          WRITE(STRING,9027) JL2
 
748
          CALL STTPUT(STRING,ISTAT)
 
749
          CALL STTPUT(' ',ISTAT)
 
750
          DO 2008 I = 1,NCM
 
751
             NUX = 50.*ISTMA(I)/MAXX
 
752
             ANC = RMAY-PAS*(I-1)
 
753
             HH  = 10.**(-.4*ANC)
 
754
 
 
755
             IF (NUX.GT.0) THEN
 
756
                XXX = ' '
 
757
                DO 2009 K13 = 1,NUX
 
758
                   XXX(K13:K13) = 'X'
 
759
 2009           CONTINUE
 
760
                WRITE(STRING,9024) ISTMA(I),ANC,HH,XXX
 
761
                CALL STTPUT(STRING,ISTAT)
 
762
             ELSE
 
763
                WRITE(STRING,9025) ISTMA(I),ANC,HH
 
764
                CALL STTPUT(STRING,ISTAT)
 
765
             END IF
 
766
 2008     CONTINUE
 
767
       END IF
 
768
 
 
769
       CALL TBTCLO(TIDINT,ISTAT)
 
770
       CALL STKWRI('LOG',ODISPL,4,1,KUN,ISTAT)
 
771
       CALL STSEPI
 
772
 
 
773
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
774
C
 
775
C      Procedure PREVET
 
776
C
 
777
C-------------------------------------------------------------------------------
 
778
30001  CONTINUE
 
779
          DO 30101 LI = 1,3
 
780
             PA(LI) = P(LI)
 
781
30101     CONTINUE
 
782
          NCOM = 0
 
783
          DO 30201 LI = 1,NCP
 
784
             IF (FLGCMP(LI).GE.1) THEN
 
785
                IDE  = (LI-1)*4+3
 
786
                NCOM = NCOM+1
 
787
                IDI  = (NCOM-1)*4+3
 
788
C
 
789
                IF (INME.EQ.'S' .OR. INME.EQ.'A') THEN
 
790
                   P(IDE+4) = SIGG
 
791
                ELSE
 
792
                   P(IDE+4) = SIGMA
 
793
                ENDIF
 
794
C
 
795
                DO 30301 LLI = 1,4
 
796
                   PA(IDI+LLI) = P(IDE+LLI)
 
797
30301           CONTINUE
 
798
                PA(IDI+1)  = PA(IDI+1)*ALT1
 
799
                FLGCMP(LI) = 1
 
800
                NONF(LI)   = 1
 
801
             END IF
 
802
30201     CONTINUE
 
803
       GO TO PREVET
 
804
 
 
805
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
806
C
 
807
C      Procedure PREMAS
 
808
C
 
809
C-------------------------------------------------------------------------------
 
810
30002  CONTINUE                                            !PREPARA MASCHERA
 
811
          DO 30102 LI = 1,NHL
 
812
             ICB = (LI-1)*3
 
813
             RR  = PB(ICB+1)**2
 
814
             L4  = AMAX1(1.,PB(ICB+3)-PB(ICB+1))
 
815
             L5  = AMIN1(FLOAT(NC),PB(ICB+3)+PB(ICB+1)+.99)
 
816
             L6  = AMAX1(1.,PB(ICB+2)-PB(ICB+1))
 
817
             L7  = AMIN1(FLOAT(NP),PB(ICB+2)+PB(ICB+1)+.99)
 
818
             DO 30202 L2 = L6,L7
 
819
                DO 30302 L3 = L4,L5
 
820
                   DELT = (FLOAT(L2)-PB(ICB+2))**2 + 
 
821
     2                    (FLOAT(L3)-PB(ICB+3))**2
 
822
                   IF (RR.GE.DELT) THEN
 
823
                      MASK(L3,L2) = 1
 
824
                   END IF
 
825
30302           CONTINUE
 
826
30202        CONTINUE
 
827
30102     CONTINUE
 
828
       GO TO PREMAS
 
829
 
 
830
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
831
C
 
832
C      Procedure FRASUZ
 
833
C
 
834
C-------------------------------------------------------------------------------
 
835
30003  CONTINUE                                            ! TO FRASER SUZUKI
 
836
          WFLAG = 1
 
837
          LFLAG = 0
 
838
          SLU   = 0.00001
 
839
          NITER = 0
 
840
          SCAP  = 10.**15
 
841
          US    = 10.**10
 
842
          KPR1  = 0
 
843
          NRFR  = 0
 
844
          IF (POSFI.EQ.'Y' .AND. SIGFI.EQ.'Y') THEN
 
845
             CALL ELMRPF(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA,
 
846
     2                   SQM,LFLAG,WEI,SIGP)
 
847
             KTP = 1
 
848
             KTS = 1
 
849
             DO 30103 K9 = 1,NCOM
 
850
                NONF(K9) = FLGCMP(K9)
 
851
30103        CONTINUE
 
852
             IF (LFLAG.EQ.0) THEN
 
853
                LFLAG = 1
 
854
             ENDIF
 
855
             NITER = 1
 
856
 
 
857
          ELSE
 
858
39803        CONTINUE
 
859
             IF (LFLAG.NE.0) GO TO 39903
 
860
39703           CONTINUE
 
861
                RTELOG = ABS(US).LE.SLU .OR. NITER.GE.TNIT .OR.
 
862
     2                   LFLAG.NE.0 
 
863
                IF (RTELOG) GO TO 39603
 
864
                   NITER = NITER+1
 
865
                   IF (MOF.EQ.'M') THEN
 
866
                      IF (NITER.LE.50.) THEN
 
867
                         D = 0.7
 
868
                      ELSE IF (NITER.GT.50. .AND. NITER.LE.80.) THEN
 
869
                         D = 0.3
 
870
                      ELSE
 
871
                         D = 0.7
 
872
                      END IF
 
873
                   END IF
 
874
                   PES(4) = 2
 
875
                   IF (NITER.LE.1.) PES(4)=.1
 
876
                   LG = 0
 
877
C
 
878
                   IF (SIGFI.EQ.'Y' .AND. SKYFI.EQ.'Y' .AND. 
 
879
     2                POSFI.NE.'Y') THEN
 
880
                      CALL ELMRF(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA,
 
881
     2                           SQM,LFLAG,WEI,SIGP)
 
882
                      KTP = 0
 
883
                      KTS = 3
 
884
 
 
885
                   ELSE IF (SIGFI.EQ.'N' .AND. SKYFI.EQ.'Y' .AND.
 
886
     2                      POSFI.NE.'Y') THEN
 
887
                      CALL ELMRFV(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA,
 
888
     2                            SQM,LFLAG,WEI,SIGP)
 
889
                      KTP = 0
 
890
                      KTS = 4
 
891
 
 
892
                   ELSE IF (SIGFI.EQ.'Y' .AND. TILT.EQ.'N' .AND.
 
893
     2                      POSFI.NE.'Y') THEN
 
894
                      CALL ELMRR(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA,
 
895
     2                           SQM,LFLAG,WEI,SIGP)
 
896
                      KTP = 1
 
897
                      KTS = 3
 
898
 
 
899
                   ELSE IF (SIGFI.EQ.'Y' .AND. TILT.EQ.'Y' .AND. 
 
900
     2                      POSFI.NE.'Y') THEN
 
901
                      CALL ELMRX(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA,
 
902
     2                            SQM,LFLAG,WEI,SIGP)
 
903
                      KTP = 3
 
904
                      KTS = 3
 
905
 
 
906
                   ELSE IF (SIGFI.EQ.'N' .AND. TILT.EQ.'N' .AND. 
 
907
     2                      POSFI.NE.'Y') THEN
 
908
                      CALL ELMRRV(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA,
 
909
     2                            SQM,LFLAG,WEI,SIGP)
 
910
                      KTP = 1
 
911
                      KTS = 4
 
912
 
 
913
                   ELSE IF (SIGFI.EQ.'N' .AND. TILT.EQ.'Y' .AND.
 
914
     2                      POSFI.NE.'Y') THEN
 
915
                      CALL ELMRV(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA,
 
916
     2                           SQM,LFLAG,WEI,SIGP)
 
917
                      KTP = 3
 
918
                      KTS = 4
 
919
 
 
920
                   ELSE IF(SIGFI.EQ.'N' .AND. POSFI.EQ.'Y') THEN
 
921
                      CALL ELMRPV(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA,
 
922
     2                            SQM,LFLAG,WEI,SIGP)
 
923
                      KTP = 1
 
924
                      KTS = 2
 
925
                   END IF
 
926
C
 
927
                   IF (LFLAG.EQ.0) THEN
 
928
                      US   = (SCAP-SQM)/SQM
 
929
                      SCAP = SQM
 
930
                   END IF
 
931
C
 
932
                   ASSIGN 30017 TO CONPAR
 
933
                   GO TO 30007
 
934
30017              CONTINUE
 
935
                   GO TO 39703
 
936
39603           CONTINUE
 
937
C
 
938
                IF (NRFR.GT.0 .AND. LFLAG.EQ.0) THEN
 
939
                   NPSG  = 0
 
940
                   HGU   = 0
 
941
                   DO 30203 K3 = 1,NCP
 
942
                      IF (KFR(K3).GT.0) THEN
 
943
                         IND = K3*4
 
944
                         IF (P(IND).GT.HGU) THEN
 
945
                            HGU  = P(IND)
 
946
                            NPSG = K3
 
947
                         END IF
 
948
                      END IF
 
949
30203              CONTINUE
 
950
 
 
951
                   NRFR=NRFR-1
 
952
C ***                                                          !calcola code
 
953
                   IND=(NPSG-1)*4+5
 
954
                   DO 30303 K7=1,NCOM
 
955
                      K6=K7*4
 
956
                      ABC = ((PA(K6+1)-P(IND))**2 +
 
957
     2                      (PA(K6+2)-P(IND+1))**2)/PA(K6+3)**2
 
958
                      IF (ABS(BETA).LT.1.0E-30) THEN
 
959
                         ABC = PA(K6)*EXP(-ABC*4*ALOG(2.))
 
960
                      ELSE
 
961
                         ABC = PA(K6)*(1.+ABC)**(-BETA)
 
962
                      END IF
 
963
                      HGU = HGU-ABC
 
964
30303              CONTINUE
 
965
                   HGU = HGU+P(3)-PA(3)
 
966
                   WRITE (STRING,9031) NITER,NPSG
 
967
                   CALL STTPUT(STRING,ISTAT)
 
968
C
 
969
                   IF (HGU.GT.SOFOT) THEN
 
970
                      FLGCMP(NPSG) = 1
 
971
                      SLU          = 0.00001
 
972
                      NITER        = 0
 
973
                      SCAP         = 10.**15
 
974
                      US           = 10.**10
 
975
                      IND          = 1
 
976
                      IND2         = 0
 
977
C
 
978
                      IF (NPSG.LE.IND) GO TO 30503
 
979
                         IF (FLGCMP(IND).GE.1) IND2=IND2+1
 
980
                         IND = IND+1
 
981
30503                 CONTINUE
 
982
C
 
983
                      IF (IND2.LT.NCOM) THEN
 
984
                         DO 30603 LI = NCOM,IND2+1,-1
 
985
                            IND = LI*4
 
986
                            DO 30703 K3 = IND,IND+3
 
987
                               PA(K3+4) = PA(K3)
 
988
30703                       CONTINUE
 
989
30603                    CONTINUE
 
990
                      END IF
 
991
                      NCOM = NCOM+1
 
992
                      IND  = (IND2+1)*4
 
993
                      IND1 = NPSG*4
 
994
                      DO 30803 K3 = IND,IND+3
 
995
                         PA(K3) = P(IND1+K3-IND)
 
996
30803                 CONTINUE
 
997
                   ELSE
 
998
                      WRITE (STRING,9032)
 
999
                      CALL STTPUT(STRING,ISTAT)
 
1000
                   END IF
 
1001
                   KFR(NPSG)=-KFR(NPSG)
 
1002
                ELSE
 
1003
                   IF (LFLAG.EQ.0) THEN
 
1004
                      LFLAG = 1
 
1005
                   ENDIF
 
1006
                ENDIF
 
1007
                GO TO 39803
 
1008
39903        CONTINUE
 
1009
          END IF
 
1010
       GO TO FRASUZ
 
1011
 
 
1012
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
1013
C
 
1014
C      Procedure RIPVET
 
1015
C
 
1016
C-------------------------------------------------------------------------------
 
1017
30004  CONTINUE                                           !RIPRISTINA VETTORE
 
1018
          DO 30104 LI=1,3
 
1019
             P(LI)=PA(LI)
 
1020
30104     CONTINUE
 
1021
          IF (GR.NE.2) THEN
 
1022
             IF (INME.EQ.'B'.OR.INME.EQ.'A') THEN
 
1023
                RKN = 1./NNF
 
1024
                NNF = NNF+1
 
1025
                FON = (FON+PA(3)*RKN)/(1+RKN)
 
1026
             END IF
 
1027
          END IF
 
1028
          NCO=0
 
1029
          DO 30204 LI = 1,NCP
 
1030
             IF (NONF(LI).GE.1) THEN
 
1031
                IDE = (LI-1)*4+3
 
1032
                NCO = NCO+1
 
1033
                IDI = (NCO-1)*4+3
 
1034
                IF (NONF(LI).EQ.1) THEN
 
1035
                   KAS  = 1./KAL
 
1036
                   KAL  = KAL+1
 
1037
                   ALT1 = (ALT1+PA(IDI+1)/P(IDE+1)*KAS)/(1.+KAS)
 
1038
                   DO 30304 LLI = 1,4
 
1039
                      P(IDE+LLI) = PA(IDI+LLI)
 
1040
30304              CONTINUE
 
1041
                   IF (SIGFI.EQ.'N' .AND.
 
1042
     2                (INME.EQ.'S'.OR.INME.EQ.'A')) THEN
 
1043
                      RKN  = 1./NNS
 
1044
                      NNS  = NNS+1
 
1045
                      SIGG = (SIGG+PA(IDI+4)*RKN)/(1+RKN)
 
1046
                   END IF
 
1047
                END IF
 
1048
             END IF
 
1049
30204     CONTINUE
 
1050
       GO TO RIPVET
 
1051
 
 
1052
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
1053
C
 
1054
C      Procedure CALINT
 
1055
C
 
1056
C-------------------------------------------------------------------------------
 
1057
30005  CONTINUE                                        ! CALCOLA INTERQUARTILE
 
1058
          IF (BETA.GT.0.0) BEBE=-1./BETA
 
1059
C
 
1060
          DO 30105 JIL = 1,NCP
 
1061
             JIL4     = JIL*4 
 
1062
             USC(JIL) = 0.
 
1063
             USE(JIL) = 0.
 
1064
             IF (FLGCMP(JIL).EQ.1) THEN
 
1065
                SIGQ = P(JIL4+3)**2
 
1066
                IF (ABS(BETA).LT.1.0E-30) THEN
 
1067
                   ALAL = 4*ALOG(2.)
 
1068
                   SMT  = AIN*EXP((1./SIGQ)*ALAL)
 
1069
                   IF (P(JIL4).GE.SMT) THEN
 
1070
                      RAG(JIL) = P(JIL4+3)**2*
 
1071
     2                             (-ALOG(AIN/P(JIL4))/ALAL)
 
1072
                   ELSE
 
1073
                      RAG(JIL) = 1.0
 
1074
                   END IF
 
1075
                ELSE
 
1076
                   BEBE = -1./BETA
 
1077
                   SMT  = AIN*((1./SIGQ)+1.)**BETA
 
1078
                   IF (P(JIL4).GE.SMT) THEN
 
1079
                      RAG(JIL) = P(JIL4+3)**2*
 
1080
     2                           ((AIN/P(JIL4))**BEBE-1.)
 
1081
                   ELSE
 
1082
                      RAG(JIL) = 1.0
 
1083
                   END IF
 
1084
                END IF
 
1085
             END IF
 
1086
30105     CONTINUE
 
1087
C
 
1088
          DO 30205 K7=1,NCP
 
1089
             K74 = K7*4
 
1090
             IF (FLGCMP(K7).EQ.1) THEN
 
1091
                SQMS = 0.
 
1092
                NPS  = 0
 
1093
                DO 30305 IQ = 1,NCA
 
1094
                   RM(IQ) = 0
 
1095
30305           CONTINUE
 
1096
                DO 30405 IQ=1,KP
 
1097
                   DPS = (IVX(IQ)-P(K74+1))**2 +
 
1098
     2                   (IVY(IQ)-P(K74+2))**2
 
1099
                   IF (DPS.LE.RAG(K7) .AND. VZ(IQ).LT.SAT) THEN
 
1100
                      NPS = NPS+1
 
1101
                      VG  = P(3)
 
1102
                      DO 30505 JQ = 1,NCP
 
1103
                         JQ4 = JQ*4
 
1104
                         IF (FLGCMP(JQ).EQ.1) THEN
 
1105
                            GIG = P(JQ4+3)**2.
 
1106
                            EEE = ((P(JQ4+1)-IVX(IQ))**2 +
 
1107
     2                             (P(JQ4+2) - IVY(IQ))**2)/GIG
 
1108
                            IF (ABS(BETA).LT.1.0E-30) THEN
 
1109
                               VG = VG+P(JQ4)*EXP(-EEE*ALAL)
 
1110
                            ELSE
 
1111
                               VG = VG+P(JQ4)*(1.+EEE)**(-BETA)
 
1112
                            END IF
 
1113
                         END IF
 
1114
30505                 CONTINUE
 
1115
C
 
1116
                      SQMS       = SQMS+(VZ(IQ)-VG)**2*WEI(IQ)
 
1117
                      VALME(NPS) = (VZ(IQ)-VG)/SQRT(VG)
 
1118
                   ENDIF
 
1119
30405           CONTINUE
 
1120
C
 
1121
                DO 30705 IQ = 2,NPS
 
1122
                   A = VALME(IQ)
 
1123
                   DO 30605 JS = IQ-1,1,-1
 
1124
                      IF (VALME(JS) .LE. A) THEN 
 
1125
                         GO TO 31605
 
1126
                      ELSE
 
1127
                         VALME(JS+1) = VALME(JS)
 
1128
                      ENDIF
 
1129
30605              CONTINUE
 
1130
                   JS = 0
 
1131
31605              CONTINUE
 
1132
                   VALME(JS+1) = A
 
1133
30705           CONTINUE
 
1134
C
 
1135
                IRINT = FLOAT(NPS)/4. + 0.5
 
1136
                RINTQ = VALME(IRINT*3)-VALME(IRINT)
 
1137
C
 
1138
                IF (NPS.LE.4) THEN
 
1139
                   NPS = 5
 
1140
                ENDIF
 
1141
                SQMS = SQMS/FLOAT(NPS-4)
 
1142
                IF (WE.NE.'N') THEN
 
1143
                   SQMS = SQMS/P(K74)
 
1144
                ENDIF
 
1145
                IF (RINTQ.LT.0) THEN
 
1146
                   RINTQ = 50
 
1147
                ENDIF
 
1148
                USC(K7) = SQMS
 
1149
                USE(K7) = RINTQ
 
1150
             END IF
 
1151
30205     CONTINUE
 
1152
       GO TO CALINT
 
1153
 
 
1154
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
1155
C
 
1156
C      Procedure AZZMA
 
1157
C
 
1158
C-------------------------------------------------------------------------------
 
1159
30006  CONTINUE                                           !AZZERA MASCHERA
 
1160
          DO 30106 LI=1,NHL
 
1161
             ICB = (LI-1)*3
 
1162
             L4  = AMAX1(1.,PB(ICB+3)-PB(ICB+1))
 
1163
             L5  = AMIN1(FLOAT(NC),PB(ICB+3)+PB(ICB+1)+.99)
 
1164
             L6  = AMAX1(1.,PB(ICB+2)-PB(ICB+1))
 
1165
             L7  = AMIN1(FLOAT(NP),PB(ICB+2)+PB(ICB+1)+.99)
 
1166
             DO 30206 L2 = L6,L7
 
1167
                DO 30306 L3 = L4,L5
 
1168
                   MASK(L3,L2) = 0
 
1169
30306           CONTINUE
 
1170
30206        CONTINUE
 
1171
30106     CONTINUE
 
1172
       GO TO AZZMAS
 
1173
 
 
1174
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
1175
C
 
1176
C      Procedure CONPAR
 
1177
C
 
1178
C-------------------------------------------------------------------------------
 
1179
30007  CONTINUE                                          !CONTROLLA PARAMETRI
 
1180
          IF(LFLAG.EQ.0) THEN
 
1181
             KG   = 0
 
1182
             IP   = 0
 
1183
             IPA  = 0
 
1184
             LFL2 = 0
 
1185
             DO 30107 K3 = 1,NCOM
 
1186
                IF (LFL2.EQ.0) THEN
 
1187
                   IPA       = IPA+1
 
1188
                   NONF(IPA) = FLGCMP(IPA)
 
1189
30207              CONTINUE
 
1190
                   IF (FLGCMP(IPA).NE.0 .OR. IPA.GE.NCP) THEN
 
1191
                      GOTO 30307
 
1192
                   ENDIF
 
1193
                   IPA       = IPA+1
 
1194
                   NONF(IPA) = FLGCMP(IPA)
 
1195
                   GOTO 30207
 
1196
30307              CONTINUE
 
1197
                   INK1 = (IPA-1)*4
 
1198
                   INK  = (K3-1)*4
 
1199
                   IF (PA(INK+4).LT.0) THEN
 
1200
                      WFLAG     = 2
 
1201
                      LFLAG     = 2
 
1202
                      NONF(IPA) = 3                      ! Altezza neagtive
 
1203
                   END IF
 
1204
                   IF (PA(INK+7).LT.0. .OR. PA(INK+7).GT.100.) THEN
 
1205
                      WFLAG     = 2
 
1206
                      LFLAG     = 2
 
1207
                      NONF(IPA) = 5                       !Sigma errato
 
1208
                   END IF
 
1209
                   IF (PA(INK+4).GT.1.E10) THEN
 
1210
                      WFLAG     = 2
 
1211
                      LFLAG     = 2
 
1212
                      NONF(IPA) = 3               !Alt. magg. di 10.000.000
 
1213
                   END IF
 
1214
                   DXY2 = (PA(INK+5) - P(INK1+5))**2 + 
 
1215
     2                    (PA(INK+6) - P(INK1+6))**2
 
1216
                   IF (DXY2.GT.PSAS) THEN
 
1217
                      WFLAG     = 2
 
1218
                      LFLAG     = 2
 
1219
                      NONF(IPA) = 4              !La posizione del massimo
 
1220
                   END IF                    !si e' spostata piu' di PSA pixels
 
1221
 
 
1222
                   IF (LFLAG.EQ.2) THEN
 
1223
                      IF (KFR(IPA).EQ.0) THEN
 
1224
                         NRFR        = NRFR+1
 
1225
                         KFR(IPA)    = K3
 
1226
                         FLGCMP(IPA) = 0
 
1227
                         LFLAG       = 0
 
1228
                         KPR1        = 1
 
1229
                      ELSE
 
1230
                         LFL2        = 1
 
1231
                      END IF
 
1232
                   END IF
 
1233
                END IF
 
1234
30107        CONTINUE
 
1235
C
 
1236
             IF (KPR1.EQ.1 .AND. LFL2.EQ.0) THEN
 
1237
                IF (NRFR.GT.1) THEN
 
1238
                   WRITE (STRING,9035) NRFR
 
1239
                   CALL STTPUT(STRING,ISTAT)
 
1240
                ELSE
 
1241
                   STRING = 'No convergence on 1 star'
 
1242
                   CALL STTPUT(STRING,ISTAT)
 
1243
                ENDIF
 
1244
C
 
1245
                KPR1  = 0
 
1246
                SLU   = 0.00001
 
1247
                NITER = 0
 
1248
                SCAP  = 10.**15
 
1249
                US    = 10.**10
 
1250
C
 
1251
                ASSIGN 30018 TO PREVET                        ! #PREPARA VETTORE
 
1252
                GO TO 30001
 
1253
30018           CONTINUE
 
1254
 
 
1255
                IF (NCOM.LE.0) THEN
 
1256
                   LFLAG = 2
 
1257
                ENDIF
 
1258
             ENDIF
 
1259
          ELSE
 
1260
             DO 30407 K3 = 1,NCOM
 
1261
                FLGCMP(K3) = 6
 
1262
                NONF(K3)   = 6
 
1263
30407        CONTINUE 
 
1264
          ENDIF
 
1265
C
 
1266
          IF (LFL2.EQ.1) THEN
 
1267
             WFLAG = 2
 
1268
             LFLAG = 2
 
1269
             WRITE(STRING,9037)
 
1270
             CALL STTPUT(STRING,ISTAT)
 
1271
             DO 30408 K3 = 1,NCOM
 
1272
                INK = (K3-1)*4
 
1273
                WRITE(STRING,9036) K3,(PA(INK+JIL),JIL=4,7)
 
1274
                CALL STTPUT(STRING,ISTAT)
 
1275
30408        CONTINUE
 
1276
          END IF
 
1277
       GO TO CONPAR
 
1278
       END
 
1279