C @(#)rfotgroup.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 GROUP
C+++
C.IDENTIFICATION: RFOTGROUP
C.PURPOSE: Group the objects found by the SEARCH command
C.AUTHOR:  R. Buonanno, G. Buscema, C. Corsi, I. Ferraro, G. Iannicola
C          Osservatorio Astronomico di Roma
C.VERSION: 16.09.87 First running version at ESO (outside MIDAS)   R. Buonanno
C          17.12.87 Installation in MIDAS                          R.H. Warmels
C          29.09.88 New version                                    R. Buonanno
C          11.01.89 Include MIDAS for background values            R.H. Warmels
C          06.06.89 MIDAS tables included; partially rewritten     R.H. Warmels
C.VERSION  910122 RHW IMPLICIT NONE added; all variables define
C----
      IMPLICIT  NONE
      INCLUDE   'MID_REL_INCL:RFOTDECL.INC'
C
      INTEGER   IGV
      REAL      VPA
      INTEGER   NDX
      INTEGER   NDY
      INTEGER   MBB
      INTEGER   NCD
      INTEGER   NOB
      INTEGER   NTM
      INTEGER   NTU
      INTEGER   NTAB
      PARAMETER (IGV=100)
      PARAMETER (VPA=.5)
      PARAMETER (NDX=55)
      PARAMETER (NDY=55)
      PARAMETER (MBB=256)
      PARAMETER (NCD=MBB-18)
      PARAMETER (NOB=NCD/3)
      PARAMETER (NTM=NCD/6.5)
      PARAMETER (NTU=100)
      PARAMETER (NTAB=30000)
C
      INTEGER      AREA(4)
      INTEGER      CERSTE,CERINT
      INTEGER      CUU
      INTEGER      EC,EL,ED
      INTEGER      ISTF(NTM),ISTB(IGV),ISTC(IGV)
      INTEGER      IAV, KUN, KNUL, ISTAT
      INTEGER      IXW, IYW
      INTEGER      IMF
      INTEGER      IAC, IPERCE
      INTEGER      IROW, IDUM3, IOBJ, IFL
      INTEGER      I, I3, ISAV, ICA, ITF, ITT, IGRP, INF, IL
      INTEGER      INFB, INFC, INFT, IFNM, INY, INX
      INTEGER      IFX, IFY
      INTEGER      IX, IY
      INTEGER      IXMA, IYMA
      INTEGER      IXMI, IYMI
      INTEGER      J, JIL, JH, JI, JEL
      INTEGER      LX0, LY0
      INTEGER      LX1, LY1
      INTEGER      LFI, LI, L, LMW, LX, LY
      INTEGER      KAA, KN, KC, KS
      INTEGER      K, KSR
      INTEGER      MADRID(1)
      INTEGER      MQQ2, MQQ3, MQQ4, MQQ5, MQQ6, MQQ7, MQQ8, MQQ9
      INTEGER      MAF, MAC, MAT
      INTEGER      NOS(NTAB)
      INTEGER      NST, NSF, NSPR, NSMQ
      INTEGER      NGRP, NOBJ
      INTEGER      NT, NDUM
      INTEGER      NIX, NIY
      INTEGER      NAXIS, NPL, NL	
      INTEGER      NTC, NO, NEX
      INTEGER      NRC, NHL, NCP, NC, NCAL, NBW, NI
      INTEGER      TIDCAT,NCCAT,NRCAT,NSCAT,NACCAT,NARCAT
      INTEGER      TIDINT,NCINT,NRINT,NSINT,NACINT,NARINT
C
C
      REAL         AIN, ANCL    
      REAL         AM1, AM2
      REAL         BETA
      REAL         DAT(7,NTAB)
      REAL         DEX, DEY, DW2
      REAL         DICI, DICS, DIS
      REAL         DWI
      REAL         FOME, FOG, FAR
      REAL         HHS, HV
      REAL         IS(NOB,4)
      REAL         PFA
      REAL         RME, RMA
      REAL         RA, RAS
      REAL         RINPUT(2)
      REAL         RDUM1, RDUM2, RDUM3, RDUM4, RDUM5, RDUM6, RDUM7
      REAL         RDUM9, RDUM10, RDUM11, RDUM12, RDUM13, RDUM
      REAL         SSLA, SUV, SATU, SIGMA, SIV
      REAL         SF(NTM,7)
      REAL         SMRA
      REAL         V(7)
      REAL         VVV(15)
      REAL         XCEN, YCEN
      REAL         ZZ, ZMRA
C
      CHARACTER*60 FRAME,INTFIL,CATFIL
      CHARACTER    TEXT*80,COO*1
      CHARACTER    CHG*40,CHB*40,XXX*60
      CHARACTER    STRING*80
C
      LOGICAL      LOG,LOG1
C
      INCLUDE      'MID_INCLUDE:ST_DEF.INC'
      COMMON       /VMR/MADRID
      INCLUDE      'MID_INCLUDE:ST_DAT.INC'
C
 9001 FORMAT('*** INFO: Reduction factor range',I3,',',I3,' %')
 9002 FORMAT('          Area in pixel units: [',
     2                  I5, ',', I5, ':', I5, ',', I5, ']')
 9003 FORMAT('          Sigma=',F6.2,          '; Beta=',F6.2)  
 9004 FORMAT('          Saturation=',F10.2,'; Threshold=',F10.2)
 9005 FORMAT('          Maximum # objects=',I5,'; Maximum # holes=',I5)
 9006 FORMAT('          Group ',I6,' written:', I6,' star(s); ',
     2                                          I6,' hole(s)')
 9007 FORMAT(80('-'))
 9008 FORMAT(' Frequency of groups with N components')
 9009 FORMAT(' Group   N')
 9010 FORMAT(1X,I5,2X,I2,1X,'I',A40)
 9011 FORMAT(1X,I5,2X,I2,1X,'I')
 9012 FORMAT(' Number of components failed to group: ',I6)
 9013 FORMAT(' # comp.   perc.   mag.    index   H')
C9014 FORMAT(2X,I5,2X,I6,2X,F5.1,2X,F9.0,2X,'I',A40)
 9015 FORMAT(2X,I5,2X,I6,2X,F5.1,2X,F9.0,2X,'I')
 9017 FORMAT(2X,I5,2X,I6,2X,F5.1,2X,F9.0,2X,'I',40('B'))
 9018 FORMAT(2X,I5,2X,I6,2X,F5.1,2X,F9.0,2X,'I',40('G'))
 9019 FORMAT(' Number of components examined: ', I6, ' (total= ', I6,
     2       ' +', I5, ' at the edge)') 
 9020 FORMAT(' Number of components grouped:  ', I6)
 9021 FORMAT(' Number of groups:              ', I6)
C
C *****
      CALL STSPRO('GROUP')
C
C *** find file name and read the data
      CALL STKRDC('IN_A',1,1,60,IAC,FRAME,KUN,KNUL,ISTAT)
      CALL STFOPN(FRAME,D_R4_FORMAT,0,F_IMA_TYPE,IMF,ISTAT)
      CALL STDRDI(IMF,'NAXIS',1,1,IAC,NAXIS,KUN,KNUL,ISTAT)
CC      CALL STFCLO(IMF,ISTAT) 
C 
C *** get the area over which the grouping should be done
      CALL STKRDC('INPUTC',1,1,80,IAC,STRING,KUN,KNUL,ISTAT)
      CALL EXTCOO(IMF,STRING,NAXIS,NDUM,AREA(1),AREA(3),ISTAT)
C
      LX0 = AREA(1)
      LY0 = AREA(2)
      LX1 = AREA(3) - AREA(1) + 1
      LY1 = AREA(4) - AREA(2) + 1
      COO = 'A'
      IF (LX1.GT.0. AND. LY1.GT.0) THEN
         COO = 'S'
      ENDIF
C
C *** get the catalogue table
      CALL STKRDC('IN_B',1,1,60,IAV,CATFIL,KUN,KNUL,ISTAT)    ! get catalogue
C
C *** open the catalogue table
      CALL STECNT('GET',EC,ED,EL)
      CALL STECNT('PUT',1,0,0)
      CALL TBTOPN(CATFIL,F_IO_MODE,TIDCAT,ISTAT)
      IF (ISTAT.NE.0) THEN
         STRING = '*** FATAL: Problems with opening the catalogue'//
     2            ' table; Try again ...'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      ENDIF
C
      CALL TBIGET(TIDCAT,NCCAT,NRCAT,NSCAT,NACCAT,NARCAT,ISTAT)
      IF (ISTAT.NE.0) THEN
         STRING = '*** FATAL: Problems with getting info for'//
     2            ' catalogue table; Try again ...'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      ENDIF
      IF (NRCAT.EQ.0) THEN
         STRING = '*** FATAL: No data in the catalogue table'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      ENDIF
C
C *** read the descriptor data
      CALL CATDRD(TIDCAT,NST,NSF,NSPR,RMA,RME,NIX,NIY,NPL,NL,
     2                   XCEN,YCEN,DICI,DICS,IXW,IYW)
C
      IROW = 1
      CALL CATTRD(TIDCAT,IROW,CUU,VVV(1),VVV(2),VVV(3),VVV(4),
     2                        VVV(5),VVV(6),VVV(7),VVV(8),VVV(9),
     3                        VVV(10),VVV(11),VVV(12),VVV(13))
      SIGMA = VVV(9)
      BETA  = VVV(10)
      SATU  = VVV(11)
      AIN   = VVV(12)
      SIV   = VVV(13)
C
C *** get the intermediate tabel
      CALL STKRDC('OUT_A',1,1,60,IAV,INTFIL,KUN,KNUL,ISTAT)  ! get interm. name
      CALL TBTOPN(INTFIL,F_IO_MODE,TIDINT,ISTAT)
      IF (ISTAT.NE.0) THEN                                         ! new table
         IROW   = 1
         NGRP   = 0
         IDNGRP = 0
         NOBJ   = 0
         CALL INTINI(INTFIL,TIDINT)

      ELSE                                                  ! append old table 
         CALL TBIGET(TIDINT,NCINT,NRINT,NSINT,NACINT,NARINT,ISTAT)
         IF (ISTAT.NE.0) THEN
            STRING = '*** FATAL: Problems with getting info '//
     2               'for intermediate table ...'
            CALL STTPUT(STRING,ISTAT)
            CALL STSEPI
         ENDIF

         IF (NRINT.EQ.0) THEN
            STRING = '*** INFO: No points in the intermediate '//
     2               'table'
            CALL STTPUT(STRING,ISTAT)
            IROW   = 1
            NGRP   = 0
            IDNGRP = 0
            NOBJ   = 0
         ELSE

C *** get the intermediate table parameters
            CALL INTDRD(TIDINT,NGRP,NOBJ,IDUM3,RDUM1,RDUM2,RDUM3,
     2                                   RDUM4,RDUM5,RDUM6,RDUM7)
            IROW   = NRINT + 1
            IDNGRP = NGRP
         ENDIF
      ENDIF
      CALL STECNT('PUT',EC,ED,EL)
C
C *** get the threshold
      SSLA = 2.*SIGMA*SQRT(2.**(1/BETA)-1.)
      DWI  = SSLA
      CALL STKRDR('INPUTR',1,1,IAV,SUV,KUN,KNUL,ISTAT)
      IF (SUV.GT.0.) THEN
         SIV = SUV
      ENDIF
C      IF (SIV.LT.1.) THEN
C         SIV = 1.
C      ENDIF
      AM1 = -2.5*ALOG10(SIV)
      AM2 = -2.5*ALOG10(SATU)
C
C *** get maximum number of objects to group
      CALL STKRDI('INPUTI',1,1,IAV,NT,KUN,KNUL,ISTAT)
C
C *** get the action radius
      CALL STKRDR('INPUTR',2,2,IAV,RINPUT,KUN,KNUL,ISTAT)
      SMRA   = RINPUT(1)
      ZMRA   = RINPUT(2)
C
C *** get the window percentage
      CALL STKRDR('INPUTR',4,1,IAV,PFA,KUN,KNUL,ISTAT)
      PFA = PFA/100
      IF (SMRA.GT.ZMRA) THEN
         KAA  = SMRA
         SMRA = ZMRA
         ZMRA = KAA
      END IF
C
      DO 40 KN=1,NST
         CALL CATTRD(TIDCAT,KN,NOS(KN),DAT(2,KN),DAT(3,KN),DAT(4,KN),
     2               DAT(5,KN),DAT(6,KN),ZZ,ZZ,DAT(7,KN),
     3               RDUM9,RDUM10,RDUM11,RDUM12,RDUM13)
         DAT(1,KN) = KN
   40 CONTINUE
 
C
      IF (NT.EQ.0) THEN
         NT = NTU
      ENDIF
      IF (NT.GT.NTM .OR. NT.LT.0) THEN
         NT = NTU
      ENDIF
      NTC = NT
      NO  = (NCD-(6.5*FLOAT(NT)))/3.
      DO 50 I3 = 1,40
        CHB(I3:I3) = 'B'
        CHG(I3:I3) = 'G'
   50 CONTINUE
      DO 60 I3 = 1,NTM
         FLGCMP(I3) = 1
         FLGHOL(I3) = 1
   60 CONTINUE
C
C *** write the info 
      WRITE(STRING,9001) INT(SMRA),INT(ZMRA)
      CALL STTPUT(STRING,ISTAT)
      WRITE(STRING,9002) LX0,LY0,LX1,LY1
      CALL STTPUT(STRING,ISTAT)
      WRITE(STRING,9003) SIGMA,BETA
      CALL STTPUT(STRING,ISTAT)
      WRITE(STRING,9004) SATU,SIV
      CALL STTPUT(STRING,ISTAT)
      WRITE(STRING,9005) NT,NO
      CALL STTPUT(STRING,ISTAT)
      WRITE(STRING,9007)
      CALL STTPUT(STRING,ISTAT)

C
C *** start the automatic selection
      MQQ2 = 0
      MQQ3 = 0
      MQQ4 = 0
      MQQ5 = 0
      MQQ6 = 0
      MQQ7 = 0
      MQQ8 = 0
      MQQ9 = 0
      ISAV = 0
      NEX  = 0
      NRC  = 0
      LFI  = 0
      I    = 0
      ICA  = 0
      NCP  = 0
      NHL  = 0

  100 CONTINUE
        IF (I.GT.NST) GO TO 101
           IF (ICA.GE.3.AND.ICA.LE.7) THEN
              FAR = FAR - 0.0999999
              IF (FAR.LT.SMRA/100.) THEN
                 ICA = 8
              ELSE
                 ICA = 0
                 NC  = 2
              END IF
           ELSE
              I   = I+1
              FAR = ZMRA/100.
              ICA = 0
              NC  = 2
           END IF
C
           IF (ICA.EQ.0) THEN        ! legge la stella da esaminare (record i)
              DO 110 LI = 2,7
                 V(LI) = DAT(LI,I)
  110         CONTINUE
C
              IF (MOD(I,100).EQ.0 .AND. FAR.EQ.ZMRA/100.
     2                            .AND. IDNGRP.GT.NGRP) THEN
                 WRITE(STRING,9006) IDNGRP,NCP,NHL
                 CALL STTPUT(STRING,ISTAT)
              ENDIF
C
              IF (V(7).NE.2) THEN
                 IF (LFI.EQ.0) THEN
                    IF (V(7).EQ.0.OR.V(7).EQ.-1) THEN
                       LFI = 1
                       ITF = V(7)
                       ITT = IABS(ITF)-1
                    END IF
                 END IF
                 DO 120 LI = 1,7
                    SF(1,LI) = DAT(LI,I)
  120            CONTINUE
                 LOG = V(2).GE.LX0 .AND. V(2).LT.LX0+LX1
                 LOG = V(3).GE.LY0 .AND. V(3).LT.LY0+LY1 .AND. LOG
                 LOG = COO.EQ.'A'.OR.LOG
C
                 IF ((V(7).EQ.ITF .OR. V(7).EQ.9) .AND.
     2             V(4).GE.SIV.AND.LOG) THEN
                   ASSIGN 10011 TO CERINT
                   GO TO 10010
10011              CONTINUE
                   IF (ICA.EQ.0) THEN
                     DO 130 JIL = 1,NC-1
                        IF (SF(JIL,7).EQ.ITT) ICA = 4
  130                CONTINUE
                     IF (ICA.EQ.0) THEN
                       ASSIGN 10021 TO CERSTE
                       GO TO 10020
10021                  CONTINUE
                       IF (ICA.EQ.0) THEN
                         IDNGRP     = IDNGRP + 1
                         IGRP       = IGRP   + 1
                         IOBJ       = IOBJ   + NC - 1
C
C *** prepare the output for the intermediate table
                         NCP        = NC-1
                         NHL        = KC
C
                         PARINT(1)  = FLOAT(IX)
                         PARINT(2)  = FLOAT(IY)
                         PARINT(3)  = 0.0
                         PARINT(4)  = 0.0
                         PARINT(5)  = 0.0
                         PARINT(6)  = FLOAT(LX)
                         PARINT(7)  = FLOAT(LY)
                         PARINT(8)  = 0.0
                         PARINT(9)  = 0.0
                         PARINT(10) = FOME
                         PARINT(11) = BETA
                         PARINT(12) = 0.0
                         PARINT(13) = 0.0
                         PARINT(14) = 0.0
                         PARINT(15) = FLOAT(NCP)
                         PARINT(16) = FLOAT(NHL)
C
                         DO 140 J = 1,NCP
                           FITCMP((J-1)*6+1) = SF(J,4)
                           FITCMP((J-1)*6+2) = SF(J,2)-FLOAT(IX)+1
                           FITCMP((J-1)*6+3) = SF(J,3)-FLOAT(IY)+1
                           FITCMP((J-1)*6+4) = SIGMA 
                           FITCMP((J-1)*6+5) = 0.0
                           FITCMP((J-1)*6+6) = 0.0
                           IDNGRP            = NOS(I)
                           IDNCMP(J)         = 100+J
  140                    CONTINUE

                         IF (NHL.GT.0) THEN
                           DO 160 J = 1, NHL
                             FITHOL((J-1)*3+1) = IS(J,3)
                             FITHOL((J-1)*3+2) = IS(J,1)-FLOAT(IX)+1
                             FITHOL((J-1)*3+3) = IS(J,2)-FLOAT(IY)+1
                             IDNGRP            = NOS(I)
                             IDNHOL(J)         = 200+J
  160                      CONTINUE
                         ENDIF
C
C *** write the info into the intermediate table !!!
                         CALL INTWWR(TIDINT,IROW,NCP,NHL)
                         IROW = IROW + NCP + NHL
C
                         DO 165 JIL = 1,NC-1
                            HV = -2.5*ALOG10(SF(JIL,4))
                            IF (HV.LT.AM2) HV = AM2
                            IF (HV.GT.AM1) HV = AM1
                            NCAL       = (AM1-HV)/VPA+1
                            ISTB(NCAL) = ISTB(NCAL)+1
  165                    CONTINUE
                         ISTF(NC-1)=ISTF(NC-1)+1
                       END IF
                     END IF
                   END IF

                 ELSE   ! Prim. gia' esam. o sotto la soglia o fuori finestra
                    ICA=2
                 END IF
              ELSE
                 ICA=10
              END IF
           END IF

           IF (ICA.EQ.10) THEN
              MQQ2 = MQQ2+1
           ELSE IF (ICA.EQ.3.AND.ISAV.NE.I) THEN
              MQQ3 = MQQ3+1
              ISAV = I
           ELSE IF (ICA.EQ.4.AND.ISAV.NE.I) THEN
              MQQ4 = MQQ4+1
              ISAV = I
           ELSE IF (ICA.EQ.5.AND.ISAV.NE.I) THEN
              MQQ5 = MQQ5+1
              ISAV = I
           ELSE IF (ICA.EQ.6.AND.ISAV.NE.I) THEN
              MQQ6 = MQQ6+1
              ISAV = I
           ELSE IF (ICA.EQ.7.AND.ISAV.NE.I) THEN
              MQQ7 = MQQ7+1
              ISAV = I
           ELSE IF (ICA.EQ.8) THEN
              MQQ8 = MQQ8+1
           ELSE IF (ICA.EQ.9) THEN
              MQQ9 = MQQ9+1
           END IF
C        
           IF ((ICA.LT.3.OR.ICA.GT.7).AND.ICA.NE.10) THEN
              DO 170 L = 1,NC-1
                 IL = SF(L,1)
                 DO 171 JH = 2,7
                    VVV(JH+1) = DAT(JH,IL)
  171            CONTINUE
                 LOG  =  VVV(3).GE.LX0 .AND. VVV(3).LT.LX0+LX1
                 LOG  =  VVV(4).GE.LY0 .AND. VVV(4).LT.LY0+LY1 .AND.
     2                   LOG
                 LOG  =  COO.EQ.'A' .OR. LOG
                 LOG1 = (VVV(8).EQ.ITF .OR. VVV(8).EQ.9) .AND.
     2                   VVV(5).GE.SIV
                 IF (LOG1.AND.LOG) NEX = NEX+1
                 IF (ICA.NE.0 .AND. VVV(8).EQ.ITF .AND. ICA.NE.2) THEN
                    IF (VVV(5).GE.SIV .AND. LOG) THEN
                       HV = -2.5*ALOG10(VVV(5))
                       IF (HV.LT.AM2) HV = AM2
                       IF (HV.GT.AM1) HV = AM1
                       NCAL       = (AM1-HV)/VPA+1
                       ISTC(NCAL) = ISTC(NCAL)+1
                       NRC        = NRC+1
                    END IF
                 END IF
                 IF (VVV(8).EQ.ITF .OR. VVV(8).EQ.9) THEN
                    IF (ICA.EQ.0) THEN
                       VVV(8) = 1
                    ELSE
                       IF (LOG) THEN
                          VVV(8) = ITT
                       ELSE
                          VVV(8) = 9
                       END IF
                    END IF
                 END IF
                 DO 172 JH = 2,7
                    DAT(JH,IL) = VVV(JH+1)
  172            CONTINUE
  170         CONTINUE
           END IF
           GO TO 100
  101   CONTINUE

        DO 210 JIL = 1,NST
           CALL CATTRD(TIDCAT,JIL,CUU,VVV(3),VVV(4),VVV(5),VVV(6),
     2                 VVV(7),VVV(8),VVV(9),VVV(10),VVV(11),VVV(12),
     3                 VVV(13),VVV(14),VVV(15))
           IF (DAT(7,JIL).EQ.ITF) THEN
              DAT(7,JIL) = ITT
           ENDIF

           DO 220 JH = 2,7
              VVV(JH+1) = DAT(JH,JIL)
  220      CONTINUE
           VVV(10) = VVV(8)
           CALL CATTWR(TIDCAT,JIL,CUU,VVV(3),VVV(4),VVV(5),VVV(6),
     2                 VVV(7),VVV(8),VVV(9),VVV(10),VVV(11),VVV(12),
     3                 VVV(13),VVV(14),VVV(15))
  210   CONTINUE

        IF (IOBJ.GT.0) THEN
           NGRP = NGRP + IGRP
           NOBJ = NOBJ + IOBJ
           RDUM = 0.0
           CALL INTDWR(TIDINT,NGRP,NOBJ,IDUM3,SATU,RDUM3,SIGMA,
     2                        BETA,SIV,AIN,FOG)
        END IF

        MAF = -32000
        MAC = MAF
        DO 230 JIL = 1,NT
           MAF = MAX0(MAF,ISTF(JIL))
  230   CONTINUE
        NCAL = (AM1-AM2)/VPA+1
        DO 240 JIL = 1,NCAL
           MAC = MAX0(MAC,ISTC(JIL))
  240   CONTINUE
C
        IF (MAF.GT.0.) THEN
           WRITE(STRING,9007)
           CALL STTPUT(STRING,ISTAT)
           WRITE(STRING,9008)
           CALL STTPUT(STRING,ISTAT)
           WRITE(STRING,9009)
           CALL STTPUT(STRING,ISTAT)
C
           DO 250 JIL = 1,NT
              INF = ISTF(JIL)*40/MAF
              IF (INF.GT.0) THEN
                 XXX = ' '
                 DO 251 JI = 1,INF
                    XXX(JI:JI) = 'X'
  251            CONTINUE            
                 WRITE(STRING,9010) ISTF(JIL),JIL,XXX(1:40)
                 CALL STTPUT(STRING,ISTAT)
              ELSE
                 WRITE(STRING,9011) ISTF(JIL),JIL
                 CALL STTPUT(STRING,ISTAT)
              END IF
  250      CONTINUE
        ELSE
           STRING = 'No groups found in the data set'
           CALL STTPUT(STRING,ISTAT)
        END IF
C
C       IF (MAC.GT.0.) THEN
C          WRITE(STRING,9007)
C          CALL STTPUT(STRING,ISTAT)
C          WRITE(STRING,9012) NRC
C          CALL STTPUT(STRING,ISTAT)
C          WRITE(STRING,9013)
C          CALL STTPUT(STRING,ISTAT)

C          DO 260 JIL = 1,NCAL
C             INF  = ISTC(JIL)*40/MAC
C             ANCL = -JIL*VPA+VPA
C             HHS  = 10.**(-.4*(ANCL+AM1))
C
C             IF (INF.GT.0) THEN
C                XXX = ' '
C                DO 261 JI = 1,INF
C                   XXX(JI:JI) = 'X'
C 261            CONTINUE
C                IPERCE = (100.*ISTC(JIL))/(ISTC(JIL)+ISTB(JIL))
C                WRITE(STRING,9014) ISTC(JIL),IPERCE,ANCL+AM1,HHS,
C    2                              XXX(1:40)
C             ELSE
C                IF (ISTC(JIL)+ISTB(JIL).NE.0) THEN
C                   IPERCE = (100*ISTC(JIL))/(ISTC(JIL)+ISTB(JIL))
C                ELSE
C                   IPERCE = 0
C                END IF
C                WRITE(STRING,9015) ISTC(JIL),IPERCE,ANCL+AM1,HHS
C             END IF
C             CALL STTPUT(STRING,ISTAT)
C 260      CONTINUE
C       END IF

        MAT  = -32000
        NCAL = (AM1-AM2)/VPA+1
        DO 270 JIL = 1,NCAL
           MAT = MAX0(ISTC(JIL)+ISTB(JIL),MAT)
  270   CONTINUE

        IF (MAT.GT.0) THEN
           WRITE(STRING,9007)
           CALL STTPUT(STRING,ISTAT)
           WRITE(STRING,9012) NRC
           CALL STTPUT(STRING,ISTAT)
           WRITE(STRING,9013)
           CALL STTPUT(STRING,ISTAT)

           DO 280 JIL = 1,NCAL
              INF  = ((ISTC(JIL)+ISTB(JIL))*40.)/MAT
              INFC = (ISTC(JIL)*40.)/MAT+0.5
              INFB = INF-INFC
              ANCL = -JIL*VPA+VPA
              HHS  = 10.**(-.4*(ANCL+AM1))
              IF (INF.GT.0) THEN
                 IPERCE = (100*ISTC(JIL))/(ISTC(JIL)+ISTB(JIL))
                 IF (INFB.GT.0 .AND. INFC.GT.0) THEN
                    STRING = ' '
                    WRITE(STRING,9015) ISTC(JIL),IPERCE,ANCL+AM1,HHS
                    INFT = 36  + INFB + INFC
                    STRING(37:36+INFB)   = CHG(1:INFB)
                    STRING(37+INFB:INFT) = CHB(1:INFC)
                    IF (INFT.GT.79) THEN
                       CALL STTPUT(STRING,ISTAT)
                    ELSE
                       TEXT(1:79) = STRING(1:79)
                       CALL STTPUT(TEXT,ISTAT)
                    ENDIF
 
                 ELSE IF(INFC.GT.0) THEN
                    STRING = ' '
                    WRITE(STRING,9017) 
                    LMW = INFC + 37
                    CALL STTPUT(STRING(1:LMW),ISTAT)
                 ELSE
                    STRING = ' '
                    WRITE(STRING,9018) ISTC(JIL),IPERCE,ANCL+AM1,HHS
                    LMW = INFB + 37
                    CALL STTPUT(STRING(1:LMW),ISTAT)
                 END IF

              ELSE
                 IF (ISTC(JIL)+ISTB(JIL).NE.0) THEN
                    IPERCE = (100.*ISTC(JIL))/(ISTC(JIL)+ISTB(JIL))
                 ELSE
                    IPERCE = 0
                 END IF
                 STRING = ' '
                 WRITE(STRING,9015) ISTC(JIL),IPERCE,ANCL+AM1,HHS
                 CALL STTPUT(STRING(1:40),ISTAT)
              END IF
  280      CONTINUE
        END IF

        WRITE(STRING,9007)
        CALL STTPUT(STRING,ISTAT)
        NSMQ = NST-MQQ2
        WRITE(STRING,9019) NEX,NSMQ,MQQ2
        CALL STTPUT(STRING,ISTAT)
        WRITE(STRING,9020) IOBJ
        CALL STTPUT(STRING,ISTAT)
        WRITE(STRING,9021) IGRP
        CALL STTPUT(STRING,ISTAT)

        CALL TBSINI(TIDINT,ISTAT)
        CALL TBTCLO(TIDINT,ISTAT)
        CALL TBTCLO(TIDCAT,ISTAT)
        CALL STSEPI

C +++
C   Procedure CERINT
C---
10010   CONTINUE

        NI  = 1
        NT  = NTC
        KSR = 0
        RA  = SF(1,6)
        ICA = 0
        CALL INTERS(SF,NTM,NTC,NC,NI,IFNM,RMA,NST,FAR,DAT)
        IF(IFNM.EQ.0) THEN  ! Trovate troppe compagne della primaria che non si
           ICA = 3          ! risolve neppure diminuendo il suo raggio d'azione
        END IF

        IF (ICA.EQ.0) THEN
           SF(1,6) = RA
           KS      = 1
           IF (NC.GT.2) THEN
10110         CONTINUE
              IF (NC-1.EQ.NI .OR. KS.NE.1) GO TO 10210
              NI = NI+1
              CALL INTERS(SF,NTM,NTC,NC,NI,IFNM,RMA,NST,FAR,DAT)
              IF (IFNM.EQ.0) THEN    ! Stelle che si intersecano superano il 
                 KS  = -1                             ! limite imposto da NT
                 ICA = 5
              END IF
              GO TO 10110
10210         CONTINUE
           END IF
        END IF
        GO TO CERINT

C +++
C   Procedure CERSTE
C---
10020   CONTINUE

        IXMA = -32000
        IXMI = -IXMA
        IYMA = IXMA
        IYMI = IXMI
        FOME = 0.

        DO 10120 J = 1,NC-1
           IXMA = MAX0(IXMA,INT(SF(J,2)+SSLA))
           IXMI = MIN0(IXMI,INT(SF(J,2)-SSLA))
           IYMA = MAX0(IYMA,INT(SF(J,3)+SSLA))
           IYMI = MIN0(IYMI,INT(SF(J,3)-SSLA))
           FOME = FOME+SF(J,5)
10120   CONTINUE
C 
        FOME = FOME/FLOAT(NC-1)
        DEX  = IXMA-IXMI+1
        DEY  = IYMA-IYMI+1
C
        IF (DEX.GT.NDX .OR. DEY.GT.NDY) THEN   ! Matrice per il fit supera i 
           ICA = 6                             ! limiti imposti da NDX,NDY
        ELSE
C          Vengono calcolate le dimensioni della matrice di dati
C          stella isolata : e' pari a 6 volte la larghezza della
C          PSF a meta altezza.
C          stelle multiple: e' pari alla somma fra la distanza delle
C          Stelle piu' estreme della finestra e 4 volte la larghezza
C          a meta' altezza della PSF scelta.
C
           DW2 = DWI
           IF (NC-1.EQ.1) THEN
              DW2 = 3.5*SSLA
           ELSE
              IF (NC-1 .LE. 6) THEN
                 DW2 = 2.*SSLA
              ENDIF
           ENDIF

           LX = DEX+1+2*DW2
           LY = DEY+1+2*DW2
           LX = LX*PFA
           LY = LY*PFA
           LX = MIN0(LX,NDX)
           LY = MIN0(LY,NDY)
           IX = IXMI-(LX-DEX)/2.
           IY = IYMI-(LY-DEY)/2.

           IF (IX.LT.NIX) THEN
              LX = LX+IX-NIX
              IX = NIX
           END IF
           IF (IY.LT.NIY) THEN
              LY = LY+IY-NIY
              IY = NIY
           END IF
 
           IF (IX+LX-1 .GT. NPL+NIX-1) LX = NPL+NIX-IX
           IF (IY+LY-1 .GT. NL+NIY-1)  LY = NL+NIY-IY
           INX = IX-RMA
           INY  = IY-RMA
           IFX  = IX+LX-1+RMA+.5
           IFY  = IY+LY-1+RMA+.5
           V(3) = SF(1,3)
           J    = SF(1,1)
           KC   = 0
           NBW  = (NCD-(6.5*(NC-1)))/3.

10220      CONTINUE
           IF (J.GE.NST .OR. V(3).GT.IFY .OR. ICA.NE.0) GO TO 10320
              J = J+1
              CALL COSTEL(K,IFL,NC,J,SF,NTM)
C
              IF (IFL.EQ.0) THEN
                 IF (KC.LT.NBW) THEN
                    DO 10420 K = 2,7
                       V(K) = DAT(K,J)
10420               CONTINUE
                    IF (V(3).LE.IFY) THEN
                       CALL LIMITX(V,IS,NOB,KC,IX,IY,LX,LY)
                    END IF
                 ELSE     !Il numero dei buchi supera il limite imposto da NBW
                    ICA = 7
                 END IF
              END IF
              GO TO 10220
C
10320      CONTINUE
           V(3) = SF(1,3)
           J    = SF(1,1)
 
10520      CONTINUE
           IF (J.LE.1 .OR. V(3).LT.INY .OR. ICA.NE.0) GO TO 10620
              J = J-1
              CALL COSTEL(K,IFL,NC,J,SF,NTM)
 
              IF (IFL.EQ.0) THEN
                 IF (KC.LT.NBW) THEN
                    DO 10521 K = 2,7
                       V(K) = DAT(K,J)
10521               CONTINUE
                    IF (V(3).GE.INY) THEN
                       CALL LIMITX(V,IS,NOB,KC,IX,IY,LX,LY)
                    END IF
                 ELSE       !Il numero dei buchi supera il limite imposto da NBW
                    ICA = 7
                 END IF
              END IF
              GO TO 10520
10620      CONTINUE
        END IF

        JIL = 0
10720   CONTINUE
        IF (JIL.GE.KC .OR. ICA.NE.0) GO TO 10820
           JIL = JIL+1
           RA  = IS(JIL,3)**2
           JEL = 0
10920      CONTINUE
           IF (JEL.GE.(NC-1) .OR. ICA.NE.0) GO TO 11020
              JEL = JEL+1
              RAS = SF(JEL,6)**2
              DIS = (IS(JIL,1)-SF(JEL,2))**2+(IS(JIL,2)-SF(JEL,3))**2
              IF (DIS.LE.RA) THEN
                 IF (IS(JIL,4).EQ.2) THEN
                    IS(JIL,3) = SQRT(DIS)-1.
                 ELSE
                    ICA = 9
                 END IF
              END IF
              IF (DIS.LE.RAS) THEN
                 IF (IS(JIL,4).NE.2) THEN 
                    ICA = 9
                 END IF
              END IF
              GO TO 10920
11020      CONTINUE
           GO TO 10720
10820   CONTINUE

        GO TO CERSTE
        END
