C @(#)rfotrfit.for 17.1.1.1 (ES0-DMD) 01/25/02 17:18:17 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C PROGRAM FIT C--- C.IDENTIFICATION: RFOTFIT C.PURPOSE: Determine characteristics of object by non-linear fitting C. Data and results are taken from and written the int. table C.AUTHOR: R. Buonanno, G. Buscema, C. Corsi, I. Ferraro, G. Iannicola C. Osservatorio Astronomico di Roma C.VERSION: 05/14/87 date of creation C.VERSION: 15/12/87 RHW implementation in MIDAS standard interfaces C.VERSION: 10/03/88 RB new version C.VERSION: 03/05/89 RHW change to ST interfaces and streamlining code C.VERSION: 22/01/91 RHW all variables declared; implicit none added C+++ IMPLICIT NONE INCLUDE 'MID_REL_INCL:RFOTDECL.INC' C INTEGER NDXY INTEGER NCA INTEGER MBB INTEGER MPA INTEGER NCD INTEGER NTM INTEGER NOB INTEGER M12 INTEGER IDP INTEGER IDB PARAMETER (NDXY=55) PARAMETER (NCA=200) PARAMETER (MBB=256) PARAMETER (MPA=10000) PARAMETER (NCD=MBB-18) PARAMETER (NTM=NCD/6.5) PARAMETER (NOB=NCD/3) PARAMETER (M12=NDXY*NDXY) PARAMETER (IDP=(NTM*4)+3) PARAMETER (IDB=NOB*3) C INTEGER ODISPL,NDISPL INTEGER EC, ED, EL INTEGER FL, GR INTEGER I, IS, IH, IQ INTEGER IPX, IPY INTEGER IOBJ, IROW, IDUM INTEGER IKKL, IFLG INTEGER*8 IPNTR INTEGER IMF INTEGER IP, IPA INTEGER IND, INK, INK1 INTEGER ISTZE(100),ISTUN(100),ISTMA(100) INTEGER IVX(M12),IVY(M12) INTEGER ISC(7), ISX, ISY INTEGER ISTAT, IAC INTEGER IDI, IND1, IND2 INTEGER IRINT INTEGER IDE, ICB INTEGER JL0, JL1, JL2, JIL4 INTEGER J5, JIL, JOLD, JQ4 INTEGER JQ, JS INTEGER KFR(NTM) INTEGER KUN, KNUL INTEGER K3, KP, KG, KKL, KL, KTP, KTS, KRO, K13 INTEGER KAL, KSF, KPR1, KAS INTEGER K6, K7, K9, K74 INTEGER MADRID(1) INTEGER MAXX INTEGER LI, LLI INTEGER L2, L3, L4, L5, L6, L7 INTEGER LJ1, LJ3 INTEGER LG, LFL2 INTEGER LFLAG INTEGER MASK(NDXY,NDXY) INTEGER NCINT,NRINT,NSINT,NACINT,NARINT INTEGER NOBJ, NGRP, NINT, NAXIS INTEGER NPIX(3) INTEGER NONF(NTM) INTEGER NPL, NL, NPS INTEGER NNS, NNF INTEGER NITER, NC, NP, NCOM, NCO INTEGER NCM, NUX, NRFR INTEGER NHL, NCP, NCOMM, NCAL INTEGER NPSG INTEGER RM(NCA) INTEGER TNIT, TINULL INTEGER TIDINT INTEGER WFLAG C INTEGER PREVET,PREMAS,FRASUZ,RIPVET,CALINT,AZZMAS,CONPAR C DOUBLE PRECISION START(3),STEP(3) DOUBLE PRECISION TDNULL,TDTRUE,TDFALS,DDUM C REAL A, ALAL, ABC, ALT1 REAL AIN, ANC REAL B, BET,BETA, BEBE REAL D, D1, D2, D3, D4, D6 REAL DELT, DPS, DXY2 REAL EEE REAL FOG, FON, FON1, FOG1, FAT REAL GIG REAL HH, HGU REAL P(IDP),PES(7),SIGP(1000),RIA(NDXY) REAL PSA, PSAS, PAS REAL PA(IDP),PB(IDB) REAL RAG(NTM) REAL RVAL(5) REAL RK4, RMIY, RKN, RR, RINTQ, RMAY REAL SAT, SIGMA, SOFOT, SIGG, SCAP, SAT1 REAL SLU, SQM , SIGQ, SQMS, SMT, SIGMA1 REAL TRNULL, TBLSEL REAL U, US REAL USC(NTM),USE(NTM) REAL V, VAL, VG REAL VALME(M12) REAL VEZE(MPA),VEUN(MPA),VEMA(MPA) REAL VZ(M12) REAL WEI(M12) C CHARACTER FITMET*20,CBETA*80,FITOPT*4,MOF*1,WE*1,AVEOPT*1 CHARACTER DIS*1 CHARACTER FRAME*60,IDENT*72,CUNIT*72 CHARACTER INTFIL*60 CHARACTER*80 STRING CHARACTER*60 XXX CHARACTER INME*1 CHARACTER SIGFI*1,SKYFI*1,POSFI*1,IUT*1,TILT*1 C LOGICAL RTELOG LOGICAL SFLAG C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID INCLUDE 'MID_INCLUDE:ST_DAT.INC' DATA PES/1.,1.,2.,2.,1.,1.,2./ DATA D/0.7/ DATA ISC/7*1/ DATA ISTZE/100*0/ DATA ISTUN/100*0/ DATA ISTMA/100*0/ 9001 FORMAT('Frame: ',A40) 9002 FORMAT('Intermediate table: ',A40) 9003 FORMAT('Sigma: ', G13.6,'; Saturation: ', G13.6) 9004 FORMAT('Tolerance: ', G13.6,'; # Iterations: ', I5) 9005 FORMAT('Fit method: MOFFET, UNIFORM weighting; BETA: ',G13.6) 9006 FORMAT('Fit method: MOFFET, 1/N weighting; BETA: ', G13.6) 9007 FORMAT('Fix sigma, fix position, fix sky, tilt sky plane: ',A) C 9011 FORMAT('Too few pixel in this window') 9012 FORMAT('Window identification: ',I6,2X, 2 ' Components: ',I2,2X,' Iterations: ',I2) 9013 FORMAT('Fit errors: ',4F9.4) 9014 FORMAT(3F9.4) C 9021 FORMAT(80('-')) 9022 FORMAT('Components not examined: ',I6) 9023 FORMAT('Total Mag. index Histogram') 9024 FORMAT(1X,I5,2X,F5.1,1X,F11.0,1X,'I',A50) 9025 FORMAT(1X,I5,2X,F5.1,1X,F11.0,1X,'I') 9026 FORMAT('Components succesfully fitted: ',I6) 9027 FORMAT('Components with no convergency: ',I6) 9031 FORMAT('Iteration: ',I6,'; Added star: ',I6) 9032 FORMAT('Component below the photometric threshold') 9035 FORMAT('N O C O N V E R G E N C E on ',I3,' components') 9036 FORMAT(1X,I6,' Component Parameters: ',F15.1,2F9.1,F7.2) 9037 FORMAT('N O C O N V E R G E N C E') C C *** Start the program ***************************************************** CALL STSPRO('FIT') CALL TBMNUL(TINULL,TRNULL,TDNULL) CALL TBMCON(TBLSEL,TDTRUE,TDFALS) JL0 = 0 JL1 = 0 JL2 = 0 RMAY = -10.**15 RMIY = -RMAY C C *** read the frame CALL STKRDC('IN_A',1,1,60,IAC,FRAME,KUN,KNUL,ISTAT) CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,3,NAXIS, 2 NPIX,START,STEP,IDENT,CUNIT,IPNTR,IMF,ISTAT) NPL = NPIX(1) NL = NPIX(2) C C *** read the intermediate file CALL STKRDC('IN_B',1,1,60,IAC,STRING,KUN,KNUL,ISTAT) ! intermediate file IFLG = INDEX(INTFIL,'/') ! look for the flag IF (IFLG.EQ.0) THEN INTFIL = STRING ELSE IUT = 'F' INTFIL = STRING(1:IFLG-1) ENDIF CALL STECNT('GET',EC,ED,EL) CALL STECNT('PUT',1,0,0) CALL TBTOPN(INTFIL,F_IO_MODE,TIDINT,ISTAT) IF (ISTAT.NE.0) THEN STRING = '*** FATAL: Problems with opening intermediate'// 2 ' table ... ' CALL STTPUT(STRING,ISTAT) CALL STSEPI ENDIF C C *** get table info CALL TBIGET(TIDINT,NCINT,NRINT,NSINT,NACINT,NARINT,ISTAT) IF (ISTAT.NE.0) THEN STRING = '*** FATAL: Problems with getting info for '// 2 ' intermediate table; Try again ... ' CALL STTPUT(STRING,ISTAT) CALL STSEPI ENDIF IF (NRINT.EQ.0) THEN STRING = '*** FATAL: No data points in intermediate table' CALL STTPUT(STRING,ISTAT) CALL STSEPI ENDIF CALL STECNT('PUT',EC,ED,EL) C C *** read table descriptor CALL INTDRD(TIDINT,NGRP,NOBJ,NINT, 2 SAT,FAT,SIGMA,BETA,SOFOT,AIN,FOG) C C *** get threshold, sky background CALL STKRDR('INPUTR',1,2,IAC,RVAL,KUN,KNUL,ISTAT) !thresh., sky backgr IF (RVAL(1) .GT. 0.0) THEN SOFOT = RVAL(1) ENDIF C FON1 = RVAL(2) FOG1 = 0.0 IF (FON1.NE.0.0) THEN FON = FON1 FOG = FOG1 ENDIF C C FON1 = RVAL(2) C ALFA = 200.0 C FOG1 = 0.0 C IF (FOG1.NE.0.0 .OR. FON1.NE.0.0 .OR. A1.NE.0.0) THEN C A1 = 200. C FON = FON1 C ALFA = A1 C FOG = FOG1 C ENDIF C C *** get sigma, sateration, tolerance, max. number of iterations CALL STKRDR('INPUTR',3,4,IAC,RVAL,KUN,KNUL,ISTAT) SIGMA1 = RVAL(1) SAT1 = RVAL(2) IF (SIGMA1.NE.0. .OR. SAT1.NE.0.) THEN SIGMA = SIGMA1 SAT = SAT1 ENDIF PSA = RVAL(3) PSAS = RVAL(3)**2 TNIT = INT(RVAL(4)) SIGG = SIGMA ISX = 0 ISY = 0 C C *** get fitting method CALL STKRDC('INPUTC',1,1,20,IAC,FITMET,KUN,KNUL,ISTAT) !fit method CALL UPCAS(FITMET,FITMET) IF (FITMET(1:2) .EQ.'MU') THEN MOF = 'M' WE = 'C' NCOMM = INDEX(FITMET,',') CBETA = FITMET(NCOMM+1:20) IF (NCOMM.GT.0) THEN CALL GENCNV(CBETA,2,1,IDUM,BETA,DDUM,ISTAT) C CALL USRINP(BETA,1,'R',CBETA) IF (ABS(BETA).LT.1.0E-30) THEN BETA = 4.0 ENDIF ENDIF ELSE IF (FITMET(1:2).EQ.'MI') THEN MOF = 'M' WE = 'N' NCOMM = INDEX(FITMET,',') CBETA = FITMET(NCOMM+1:20) IF (NCOMM.GT.0) THEN CALL GENCNV(CBETA,2,1,IDUM,BETA,DDUM,ISTAT) C CALL USRINP(BETA,1,'R',CBETA) IF (ABS(BETA).LT.1.0E-30) THEN BETA = 4.0 ENDIF ENDIF ELSE IF (FITMET(1:2).EQ.'GU') THEN MOF = 'G' WE = 'C' BETA = 0.0 ELSE IF (FITMET(1:2).EQ.'GI') THEN MOF = 'G' WE = 'N' BETA = 0.0 ELSE MOF = 'M' WE = 'C' BETA = 4.0 ENDIF C C *** get fit option CALL STKRDC('INPUTC',1,21,4,IAC,FITOPT,KUN,KNUL,ISTAT) !fit option CALL UPCAS(FITOPT,FITOPT) SIGFI = FITOPT(1:1) SKYFI = FITOPT(2:2) POSFI = FITOPT(3:3) TILT = FITOPT(4:4) C C *** get the fix option, sigma, position, sky background, sky tilt 1000 CONTINUE IF (SKYFI.EQ.'Y' .AND. POSFI.EQ.'Y') THEN CALL STTPUT('*** WARNING: Illegal combination, '// 2 '... try again',ISTAT) FITOPT = 'YNNN' CALL STKWRC('INPUTC',1,FITOPT,21,4,KUN,ISTAT) CALL STKPRC(STRING,'INPUTC',1,21,4,IAC,FITOPT,KUN,KNUL,ISTAT) CALL UPCAS(FITOPT,FITOPT) SIGFI = FITOPT(1:1) SKYFI = FITOPT(2:2) POSFI = FITOPT(3:3) TILT = FITOPT(4:4) GO TO 1000 END IF C IF (SIGFI.NE.'N') SIGFI = 'Y' IF (POSFI.NE.'Y') POSFI = 'N' IF (SKYFI.NE.'Y') THEN SKYFI = 'N' IF (POSFI.NE.'Y') THEN TILT = FITOPT(4:4) END IF IF (TILT.NE.'Y') THEN TILT = 'N' ENDIF END IF C CALL STKRDC('INPUTC',1,41,1,IAC,AVEOPT,KUN,KNUL,ISTAT) ! aver. option CALL UPCAS(AVEOPT,AVEOPT) INME = AVEOPT C CALL STKRDC('INPUTC',1,61,1,IAC,DIS,KUN,KNUL,ISTAT) ! display CALL UPCAS(DIS,DIS) IF (DIS.EQ.'Y') THEN NDISPL = 0 ELSE NDISPL = 1 ENDIF C C *** write the info to the screen WRITE(STRING,9001) FRAME CALL STTPUT(STRING,ISTAT) WRITE(STRING,9002) INTFIL CALL STTPUT(STRING,ISTAT) WRITE(STRING,9003) SIGMA, SAT CALL STTPUT(STRING,ISTAT) WRITE(STRING,9004) PSA,TNIT CALL STTPUT(STRING,ISTAT) C IF (FITMET(1:2).EQ.'MU') THEN WRITE(STRING,9005) BETA ELSE IF (FITMET(1:2).EQ.'MI') THEN WRITE(STRING,9006) BETA ELSE IF (FITMET(1:2).EQ.'GU') THEN STRING = 'Fit method: GAUSS, UNIFORM weighting' ELSE IF (FITMET(1:2).EQ.'GI') THEN STRING = 'Fit method: GAUSS, 1/N weighting' ENDIF CALL STTPUT(STRING,ISTAT) C WRITE(STRING,9007) FITOPT CALL STTPUT(STRING,ISTAT) C IF (AVEOPT.EQ.'S') THEN STRING = 'Trial values: mean sigma' ELSE IF (AVEOPT.EQ.'B') THEN STRING = 'Trial values: mean sky background' ELSE IF (AVEOPT.EQ.'A') THEN STRING = 'Trial values: mean sigma and mean sky background' ELSE STRING = 'Trial values: none' ENDIF CALL STTPUT(STRING,ISTAT) C C *** set the display flag CALL STKRDI('LOG',4,1,IAC,ODISPL,KUN,KNUL,ISTAT) CALL STKWRI('LOG',NDISPL,4,1,KUN,ISTAT) WRITE(STRING,9021) CALL STTPUT(STRING,ISTAT) C *** do the work ***************************************************** C NNS = 1 NNF = 1 KAL = 1 ALT1 = 1. KSF = 0. IOBJ = 0 IROW = 1 C 1001 CONTINUE CALL INTWRD(TIDINT,IROW,NCP,NHL) ! read the window CALL TBSGET(TIDINT,IROW,SFLAG,ISTAT) IF (SFLAG) THEN D1 = PARINT(1) D2 = PARINT(2) V = PARINT(3) B = PARINT(4) U = PARINT(5) D3 = PARINT(6) D4 = PARINT(7) P(1) = PARINT(8) P(2) = PARINT(9) P(3) = PARINT(10) BET = PARINT(11) SCAP = PARINT(12) D6 = PARINT(13) FL = INT(PARINT(14)) C DO 1011 IS = 1,NCP P((IS-1)*4+4) = FITCMP((IS-1)*6+1) P((IS-1)*4+5) = FITCMP((IS-1)*6+2) P((IS-1)*4+6) = FITCMP((IS-1)*6+3) P((IS-1)*4+7) = FITCMP((IS-1)*6+4) USC(IS) = FITCMP((IS-1)*6+5) USE(IS) = FITCMP((IS-1)*6+6) 1011 CONTINUE C DO 1012 IH = 1,NHL PB((IH-1)*3+1) = FITHOL((IH-1)*3+1) PB((IH-1)*3+2) = FITHOL((IH-1)*3+2) PB((IH-1)*3+3) = FITHOL((IH-1)*3+3) 1012 CONTINUE C DO 1013 K3 = 1,NTM ! initialize the flag array IF (FLGCMP(K3).GE.2) THEN FLGCMP(K3) = 1 ENDIF NONF(K3) = FLGCMP(K3) KFR(K3) = 0 1013 CONTINUE C FL = 0 IF (TILT.EQ.'N') THEN P(1) = 0. P(2) = 0. ENDIF IF (FON1.GT.0 .OR. (INME.EQ.'B' .OR. INME.EQ.'A')) THEN IF (FON.NE.0) THEN P(3) = FON ELSE FON = P(3) END IF P(1) = 0. P(2) = 0. END IF C D1 = D1+ISX D2 = D2+ISY IPX = D1 IPY = D2 NP = D3 NC = D4 C ASSIGN 30011 TO PREVET ! PREPARA VETTORE GO TO 30001 30011 CONTINUE C IF (NCOM.GT.0 .AND. FL.EQ.0) THEN IF (NHL.GT.0) THEN ASSIGN 30012 TO PREMAS ! PREPARA MASCHERA GO TO 30002 30012 CONTINUE END IF C KP = 0 LJ1 = IPX-START(1)+1 LJ3 = IPY-START(2)+1 DO 1014 J5 = LJ3,LJ3+NC-1 JOLD = J5-LJ3+1 CALL REALIN(NPL,NL,J5,LJ1,NP,MADRID(IPNTR),RIA) DO 1015 I = 1,NP IF (MASK(JOLD,I).EQ.0) THEN VAL = RIA(I) - FOG IF (VAL.LT.SAT.AND.VAL.GT.0.) THEN KP = KP+1 IVX(KP) = I IVY(KP) = JOLD VZ(KP) = VAL IF (WE.EQ.'N') THEN IF (VAL.GE.1) THEN WEI(KP) = 1./VAL ENDIF ELSE WEI(KP) = 1 ENDIF END IF END IF 1015 CONTINUE 1014 CONTINUE IF (KP.GT.NCOM*4+3) THEN ASSIGN 30013 TO FRASUZ !FRASER SUZUKI GO TO 30003 30013 CONTINUE IF (WFLAG.EQ.1) THEN GR = 1 ELSE GR = 2 END IF WRITE (STRING,9012) IDNGRP,NCP,NITER CALL STTPUT(STRING,ISTAT) C ASSIGN 30014 TO RIPVET !ripristina vettore GO TO 30004 30014 CONTINUE C C IF (GR.EQ.1) THEN DO 1016 K3 = 1,NTM FLGCMP(K3) = NONF(K3) 1016 CONTINUE C END IF C IF (AIN.GT.0) THEN ASSIGN 30015 TO CALINT !calcola interquartile GO TO 30005 30015 CONTINUE END IF DO 1017 KKL = 1,NCOM IKKL = KTP+(KKL-1)*KTS WRITE(STRING,9013) 2 (SIGP(KL),KL=IKKL+1,IKKL+KTS) IF (KKL.EQ.1) THEN IF (KTP.GT.0) THEN WRITE(STRING(53:),9014) 2 (SIGP(KL),KL=1,KTP) ENDIF ENDIF CALL STTPUT(STRING,ISTAT) 1017 CONTINUE CALL STTPUT(' ',ISTAT) ELSE GR = 2 WRITE(STRING,9011) CALL STTPUT(STRING,ISTAT) END IF C C *** write the output row for this window PARINT(1) = D1 PARINT(2) = D2 PARINT(3) = V PARINT(4) = B PARINT(5) = U PARINT(6) = FLOAT(NP) PARINT(7) = FLOAT(NC) PARINT(8) = P(1) PARINT(9) = P(2) PARINT(10) = P(3) PARINT(11) = BETA PARINT(12) = SCAP PARINT(13) = FLOAT(NITER) PARINT(14) = FLOAT(GR) C DO 1018 IS = 1,NCP FITCMP((IS-1)*6+1) = P((IS-1)*4+4) FITCMP((IS-1)*6+2) = P((IS-1)*4+5) FITCMP((IS-1)*6+3) = P((IS-1)*4+6) FITCMP((IS-1)*6+4) = P((IS-1)*4+7) FITCMP((IS-1)*6+5) = USC(IS) FITCMP((IS-1)*6+6) = USE(IS) 1018 CONTINUE C DO 1019 IH = 1,NHL FITHOL((IH-1)*3+1) = PB((IH-1)*3+1) FITHOL((IH-1)*3+3) = PB((IH-1)*3+2) FITHOL((IH-1)*3+3) = PB((IH-1)*3+3) 1019 CONTINUE C CALL INTWWR(TIDINT,IROW,NCP,NHL) DO 1020 KRO = 1,NCP RK4 = P(4*KRO) IF (FLGCMP(KRO).EQ.0) THEN JL0 = JL0+1 IF (RK4.GT.1) THEN VEZE(JL0) = -2.5*ALOG10(RK4) ELSE VEZE(JL0) = 0 END IF RMAY = AMAX1(RMAY,VEZE(JL0)) RMIY = AMIN1(RMIY,VEZE(JL0)) C ELSE IF (FLGCMP(KRO).EQ.1) THEN JL1 = JL1+1 IF (RK4.GT.1) THEN VEUN(JL1) = -2.5*ALOG10(RK4) ELSE VEUN(JL1) = 0 END IF RMAY = AMAX1(RMAY,VEUN(JL1)) RMIY = AMIN1(RMIY,VEUN(JL1)) C ELSE IF (FLGCMP(KRO).GT.1)THEN JL2 = JL2+1 IF (RK4.GT.1) THEN VEMA(JL2) = -2.5*ALOG10(RK4) ELSE VEMA(JL2) = 0 END IF RMAY = AMAX1(RMAY,VEMA(JL2)) RMIY = AMIN1(RMIY,VEMA(JL2)) END IF 1020 CONTINUE C CALL INTDWR(TIDINT,NGRP,NOBJ,NINT,SAT,FAT,SIGMA, 2 BETA,SOFOT,AIN,FOG) IF (NHL.GT.0) THEN ASSIGN 30016 TO AZZMAS !AZZERA MASCHERA GO TO 30006 30016 CONTINUE END IF END IF END IF IROW = IROW + NCP + NHL IF (IROW.LE.NRINT) THEN GO TO 1001 ENDIF C C *** display the histogram PAS = .5 MAXX = 0 DO 2001 I = 1,JL0 NCAL = (RMAY-VEZE(I))/PAS+1 IF (NCAL.GT.100) NCAL=100 ISTZE(NCAL) = ISTZE(NCAL)+1 MAXX = MAX0(MAXX,ISTZE(NCAL)) 2001 CONTINUE DO 2002 I = 1,JL1 NCAL = (RMAY-VEUN(I))/PAS+1 IF (NCAL.GT.100) NCAL=100 ISTUN(NCAL) = ISTUN(NCAL)+1 MAXX = MAX0(MAXX,ISTUN(NCAL)) 2002 CONTINUE DO 2003 I = 1,JL2 NCAL = (RMAY-VEMA(I))/PAS+1 IF (NCAL.GT.100) NCAL=100 ISTMA(NCAL) = ISTMA(NCAL)+1 MAXX = MAX0(MAXX,ISTMA(NCAL)) 2003 CONTINUE NCM = (RMAY-RMIY)/PAS+1 IF (JL0.GT.0) THEN WRITE(STRING,9021) CALL STTPUT(STRING,ISTAT) WRITE(STRING,9022) JL0 CALL STTPUT(STRING,ISTAT) CALL STTPUT(' ',ISTAT) WRITE (STRING,9023) CALL STTPUT(STRING,ISTAT) CALL STTPUT(' ',ISTAT) DO 2004 I = 1,NCM NUX = 50.*ISTZE(I)/MAXX ANC = RMAY-PAS*(I-1) HH = 10.**(-.4*ANC) IF (NUX.GT.0) THEN XXX = ' ' DO 2005 K13 = 1,NUX XXX(K13:K13) = 'X' 2005 CONTINUE WRITE(STRING,9024) ISTZE(I),ANC,HH,XXX CALL STTPUT(STRING,ISTAT) ELSE WRITE(STRING,9025) ISTZE(I),ANC,HH CALL STTPUT(STRING,ISTAT) END IF 2004 CONTINUE END IF IF (JL1.GT.0) THEN WRITE(STRING,9021) CALL STTPUT(STRING,ISTAT) WRITE(STRING,9026) JL1 CALL STTPUT(STRING,ISTAT) CALL STTPUT(' ',ISTAT) DO 2006 I = 1,NCM NUX = 50.*ISTUN(I)/MAXX ANC = RMAY-PAS*(I-1) HH = 10.**(-.4*ANC) IF (NUX.GT.0) THEN XXX = ' ' DO 2007 K13 = 1,NUX XXX(K13:K13) = 'X' 2007 CONTINUE WRITE(STRING,9024) ISTUN(I),ANC,HH,XXX CALL STTPUT(STRING,ISTAT) ELSE WRITE(STRING,9025) ISTUN(I),ANC,HH CALL STTPUT(STRING,ISTAT) END IF 2006 CONTINUE END IF C IF (JL2.GT.0) THEN WRITE(STRING,9021) CALL STTPUT(STRING,ISTAT) WRITE(STRING,9027) JL2 CALL STTPUT(STRING,ISTAT) CALL STTPUT(' ',ISTAT) DO 2008 I = 1,NCM NUX = 50.*ISTMA(I)/MAXX ANC = RMAY-PAS*(I-1) HH = 10.**(-.4*ANC) IF (NUX.GT.0) THEN XXX = ' ' DO 2009 K13 = 1,NUX XXX(K13:K13) = 'X' 2009 CONTINUE WRITE(STRING,9024) ISTMA(I),ANC,HH,XXX CALL STTPUT(STRING,ISTAT) ELSE WRITE(STRING,9025) ISTMA(I),ANC,HH CALL STTPUT(STRING,ISTAT) END IF 2008 CONTINUE END IF CALL TBTCLO(TIDINT,ISTAT) CALL STKWRI('LOG',ODISPL,4,1,KUN,ISTAT) CALL STSEPI C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Procedure PREVET C C------------------------------------------------------------------------------- 30001 CONTINUE DO 30101 LI = 1,3 PA(LI) = P(LI) 30101 CONTINUE NCOM = 0 DO 30201 LI = 1,NCP IF (FLGCMP(LI).GE.1) THEN IDE = (LI-1)*4+3 NCOM = NCOM+1 IDI = (NCOM-1)*4+3 C IF (INME.EQ.'S' .OR. INME.EQ.'A') THEN P(IDE+4) = SIGG ELSE P(IDE+4) = SIGMA ENDIF C DO 30301 LLI = 1,4 PA(IDI+LLI) = P(IDE+LLI) 30301 CONTINUE PA(IDI+1) = PA(IDI+1)*ALT1 FLGCMP(LI) = 1 NONF(LI) = 1 END IF 30201 CONTINUE GO TO PREVET C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Procedure PREMAS C C------------------------------------------------------------------------------- 30002 CONTINUE !PREPARA MASCHERA DO 30102 LI = 1,NHL ICB = (LI-1)*3 RR = PB(ICB+1)**2 L4 = AMAX1(1.,PB(ICB+3)-PB(ICB+1)) L5 = AMIN1(FLOAT(NC),PB(ICB+3)+PB(ICB+1)+.99) L6 = AMAX1(1.,PB(ICB+2)-PB(ICB+1)) L7 = AMIN1(FLOAT(NP),PB(ICB+2)+PB(ICB+1)+.99) DO 30202 L2 = L6,L7 DO 30302 L3 = L4,L5 DELT = (FLOAT(L2)-PB(ICB+2))**2 + 2 (FLOAT(L3)-PB(ICB+3))**2 IF (RR.GE.DELT) THEN MASK(L3,L2) = 1 END IF 30302 CONTINUE 30202 CONTINUE 30102 CONTINUE GO TO PREMAS C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Procedure FRASUZ C C------------------------------------------------------------------------------- 30003 CONTINUE ! TO FRASER SUZUKI WFLAG = 1 LFLAG = 0 SLU = 0.00001 NITER = 0 SCAP = 10.**15 US = 10.**10 KPR1 = 0 NRFR = 0 IF (POSFI.EQ.'Y' .AND. SIGFI.EQ.'Y') THEN CALL ELMRPF(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA, 2 SQM,LFLAG,WEI,SIGP) KTP = 1 KTS = 1 DO 30103 K9 = 1,NCOM NONF(K9) = FLGCMP(K9) 30103 CONTINUE IF (LFLAG.EQ.0) THEN LFLAG = 1 ENDIF NITER = 1 ELSE 39803 CONTINUE IF (LFLAG.NE.0) GO TO 39903 39703 CONTINUE RTELOG = ABS(US).LE.SLU .OR. NITER.GE.TNIT .OR. 2 LFLAG.NE.0 IF (RTELOG) GO TO 39603 NITER = NITER+1 IF (MOF.EQ.'M') THEN IF (NITER.LE.50.) THEN D = 0.7 ELSE IF (NITER.GT.50. .AND. NITER.LE.80.) THEN D = 0.3 ELSE D = 0.7 END IF END IF PES(4) = 2 IF (NITER.LE.1.) PES(4)=.1 LG = 0 C IF (SIGFI.EQ.'Y' .AND. SKYFI.EQ.'Y' .AND. 2 POSFI.NE.'Y') THEN CALL ELMRF(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA, 2 SQM,LFLAG,WEI,SIGP) KTP = 0 KTS = 3 ELSE IF (SIGFI.EQ.'N' .AND. SKYFI.EQ.'Y' .AND. 2 POSFI.NE.'Y') THEN CALL ELMRFV(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA, 2 SQM,LFLAG,WEI,SIGP) KTP = 0 KTS = 4 ELSE IF (SIGFI.EQ.'Y' .AND. TILT.EQ.'N' .AND. 2 POSFI.NE.'Y') THEN CALL ELMRR(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA, 2 SQM,LFLAG,WEI,SIGP) KTP = 1 KTS = 3 ELSE IF (SIGFI.EQ.'Y' .AND. TILT.EQ.'Y' .AND. 2 POSFI.NE.'Y') THEN CALL ELMRX(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA, 2 SQM,LFLAG,WEI,SIGP) KTP = 3 KTS = 3 ELSE IF (SIGFI.EQ.'N' .AND. TILT.EQ.'N' .AND. 2 POSFI.NE.'Y') THEN CALL ELMRRV(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA, 2 SQM,LFLAG,WEI,SIGP) KTP = 1 KTS = 4 ELSE IF (SIGFI.EQ.'N' .AND. TILT.EQ.'Y' .AND. 2 POSFI.NE.'Y') THEN CALL ELMRV(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA, 2 SQM,LFLAG,WEI,SIGP) KTP = 3 KTS = 4 ELSE IF(SIGFI.EQ.'N' .AND. POSFI.EQ.'Y') THEN CALL ELMRPV(IVX,IVY,VZ,KP,PA,D,PES,NCOM,BETA, 2 SQM,LFLAG,WEI,SIGP) KTP = 1 KTS = 2 END IF C IF (LFLAG.EQ.0) THEN US = (SCAP-SQM)/SQM SCAP = SQM END IF C ASSIGN 30017 TO CONPAR GO TO 30007 30017 CONTINUE GO TO 39703 39603 CONTINUE C IF (NRFR.GT.0 .AND. LFLAG.EQ.0) THEN NPSG = 0 HGU = 0 DO 30203 K3 = 1,NCP IF (KFR(K3).GT.0) THEN IND = K3*4 IF (P(IND).GT.HGU) THEN HGU = P(IND) NPSG = K3 END IF END IF 30203 CONTINUE NRFR=NRFR-1 C *** !calcola code IND=(NPSG-1)*4+5 DO 30303 K7=1,NCOM K6=K7*4 ABC = ((PA(K6+1)-P(IND))**2 + 2 (PA(K6+2)-P(IND+1))**2)/PA(K6+3)**2 IF (ABS(BETA).LT.1.0E-30) THEN ABC = PA(K6)*EXP(-ABC*4*ALOG(2.)) ELSE ABC = PA(K6)*(1.+ABC)**(-BETA) END IF HGU = HGU-ABC 30303 CONTINUE HGU = HGU+P(3)-PA(3) WRITE (STRING,9031) NITER,NPSG CALL STTPUT(STRING,ISTAT) C IF (HGU.GT.SOFOT) THEN FLGCMP(NPSG) = 1 SLU = 0.00001 NITER = 0 SCAP = 10.**15 US = 10.**10 IND = 1 IND2 = 0 C IF (NPSG.LE.IND) GO TO 30503 IF (FLGCMP(IND).GE.1) IND2=IND2+1 IND = IND+1 30503 CONTINUE C IF (IND2.LT.NCOM) THEN DO 30603 LI = NCOM,IND2+1,-1 IND = LI*4 DO 30703 K3 = IND,IND+3 PA(K3+4) = PA(K3) 30703 CONTINUE 30603 CONTINUE END IF NCOM = NCOM+1 IND = (IND2+1)*4 IND1 = NPSG*4 DO 30803 K3 = IND,IND+3 PA(K3) = P(IND1+K3-IND) 30803 CONTINUE ELSE WRITE (STRING,9032) CALL STTPUT(STRING,ISTAT) END IF KFR(NPSG)=-KFR(NPSG) ELSE IF (LFLAG.EQ.0) THEN LFLAG = 1 ENDIF ENDIF GO TO 39803 39903 CONTINUE END IF GO TO FRASUZ C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Procedure RIPVET C C------------------------------------------------------------------------------- 30004 CONTINUE !RIPRISTINA VETTORE DO 30104 LI=1,3 P(LI)=PA(LI) 30104 CONTINUE IF (GR.NE.2) THEN IF (INME.EQ.'B'.OR.INME.EQ.'A') THEN RKN = 1./NNF NNF = NNF+1 FON = (FON+PA(3)*RKN)/(1+RKN) END IF END IF NCO=0 DO 30204 LI = 1,NCP IF (NONF(LI).GE.1) THEN IDE = (LI-1)*4+3 NCO = NCO+1 IDI = (NCO-1)*4+3 IF (NONF(LI).EQ.1) THEN KAS = 1./KAL KAL = KAL+1 ALT1 = (ALT1+PA(IDI+1)/P(IDE+1)*KAS)/(1.+KAS) DO 30304 LLI = 1,4 P(IDE+LLI) = PA(IDI+LLI) 30304 CONTINUE IF (SIGFI.EQ.'N' .AND. 2 (INME.EQ.'S'.OR.INME.EQ.'A')) THEN RKN = 1./NNS NNS = NNS+1 SIGG = (SIGG+PA(IDI+4)*RKN)/(1+RKN) END IF END IF END IF 30204 CONTINUE GO TO RIPVET C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Procedure CALINT C C------------------------------------------------------------------------------- 30005 CONTINUE ! CALCOLA INTERQUARTILE IF (BETA.GT.0.0) BEBE=-1./BETA C DO 30105 JIL = 1,NCP JIL4 = JIL*4 USC(JIL) = 0. USE(JIL) = 0. IF (FLGCMP(JIL).EQ.1) THEN SIGQ = P(JIL4+3)**2 IF (ABS(BETA).LT.1.0E-30) THEN ALAL = 4*ALOG(2.) SMT = AIN*EXP((1./SIGQ)*ALAL) IF (P(JIL4).GE.SMT) THEN RAG(JIL) = P(JIL4+3)**2* 2 (-ALOG(AIN/P(JIL4))/ALAL) ELSE RAG(JIL) = 1.0 END IF ELSE BEBE = -1./BETA SMT = AIN*((1./SIGQ)+1.)**BETA IF (P(JIL4).GE.SMT) THEN RAG(JIL) = P(JIL4+3)**2* 2 ((AIN/P(JIL4))**BEBE-1.) ELSE RAG(JIL) = 1.0 END IF END IF END IF 30105 CONTINUE C DO 30205 K7=1,NCP K74 = K7*4 IF (FLGCMP(K7).EQ.1) THEN SQMS = 0. NPS = 0 DO 30305 IQ = 1,NCA RM(IQ) = 0 30305 CONTINUE DO 30405 IQ=1,KP DPS = (IVX(IQ)-P(K74+1))**2 + 2 (IVY(IQ)-P(K74+2))**2 IF (DPS.LE.RAG(K7) .AND. VZ(IQ).LT.SAT) THEN NPS = NPS+1 VG = P(3) DO 30505 JQ = 1,NCP JQ4 = JQ*4 IF (FLGCMP(JQ).EQ.1) THEN GIG = P(JQ4+3)**2. EEE = ((P(JQ4+1)-IVX(IQ))**2 + 2 (P(JQ4+2) - IVY(IQ))**2)/GIG IF (ABS(BETA).LT.1.0E-30) THEN VG = VG+P(JQ4)*EXP(-EEE*ALAL) ELSE VG = VG+P(JQ4)*(1.+EEE)**(-BETA) END IF END IF 30505 CONTINUE C SQMS = SQMS+(VZ(IQ)-VG)**2*WEI(IQ) VALME(NPS) = (VZ(IQ)-VG)/SQRT(VG) ENDIF 30405 CONTINUE C DO 30705 IQ = 2,NPS A = VALME(IQ) DO 30605 JS = IQ-1,1,-1 IF (VALME(JS) .LE. A) THEN GO TO 31605 ELSE VALME(JS+1) = VALME(JS) ENDIF 30605 CONTINUE JS = 0 31605 CONTINUE VALME(JS+1) = A 30705 CONTINUE C IRINT = FLOAT(NPS)/4. + 0.5 RINTQ = VALME(IRINT*3)-VALME(IRINT) C IF (NPS.LE.4) THEN NPS = 5 ENDIF SQMS = SQMS/FLOAT(NPS-4) IF (WE.NE.'N') THEN SQMS = SQMS/P(K74) ENDIF IF (RINTQ.LT.0) THEN RINTQ = 50 ENDIF USC(K7) = SQMS USE(K7) = RINTQ END IF 30205 CONTINUE GO TO CALINT C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Procedure AZZMA C C------------------------------------------------------------------------------- 30006 CONTINUE !AZZERA MASCHERA DO 30106 LI=1,NHL ICB = (LI-1)*3 L4 = AMAX1(1.,PB(ICB+3)-PB(ICB+1)) L5 = AMIN1(FLOAT(NC),PB(ICB+3)+PB(ICB+1)+.99) L6 = AMAX1(1.,PB(ICB+2)-PB(ICB+1)) L7 = AMIN1(FLOAT(NP),PB(ICB+2)+PB(ICB+1)+.99) DO 30206 L2 = L6,L7 DO 30306 L3 = L4,L5 MASK(L3,L2) = 0 30306 CONTINUE 30206 CONTINUE 30106 CONTINUE GO TO AZZMAS C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Procedure CONPAR C C------------------------------------------------------------------------------- 30007 CONTINUE !CONTROLLA PARAMETRI IF(LFLAG.EQ.0) THEN KG = 0 IP = 0 IPA = 0 LFL2 = 0 DO 30107 K3 = 1,NCOM IF (LFL2.EQ.0) THEN IPA = IPA+1 NONF(IPA) = FLGCMP(IPA) 30207 CONTINUE IF (FLGCMP(IPA).NE.0 .OR. IPA.GE.NCP) THEN GOTO 30307 ENDIF IPA = IPA+1 NONF(IPA) = FLGCMP(IPA) GOTO 30207 30307 CONTINUE INK1 = (IPA-1)*4 INK = (K3-1)*4 IF (PA(INK+4).LT.0) THEN WFLAG = 2 LFLAG = 2 NONF(IPA) = 3 ! Altezza neagtive END IF IF (PA(INK+7).LT.0. .OR. PA(INK+7).GT.100.) THEN WFLAG = 2 LFLAG = 2 NONF(IPA) = 5 !Sigma errato END IF IF (PA(INK+4).GT.1.E10) THEN WFLAG = 2 LFLAG = 2 NONF(IPA) = 3 !Alt. magg. di 10.000.000 END IF DXY2 = (PA(INK+5) - P(INK1+5))**2 + 2 (PA(INK+6) - P(INK1+6))**2 IF (DXY2.GT.PSAS) THEN WFLAG = 2 LFLAG = 2 NONF(IPA) = 4 !La posizione del massimo END IF !si e' spostata piu' di PSA pixels IF (LFLAG.EQ.2) THEN IF (KFR(IPA).EQ.0) THEN NRFR = NRFR+1 KFR(IPA) = K3 FLGCMP(IPA) = 0 LFLAG = 0 KPR1 = 1 ELSE LFL2 = 1 END IF END IF END IF 30107 CONTINUE C IF (KPR1.EQ.1 .AND. LFL2.EQ.0) THEN IF (NRFR.GT.1) THEN WRITE (STRING,9035) NRFR CALL STTPUT(STRING,ISTAT) ELSE STRING = 'No convergence on 1 star' CALL STTPUT(STRING,ISTAT) ENDIF C KPR1 = 0 SLU = 0.00001 NITER = 0 SCAP = 10.**15 US = 10.**10 C ASSIGN 30018 TO PREVET ! #PREPARA VETTORE GO TO 30001 30018 CONTINUE IF (NCOM.LE.0) THEN LFLAG = 2 ENDIF ENDIF ELSE DO 30407 K3 = 1,NCOM FLGCMP(K3) = 6 NONF(K3) = 6 30407 CONTINUE ENDIF C IF (LFL2.EQ.1) THEN WFLAG = 2 LFLAG = 2 WRITE(STRING,9037) CALL STTPUT(STRING,ISTAT) DO 30408 K3 = 1,NCOM INK = (K3-1)*4 WRITE(STRING,9036) K3,(PA(INK+JIL),JIL=4,7) CALL STTPUT(STRING,ISTAT) 30408 CONTINUE END IF GO TO CONPAR END