C @(#)rfotexamin.for 17.1.1.1 (ES0-DMD) 01/25/02 17:18:16 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 EXAMIN C+++ C.IDENTIFICATION: EXAMIN C.PURPOSE: The examine the quality of the fitted object C.AUTHOR: R. Buonanno, G. Buscema, C. Corsi, I. Ferraro, G. Iannicola C Osservatorio Astronomico di Roma C.VERSION: 870328 Created C.VERSION: 870417 Modified for idi interfaces M.Pucillo, OAT C.VERSION: 871216 Modified for integration into MIDAS Rein Warmels, ESO C.VERSION: 881028 Modifications for portable MIDAS R. Buonanno C.VERSION: 890320 Display part based on MIDAS plot software; C. modifications for portability Rein H. Warmels C.VERSION: 890406 Rewritten for portable MIDAS; C Using MIDAS table interfaces C Based on MIDAS graphics interfaces Rein H. Warmels C--- IMPLICIT NONE INCLUDE 'MID_REL_INCL:RFOTDECL.INC' C INTEGER NHIST INTEGER NTO PARAMETER (NHIST=512) PARAMETER (NTO=10000) C INTEGER MADRID(1) INTEGER IFREQ(NHIST) INTEGER IND(NTO,2) INTEGER ED, EC, EL INTEGER IS, IROW, ID, IDUM INTEGER IHN INTEGER IM, IT, IB INTEGER IX1, IX2 INTEGER I INTEGER IX, IXMAX, IYMAX INTEGER ISTAT,IAC INTEGER ICOUNT, ICOL INTEGER K, KWR, KK INTEGER NCP,NHL, NBINS INTEGER NS, NGRP, NONC INTEGER NNN, NACT, NCHAR INTEGER NPIX, NPOS INTEGER KUN,KNUL INTEGER TIDINT INTEGER NCINT,NRINT,NSINT,NACINT,NARINT INTEGER PLMODE,ACCESS,ILOG INTEGER MAX INTEGER GRE INTEGER TINULL INTEGER NORDAT,PLTHST, PLTTBL DOUBLE PRECISION TDNULL,DDUM REAL TRNULL REAL RFREQ(NHIST) REAL SQ(NTO),SIQ(NTO),H(NTO) REAL SSIQ(NHIST), USC(NHIST), QP(NHIST) REAL SIGMA(2), RMX(2) REAL XFIT(NHIST),YFIT(NHIST) REAL XPOS(NTO), YPOS(NTO) REAL RINPUT(3), MNMX(2) REAL WIND(8), SCALES(3), OFFSET(2) REAL CINT(NHIST) REAL DINT(NHIST) REAL PAR(3) REAL FOPT(3) REAL XMIN, XMAX, AXMIN, AXMAX, X1, X0 REAL STEP REAL XMEAN, XSTD REAL SQR, SL, SQM, FG REAL XCUR, YCUR REAL SQMAX, SMAX, SMIN, SQMIN REAL HMIN, HMAX REAL YOFF, DELTX REAL RSMAX, RSQMAX REAL XT C CHARACTER*60 INTFIL CHARACTER*80 STRING CHARACTER*72 CINPUT CHARACTER*17 COLUMN CHARACTER*60 LABEL1,LABEL2,LABEL3 CHARACTER*4 XFRAME, YFRAME CHARACTER*1 OPTION*1, RET*1, KEY*1 C INCLUDE 'MID_INCLUDE:TABLES.INC' COMMON /VMR/MADRID INCLUDE 'MID_INCLUDE:TABLED.INC' C DATA FOPT/0.0,999,0.0/ DATA PLMODE/1/ DATA ACCESS/0/ DATA ILOG/0/ DATA YOFF/0.0/ DATA XFRAME/'AUTO'/ DATA YFRAME/'AUTO'/ DATA SCALES/0.0,0.0,0.0/ DATA OFFSET/-999., -999./ C 9001 FORMAT('*** WARNING: Something wrong with component ', 2 'in group: ',I8) 9002 FORMAT('*** INFO: Number of groups found: ',I6) 9003 FORMAT(' Number of groups used: ',I6) C C *** start of executable code *************************************** CALL STSPRO('EXAMINE') CALL TBMNUL(TINULL,TRNULL,TDNULL) C CALL STKRDC('P1',1,1,60,IAC,INTFIL,KUN,KNUL,ISTAT) 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: Error opening intermediate table file' 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 CALL STKRDR('INPUTR',1,2,IAC,MNMX,KUN,KNUL,ISTAT) HMIN = MNMX(1) HMAX = MNMX(2) IF (HMAX.EQ.0.0) THEN HMIN = 0.0 HMAX = 10.**15 ENDIF C RET = 'Y' NS = 0 IROW = 1 NGRP = 0 C C *** start loop for reading components 1001 CONTINUE CALL INTWRD(TIDINT,IROW,NCP,NHL) QP(1) = PARINT(8) QP(2) = PARINT(9) QP(3) = PARINT(10) GRE = INT(PARINT(14)) DO 1011 IS = 1,NCP QP((IS-1)*4+4) = FITCMP((IS-1)*6+1) QP((IS-1)*4+5) = FITCMP((IS-1)*6+2) QP((IS-1)*4+6) = FITCMP((IS-1)*6+3) QP((IS-1)*4+7) = FITCMP((IS-1)*6+4) USC(IS) = FITCMP((IS-1)*6+5) SSIQ(IS) = FITCMP((IS-1)*6+6) 1011 CONTINUE C IF (GRE.EQ.1 .OR. GRE.EQ.3) THEN DO 1021 K = 1,NCP IHN = 4*K IF (QP(IHN).GE.HMIN .AND. QP(IHN).LT.HMAX) THEN NONC = FLGCMP(K) IF (NONC.EQ.1 .OR. 2 (NONC.EQ.2 .AND. RET.EQ.'N')) THEN NS = NS+1 IND(NS,1) = IROW IND(NS,2) = K SQ(NS) = USC(K) SIQ(NS) = SSIQ(K) ID = 4*K H(NS) = QP(ID) IF (H(NS).LE.0) THEN WRITE(STRING,9001) IDNGRP CALL STTPUT(STRING,ISTAT) NS = NS - 1 ENDIF IF (NONC.EQ.2 .AND. RET.EQ.'N') THEN KWR = 1 FLGCMP(K) = 1 ENDIF ENDIF ENDIF 1021 CONTINUE ENDIF C IF (KWR.EQ.1) THEN CALL INTWRD(TIDINT,IROW,NCP,NHL) ENDIF NGRP = NGRP + 1 IROW = IROW + NCP + NHL IF (IROW.LE.NRINT) THEN GO TO 1001 ELSE WRITE (STRING,9002) NGRP CALL STTPUT(STRING,ISTAT) WRITE (STRING,9003) NS CALL STTPUT(STRING,ISTAT) ENDIF C ASSIGN 30011 TO NORDAT GO TO 30001 30011 CONTINUE ASSIGN 30012 TO PLTHST GO TO 30002 30012 CONTINUE ASSIGN 30013 TO PLTTBL GO TO 30003 30013 CONTINUE C C *** terminate AGL operations CALL PTCLOS() CALL STSEPI C C+++ 30001 CONTINUE ! Procedure NORDAT C--- SQMAX = -10.**15 SQMIN = -SQMAX SMAX = SQMAX SMIN = SQMIN DO 30101 I = 1,NS C YY = SQRT(H(I)) C SIQ(I) = SIQ(I)/YY SQMAX = AMAX1(SQMAX,SQ(I)) SQMIN = AMIN1(SQMIN,SQ(I)) SMAX = AMAX1(SMAX,SIQ(I)) SMIN = AMIN1(SMIN,SIQ(I)) 30101 CONTINUE C RSMAX = SMAX RSQMAX = SQMAX GO TO NORDAT C C+++ 30002 CONTINUE ! Procedure PLTHST C--- 9004 FORMAT(' Position of maximum and standard deviation: ', 2 F10.4, '; ', F10.4) C DO 30102 KK = 1,2 STEP = 0.0 ICOUNT = 0 IF (KK.EQ.1) THEN COLUMN = ':CHI_SQ' ELSE COLUMN = ':SIQ' ENDIF C C *** read the table columns CALL TBCSER(TIDINT,COLUMN,ICOL,ISTAT) IF (ISTAT.NE.0 .OR. ICOL.EQ.-1) THEN STRING = '*** FATAL: Problems with finding column'// 2 COLUMN CALL STTPUT(STRING,ISTAT) CALL STSEPI ENDIF C C *** get the minimum and maximum CALL TDUMNX(TIDINT,ICOL,NRINT,0,XMIN,XMAX) IF (XMIN.EQ.XMAX) THEN XMAX = XMIN + 1. ENDIF CALL TDSCAL(XMIN,XMAX,1.,X0,X1,IX,XT,NNN) IF (AXMIN.EQ.TRNULL) THEN AXMIN = X0*10.**IX ELSE AXMIN = XMIN ENDIF IF (AXMAX.EQ.TRNULL) THEN AXMAX = X1*10.**IX ELSE AXMAX = XMAX ENDIF IF (STEP.EQ.0.0) THEN STEP = XT ENDIF C C *** determine the frequency distribution 39802 CONTINUE NBINS = (AXMAX-AXMIN)/STEP+1 NBINS = MIN(NHIST,NBINS) CINT(1) = AXMIN DO 30202 I = 1,NBINS CINT(I+1) = CINT(I) + STEP 30202 CONTINUE CALL TDRSTA(TIDINT,ICOL,NRINT,NBINS,1,CINT, 2 IFREQ,NACT,XMEAN,XSTD) DO 101 I = 1,NBINS RFREQ(I) = FLOAT(IFREQ(I)) 101 CONTINUE IF (NACT.EQ.0) THEN CALL STTPUT(' *** FATAL: No selected entries '// 2 ' or empty column', ISTAT) CALL STSEPI ENDIF C C *** statistics determined; start the display CALL PTKWRR('SCALE',2,SCALES) CALL PTKWRR('OFFS',2,OFFSET) C C *** calculate frame WIND(1) = CINT(1) WIND(2) = CINT(NBINS-1) WIND(3) = 0.0 WIND(4) = 0.0 WIND(5) = 0.0 WIND(6) = 0.0 WIND(7) = 0.0 WIND(8) = 0.0 DO 30302 I = 1,NBINS IYMAX = IFREQ(I) IF (IYMAX.GT.WIND(6)) THEN WIND(6) = IYMAX IXMAX = I ENDIF 30302 CONTINUE IF (WIND(5).EQ.WIND(6)) THEN WIND(6) = WIND(5) + 1 ENDIF C C *** get the frame ranges in x and y CALL GETAXS(XFRAME,WIND(1)) CALL GETAXS(YFRAME,WIND(5)) CALL PTKWRR('XWNDL',4,WIND(1)) CALL PTKWRR('YWNDL',4,WIND(5)) C C *** do the plot setup IF (ICOUNT.GT.0) THEN CALL AGVERS ! reset graphics window CALL AGWDEF(WIND(1),WIND(2),WIND(5),WIND(6)) ELSE CALL PTOPEN(' ','none',ACCESS,PLMODE) ! open the device ENDIF C C *** do the work DO 102 IB = 1,NBINS DINT(IB) = CINT(IB) - STEP/2. 102 CONTINUE CALL PTHIST(NBINS,DINT,RFREQ,FOPT) IF (KK.EQ.1) THEN LABEL1 = 'CHI**2' ELSE LABEL1 = 'SIQ' ENDIF LABEL2 = 'Frequency' LABEL3 = 'TITLE=Table: '//INTFIL CALL PTAXES(WIND(1),WIND(5),LABEL1,LABEL2,LABEL3) C C *** fit the frequency distribution with a gaussian NPIX = NBINS-1 MAX = -32000 DO 30402 IM = 1,NPIX IF (IFREQ(IM).GT.MAX) THEN MAX = IFREQ(IM) IX1 = IM ENDIF 30402 CONTINUE IX2 = IX1 30502 CONTINUE IF (IFREQ(IX2).LE.MAX/2. .OR. IX2 .LE.1) THEN GO TO 30602 ELSE IX2 = IX2 - 1 GO TO 30502 ENDIF C 30602 CONTINUE PAR(1) = MAX PAR(2) = IX1 PAR(3) = IX2 SQR = 10.**15 SL = SQR IT = 0 30702 CONTINUE IF (SL.LE.0.00001 .OR. IT.GT.100) GO TO 30802 IT = IT+1 CALL FIST(IFREQ,NPIX,PAR) SQM = 0. DO 30902 IB = 1,NPIX FG = PAR(1)*EXP(-4.*ALOG(2.)*((IB-PAR(2))/PAR(3))**2) SQM = SQM +(IFREQ(IB)-FG)**2 30902 CONTINUE SQM = SQRT(SQM/NPIX) SL = ABS(SQR-SQM)/SQR SQR = SQM GO TO 30702 30802 CONTINUE RMX(KK) = (PAR(2)-1.5)*STEP+AXMIN SIGMA(KK) = PAR(3)*STEP/SQRT(8*ALOG(2.)) DELTX = (AXMAX-AXMIN)/NHIST C DO 31002 IB = 1,NHIST XFIT(IB) = AXMIN + (IB-1)*DELTX YFIT(IB) = PAR(1) * EXP(-4.*ALOG(2.)* 2 ((XFIT(IB)-RMX(KK))/(PAR(3)*STEP))**2) 31002 CONTINUE CALL PTDATA(0,1,1,XFIT,YFIT,0,NHIST) C *** satisfied or not ? ICOUNT = ICOUNT + 1 IF (ICOUNT.GE.1) THEN IF (KK.EQ.1) THEN STRING = '*** INFO: Histogram of chi square values:' CALL STTPUT(STRING,ISTAT) ELSE STRING = '*** INFO: Histogram of semi interq. interval:' CALL STTPUT(STRING,ISTAT) ENDIF WRITE(STRING,9004) RMX(KK),SIGMA(KK) CALL STTPUT(STRING,ISTAT) STRING = ' Enter new min.,max.,step for '// 2 'histogram [fit ok]: ' NCHAR = INDEX(STRING,':')+1 CALL STKPRC(STRING(1:NCHAR),'INPUTC',1,1,60,IAC,CINPUT, 2 KUN,KNUL,ISTAT) IF (IAC.EQ.0 .OR. CINPUT(1:2).EQ.'OK') THEN GO TO 39902 ELSE CALL GENCNV(CINPUT,2,3,IDUM,RINPUT,DDUM,ISTAT) C CALL USRINP(RINPUT,3,'R',CINPUT) AXMIN = RINPUT(1) AXMAX = RINPUT(2) STEP = RINPUT(3) GO TO 39802 ENDIF ENDIF 39902 CONTINUE 30102 CONTINUE GO TO PLTHST C C+++ 30003 CONTINUE ! Procedure PLTTBL C--- ICOUNT = 0 39903 CONTINUE C C C *** do the display part here WIND(1) = SQMIN/SIGMA(1) - 0.5 WIND(2) = SQMAX/SIGMA(1) + 0.5 WIND(3) = 0.0 WIND(4) = 0.0 WIND(5) = SMIN/SIGMA(2) - 0.5 WIND(6) = SMAX/SIGMA(2) + 0.5 WIND(7) = 0.0 WIND(8) = 0.0 C C *** statistics determined; start the display CALL PTKWRR('SCALE',2,SCALES) CALL PTKWRR('OFFS',2,OFFSET) CALL GETAXS(XFRAME,WIND(1)) CALL GETAXS(YFRAME,WIND(5)) CALL PTKWRR('XWNDL',4,WIND(1)) ! store the axes def.'s CALL PTKWRR('YWNDL',4,WIND(5)) C C *** do the plot setup CALL AGVERS CALL AGWDEF(WIND(1),WIND(2),WIND(5),WIND(6)) C C *** plot the data NPOS = 0 DO 30103 IS = 1,NS IF (IND(IS,1).GT.0) THEN NPOS = NPOS + 1 XPOS(NPOS) = SQ(IS)/SIGMA(1) YPOS(NPOS) = SIQ(IS)/SIGMA(2) ENDIF 30103 CONTINUE C C *** set the window CALL PTDATA(5,0,0,XPOS,YPOS,0.0,NPOS) C C *** plot the frame LABEL1 = 'CHI**2' LABEL2 = 'SIQ' LABEL3 = 'TITLE=Table: '//INTFIL CALL PTAXES(WIND(1),WIND(5),LABEL1,LABEL2,LABEL3) C 39803 CONTINUE STRING = '*** Enter flag option (X, Y, A, R or E) [H]: ' NCHAR = INDEX(STRING,':')+1 OPTION = 'H' CALL STKWRC('INPUTC',1,OPTION,1,1,KUN,ISTAT) ! store the scales CALL STKPRC(STRING(1:NCHAR),'INPUTC',1,1,1,IAC,OPTION, 2 KUN,KNUL,ISTAT) CALL UPCAS(OPTION,OPTION) C IF (OPTION.EQ.'F' .OR. OPTION.EQ.'E') THEN DO 30203 IS = 1,NS ! go through the objects IF (IND(IS,1).LT.0) THEN ! is the object selected IROW = IABS(IND(IS,1)) ! row number of object CALL INTWRD(TIDINT,IROW,NCP,NHL) ! read object components FLGCMP(IND(I,2)) = 2 ! set the component flag CALL INTWWR(TIDINT,IROW,NCP,NHL) ! write the object back ENDIF 30203 CONTINUE GO TO PLTTBL ELSE IF (OPTION.EQ.'X') THEN CALL PTGCUR(XCUR,YCUR,KEY,ISTAT) ! get the cursor XCUR = XCUR*SIGMA(1) DO 30303 IS = 1,NS ! run through objects IF (IND(IS,1).GT.0) THEN ! selected ? IF (SQ(IS).GT.XCUR) THEN ! SQ > x_cursor IND(IS,1) = -IND(IS,1) ! flag the object ENDIF ENDIF 30303 CONTINUE SQMAX = XCUR ! redefined maximum ICOUNT = ICOUNT + 1 ! increase graph counter GO TO 39903 ELSE IF (OPTION.EQ.'Y') THEN CALL PTGCUR(XCUR,YCUR,KEY,ISTAT) ! get the cursor YCUR = YCUR*SIGMA(2) DO 30403 IS = 1,NS IF (IND(IS,1).GT.0) THEN IF (SIQ(IS).GT.YCUR) THEN IND(IS,1) = -IND(IS,1) ENDIF ENDIF 30403 CONTINUE SMAX = YCUR ICOUNT = ICOUNT + 1 GO TO 39903 ELSE IF (OPTION.EQ.'A') THEN CALL PTGCUR(XCUR,YCUR,KEY,ISTAT) ! get the cursor XCUR = XCUR*SIGMA(1) YCUR = YCUR*SIGMA(2) DO 30503 IS = 1,NS ! run through objects IF (IND(IS,1).GT.0) THEN ! selected ? IF (SQ(IS).GT.XCUR) THEN ! SQ > x_cursor IND(IS,1) = -IND(IS,1) ! flag the object ENDIF IF (SIQ(IS).GT.YCUR) THEN IND(IS,1) = -IND(IS,1) ENDIF ENDIF 30503 CONTINUE SQMAX = XCUR ! redefined maximum SMAX = YCUR ICOUNT = ICOUNT + 1 ! increase graph counter GO TO 39903 ELSE IF (OPTION.EQ.'R') THEN DO 30603 IS = 1,NS IF (IND(IS,1).LT.0) THEN IND(IS,1) = -IND(IS,1) ENDIF 30603 CONTINUE SMAX = RSMAX SQMAX = RSQMAX ICOUNT = ICOUNT + 1 GO TO 39903 ELSE IF (OPTION.EQ.'H' .OR. IAC.EQ.0) THEN CALL STTPUT(' A: Flag objects at the right OR '// 2 'above cursor', ISTAT) CALL STTPUT(' X: Flag objects at the right of '// 2 'vertical cursor', ISTAT) CALL STTPUT(' Y: Flag objects above of '// 2 'horizontal cursor', ISTAT) CALL STTPUT(' R: Restore all flags', ISTAT) CALL STTPUT(' H: Display this help information', 2 ISTAT) CALL STTPUT(' F or E: EXIT', ISTAT) GO TO 39803 ELSE CALL STTPUT('*** WARNING: Unknown option', ISTAT) GO TO 39803 ENDIF C END