C @(#)rfotexamin.for	19.1 (ES0-DMD) 02/25/03 13:30:14
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
