C @(#)rfotsearch.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 SEARCH   
C+++
C.IDENTIFICATION: RFOTSEARCH
C.PURPOSE: Search for objects above a certain threshold in a frame
C.AUTHOR:  R. Buonanno, G. Buscema, C. Corsi, I. Ferraro, G. Iannicola
C          Osservatorio Astronomico di Roma
C.VERSION: 870906 RXB First running version at ESO (outside MIDAS)
C          871217 RHW Installation in MIDAS
C          880929 RHW New version
C          880104 RHW Include MIDAS for background values
C          8901?? RHW Rewritten for the portable MIDAS version 
C                     Inclusion of MIDAS tables 
C.VERSION: 910122 RHW IMPLICIT NONE added; all variables defined
C----
       IMPLICIT    NONE
C
       INTEGER     NX
       INTEGER     NDSAT
       INTEGER     IDFS
       PARAMETER   (NX=10000)
       PARAMETER   (NDSAT=50000)
       PARAMETER   (IDFS=1000)

       DOUBLE PRECISION BEGIN(2),STEP(2)

       INTEGER      AREA(6)
       INTEGER      EC, ED, EL
       INTEGER      I, I2
       INTEGER      IXW, IYW
       INTEGER      IX0, IY0
       INTEGER      IXR, IYR
       INTEGER      IST(2000)
       INTEGER      ISTAT,IAC
       INTEGER*8    IPNTR
       INTEGER      IMF
       INTEGER      IACOL, IAROW
       INTEGER      IDX, IDY
       INTEGER      IMX, IMY
       INTEGER      I1, IB, IA1, IA5, IA6
       INTEGER      IM, IK, IW, INL, ILA
       INTEGER      IDUM, IGRP
       INTEGER      J, J1
       INTEGER      KUN, KNUL
       INTEGER      KOSTE, KLIN, KSAT
       INTEGER      KR, KN, KR1, KK
       INTEGER      LN(3)
       INTEGER      LK, LF
       INTEGER      L, LL
       INTEGER      MADRID(1)
       INTEGER      MTS(NDSAT,4)
       INTEGER      NLN(3)
       INTEGER      NPIX(2)
       INTEGER      NAXIS, NPL, NL
       INTEGER      NLI, NI
       INTEGER      NCX, NCY
       INTEGER      NDUM, NCO, NCHAR, NTRY
       INTEGER      NCSKY, NRSKY, NN, NSF
       INTEGER      N3, NYM, NY2
       INTEGER      SKYCOL(256)
       INTEGER      TIDCAT,TIDSKY
C       
       REAL         A, A1, A2, A3, A4, A5, A6, A7
       REAL         AIN, AB1, RA
       REAL         BETA
       REAL         BACKGR(256)
       REAL         DICI, DICS
       REAL         FON(IDFS)
       REAL         FLB, FONM, FL, F
       REAL         H 
       REAL         P(5000,6)
       REAL         PRX, PRY, PR, PESO
       REAL         RINPUT(3)
       REAL         RIJ, RISI, RAG, R
       REAL         RAMA, RAME
       REAL         SGL(IDFS)
       REAL         SINC, SINR
       REAL         SIGMA, SAT, SUPE
       REAL         XCEN, YCEN
       REAL         XBAR, YBAR
       REAL         X,Y
       REAL         ZMA(NX,3)
C
       DOUBLE PRECISION    DDUM
C
       LOGICAL      FDX,FDY,VLO,NULL(256)
C
       CHARACTER    IDENT*72,CUNIT*80
       CHARACTER    FRAME*60, SKYFIL*60, CATFIL*60
       CHARACTER    SPF*1
       CHARACTER    STRING*80, CVALUE*80
C	
       INCLUDE      'MID_INCLUDE:ST_DEF.INC/NOLIST'
       COMMON      /PRONE/KR,DICI,DICS,RIJ,ZMA,RISI,IXW,IYW,SGL,FON,
     *                    SAT,XCEN,YCEN,IX0,IY0,IM,NCX,NCY,IDX,IDY,
     *                    AIN,RAMA,RAME,FONM,MTS,KSAT,I1,NYM,LN,
     *                    BETA,SIGMA,KN,PR,NY2
       COMMON       /VMR/MADRID
       INCLUDE      'MID_INCLUDE:ST_DAT.INC/NOLIST'
C
       DATA         NLN/2,3,1/
C
 9001  FORMAT('*** INFO: Area in pixel units: [', 
     2            I5, ',', I5, ':', I5, ',', I5, ']')
 9010  FORMAT('          Minimum threshold for star detection:',E12.5)
 9020  FORMAT('          Number of objects found: ',I6,
     2        '; saturated: ',I6)
C
C *** START 
       CALL STSPRO('SEARCH')
       LN(1) = 1
       LN(2) = 2
       LN(3) = 3
       KSAT  = 0
C
C *** input parameters
       CALL STKRDC('IN_A',1,1,60,IAC,FRAME,KUN,KNUL,ISTAT)      ! name of image
       CALL STKRDC('IN_B',1,1,60,IAC,SKYFIL,KUN,KNUL,ISTAT)  ! sky backgr table
       CALL STKRDC('OUTPUTC',1,1,60,IAC,CATFIL,KUN,KNUL,ISTAT)  ! catalog table
C
       CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     2             2,NAXIS,NPIX,BEGIN,STEP,IDENT,CUNIT,IPNTR,IMF,ISTAT)
       NPL = NPIX(1)
       NL  = NPIX(2)
       IXW = BEGIN(1)
       IYW = BEGIN(2)
C
C *** read the background values from the sky table
      CALL STECNT('GET',EC,ED,EL)
      CALL STECNT('PUT',1,0,0)
      CALL TBTOPN(SKYFIL,F_I_MODE,TIDSKY,ISTAT)
      IF (ISTAT.NE.0) THEN
         STRING = '*** FATAL: problems with opening the sky '//
     2            'background table; try again ...'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      ENDIF
C
      CALL STDRDI(IMF,'SKY_AREA',1,6,IAC,AREA,KUN,KNUL,ISTAT)
      IF (ISTAT.NE.0) THEN
         STRING = '*** FATAL: No information about the background'//
     2            ' determination'
         CALL STTPUT(STRING,ISTAT)
         STRING = '           Did you run SKY/ROMAFOT before ?'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      ENDIF
      CALL STECNT('PUT',EC,ED,EL)
C
      CALL STKRDC('INPUTC',1,1,40,IAC,STRING,KUN,KNUL,ISTAT)        ! user area
      IF (STRING(1:3).EQ.'SKY') THEN
         IX0  = AREA(1)
         IY0  = AREA(2)
         IDX  = AREA(3) - AREA(1) + 1
         IDY  = AREA(4) - AREA(2) + 1
      ELSE
         CALL EXTCOO(IMF,STRING,NAXIS,NDUM,AREA(1),AREA(3),ISTAT)
         IX0  = AREA(1)
         IY0  = AREA(2)
         IDX  = AREA(3) - AREA(1) + 1
         IDY  = AREA(4) - AREA(2) + 1
      ENDIF

      NCX = AREA(5)
      NCY = AREA(6)
      NCO = NCX*NCY
      IF (IX0.LT.1) IX0 = 1
      IF (IY0.LT.1) IY0 = 1
      IF ((IX0+IDX-1).GT.NPL) IDX = NPL-IX0+1
      IF ((IY0+IDY-1).GT.NL)  IDY = NL-IY0+1
C
      WRITE (STRING,9001) AREA(1),AREA(2),AREA(3),AREA(4)
      CALL STTPUT(STRING,ISTAT)
C
C *** get the PSF values
      CALL STKRDR('INPUTR',1,3,IAC,RINPUT,KUN,KNUL,ISTAT)
      SIGMA = RINPUT(1)
      BETA  = RINPUT(2)
      SAT   = RINPUT(3)
C
C *** get photometric option, threshold
      STRING = ' '
      CALL STKRDC('INPUTC',1,41,40,IAC,STRING,KUN,KNUL,ISTAT)
      CALL UPCAS(STRING,STRING)
      IF (STRING(1:3).EQ.'PRO') THEN
         CALL STKPRC('Enter photom. option and threshold: ',
     2               'INPUTC',1,41,40,IAC,STRING,KUN,KNUL,ISTAT)
      ENDIF
      SPF    = STRING(1:1)
      NCHAR  = INDEX(STRING,',')
      CVALUE = STRING(NCHAR+1:)
      CALL GENCNV(CVALUE,2,1,IDUM,RINPUT,DDUM,ISTAT)
C     CALL USRINP(RINPUT,1,'R',CVALUE)       
      IF (SPF.EQ.'C') THEN                       !constant photometric treshold
         SINC = RINPUT(1)
      ELSE IF (SPF.EQ.'R') THEN                      !N times the Sky Stan. Dev
         SINR = INT(RINPUT(1))
      ELSE
         STRING = '*** FATAL: Unknown option; try again ...'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      END IF
C
C *** min. height to determine influence of star
      CALL STKRDR('INPUTR',4,1,IAC,AIN,KUN,KNUL,ISTAT) 
      NTRY = 0
   12 CONTINUE
      IF (AIN.EQ.0 .AND. NTRY.EQ.0) THEN
         CALL STKPRR('Enter minimum height: ',
     2               'INPUTR',4,1,IAC,AIN,KUN,KNUL,ISTAT)
         NTRY = NTRY + 1
         GO TO 12

      ELSE IF (AIN.EQ.0 .AND. NTRY.GE.1) THEN
         STRING = '*** WARNING: Minimum can not be 0.0; try again ...'
         CALL STTPUT(STRING,ISTAT)
         CALL STKPRR('Enter minimum height: ',
     2               'INPUTR',4,1,IAC,AIN,KUN,KNUL,ISTAT)
         NTRY = NTRY + 1
         GO TO 12
      ENDIF
C
C *** get the info from the sky table
      CALL TBIGET(TIDSKY, NCSKY, NRSKY, NN, IACOL, IAROW, ISTAT)
      DO 19 I = 1,NCSKY
         SKYCOL(I) = I
   19 CONTINUE
      DO 20 J = 1,NRSKY
         CALL TBRRDR(TIDSKY,J,NCSKY,SKYCOL,BACKGR,NULL,ISTAT)
         DO 21 I = 1,NCSKY
            FON((J-1)*NCSKY+I) = BACKGR(I)
   21    CONTINUE
   20 CONTINUE
C
C *** create the catalogue table
      CALL CATINI(CATFIL,TIDCAT)
C
C *** get the threshold
      FONM = 2000000000.
      DO 50 LL = 1,NCO
         IF (SPF.EQ.'C' .OR. SPF.EQ.'c') THEN
            SGL(LL) = SINC
         ELSE
            SGL(LL) = SINR*SQRT(FON(LL))
         END IF
         FONM = MIN(FONM,SGL(LL))
   50 CONTINUE

      RAMA = 0.
      RAME = 0.
      NYM  = 3
      NY2  = NYM/2
      N3   = 1
      KN   = 0
      PR   = (DICS-DICI)/NCO
      WRITE (STRING,9010) FONM
      CALL STTPUT(STRING,ISTAT)
C
C *** search for objects
       DO 110 I = 1,NYM
          INL = I
          KLIN = I + IY0 - 1
          CALL REALIN(NPL,NL,KLIN,IX0,IDX,MADRID(IPNTR),ZMA(1,LN(INL)))
  110  CONTINUE
       KR = 0
       DO 120 I1 = 1,NY2+1
          CALL PROONE(TIDCAT)
  120  CONTINUE
C
       I1  = NY2+1
       INL = 3
       DO 130 I = 4,IDY
C         IF (MOD(I,100).EQ.0) THEN
C            WRITE(STRING,9025) I
C            CALL STTPUT(STRING,ISTAT)
C         END IF
C
          DO 131 L = 1,NYM
             LN(L) = NLN(LN(L))
  131     CONTINUE
          KLIN = I + IY0 - 1
          CALL REALIN(NPL,NL,KLIN,IX0,IDX,MADRID(IPNTR),ZMA(1,LN(INL)))
          CALL PROONE(TIDCAT)
  130  CONTINUE
C
       DO 140 I = 1,NYM-NY2-1
          I1 = I1+1
          CALL PROONE(TIDCAT)
  140  CONTINUE
C
       NSF  = 0.
       RAME = RAME/KN
       IXR  = IX0+IXW-1
       IYR  = IY0+IYW-1
       IA1  = 0
       AB1  = 0.
       CALL CATDWR(TIDCAT,KN,NSF,IA1,RAMA,RAME,IXR,IYR,IDX,IDY,
     2             AB1,AB1,AB1,AB1,IXW,IYW)
C
C ***
       IF (KSAT.GT.0) THEN             
          KOSTE = 0
          RA    = 0
C
          DO 210 I = 1,KSAT
             IF (MTS(I,4).EQ.0) THEN
                NLI      = 1
                MTS(I,4) = 1
                NI       = 1
                IST(NI)  = I
  211           CONTINUE
                IF (NI.GT.NLI) GO TO 212
                CALL INTSAT(MTS,NDSAT,4,NI,IST,NLI,KSAT)
                NI = NI+1
                GO TO 211
  212           CONTINUE
                KOSTE = KOSTE+1
C
C *** start procedure PARSAT
                XBAR = 0
                YBAR = 0
                SUPE = 0
     
                DO 213 I1 = 1,NLI
                   PESO = MTS(IST(I1),3)-MTS(IST(I1),1)+1
                   XBAR = XBAR+(MTS(IST(I1),3)+MTS(IST(I1),1))/2*PESO
                   YBAR = YBAR+MTS(IST(I1),2)*PESO
                   SUPE = SUPE+PESO
  213           CONTINUE

                P(KOSTE,1) = XBAR/SUPE+IX0+IXW-2
                P(KOSTE,2) = YBAR/SUPE+IY0+IYW-2

                IF (BETA.GT.0) THEN
                  P(KOSTE,3) = SAT*(1+SUPE/(3.14159*SIGMA**2))**BETA
                ELSE
                  P(KOSTE,3) = SAT*(EXP(4*ALOG(2.)*SUPE/
     2                             (3.14159*SIGMA**2)))
                END IF

                J1  = NINT(P(KOSTE,1))+0.5
                KR1 = NINT(P(KOSTE,2))+0.5
                J   = J1-IX0-IXW+2
                KR  = KR1-IY0-IYW+2
C 
C *** start procedure CALCOR
                PRX = FLOAT(IDX)/FLOAT(NCX)
                PRY = FLOAT(IDY)/FLOAT(NCY)
                IF (PRX-INT(PRX).EQ.0.) THEN
                   IMX = (J-1)/INT(PRX)+1
                ELSE
                   IMX = (J-1)/PRX+1
                END IF
                IF (PRY-INT(PRY).EQ.0.) THEN
                   IMY = (KR-1)/INT(PRY)+1
                ELSE
                   IMY = (KR-1)/PRY+1
                END IF
                IM = NCX*(IMY-1)+IMX
C *** end procedure CALCOR
C
                P(KOSTE,4) = FON(IM)
                IF (BETA.GT.0) THEN
                   RAG = SIGMA*SQRT((AIN/P(KOSTE,3))**(-1./BETA)-1.)
                ELSE
                   RAG = SIGMA*SQRT(ALOG(P(KOSTE,3)/AIN)/(4*ALOG(2.)))
                END IF
                RAMA       = AMAX1(RAMA,RAG)
                RA         = RA+RAG
                P(KOSTE,5) = RAG
                FLB        = 0
                FDX        = J.EQ.1.OR.J.EQ.IDX
                FDY        = KR.EQ.1.OR.KR.EQ.IDY
                IF (FDX.OR.FDY) FLB = 2
                P(KOSTE,6) = FLB
C *** end procedure PARSAT
C
             END IF
  210     CONTINUE
          IK = KOSTE-1
  215     CONTINUE    
          LF = 0
C
          DO 220 I = 1,IK
             VLO = P(I,2).GT.P(I+1,2)
             IF (VLO .OR. (P(I,2).EQ.P(I+1,2) .AND. 
     2          P(I,1).GT.P(I+1,1))) THEN
                DO 221 KK = 1,6
                   A         = P(I,KK)
                   P(I,KK)   = P(I+1,KK)
                   P(I+1,KK) = A
 221            CONTINUE
                LK = I
                LF = 1
             END IF
 220      CONTINUE
          IK = LK
          IF (LF.EQ.1) THEN
             GO TO 215
          ENDIF
C
C ***  ordina stelle
          IF (KOSTE.GT.0) THEN
             ILA = 0
             J   = KOSTE
             CALL CATDRD(TIDCAT,KN,NSF,IB,A,RAME,IXR,IYR,IDX,IDY,A1,
     2                   A2,A3,A4,IA5,IA6)
             I   = KN
             IW  = I+J
             WRITE(STRING,9020) IW,J 
             CALL STTPUT(STRING,IST)
C
  301        CONTINUE
             IF (I.LT.1.OR.ILA.NE.0) GO TO 300
                CALL CATTRD(TIDCAT,I,IGRP,X,Y,H,F,R,FL,A1,A2,A3,A4,
     2                      A5,A6,A7)
                I = I-1
                IF (Y.GT.P(J,2).OR.(Y.EQ.P(J,2).AND.X.GT.P(J,1))) THEN
                   IGRP = IW
                   IF (H.GT.999999.) H = 999999.
                   CALL CATTWR(TIDCAT,IW,IGRP,X,Y,H,F,R,FL,F,AB1,
     2                         SIGMA,BETA,SAT,AIN,FONM)
                   IW = IW-1

                ELSE
  302              CONTINUE
                      IF (P(J,2).LE.Y .AND. (P(J,2).NE.Y .OR. 
     2                   P(J,1).LE.X) .OR. ILA.NE.0) GO TO 303
                      IGRP = IW
                      IF (P(J,3).GT.999999.) P(J,3) = 999999.
                      CALL CATTWR(TIDCAT,IW,IGRP,P(J,1),P(J,2),P(J,3),
     2                            P(J,4),P(J,5),P(J,6),P(J,4),AB1,
     2                            SIGMA,BETA,SAT,AIN,FONM)
                      IW = IW-1
                      J  = J-1
                      IF (J.EQ.0) ILA = 1
                   GO TO 302
C
  303              CONTINUE
                   IGRP = IW
                   IF (H.GT.999999.) H = 999999.
                   CALL CATTWR(TIDCAT,IW,IGRP,X,Y,H,F,R,FL,F,AB1,
     2                         SIGMA,BETA,SAT,AIN,FONM)
                   IW = IW-1
                END IF
             GO TO 301

  300        CONTINUE
             IF (J.GT.0) THEN
                DO 310 I2 = 1,J
                   IGRP = I2
                   IF (P(I2,3).GT.999999.) P(I2,3)=999999.
                   CALL CATTWR(TIDCAT,I2,IGRP,P(I2,1),P(I2,2),P(I2,3),
     2                         P(I2,4),P(I2,5),P(I2,6),P(I2,4),AB1,
     3                         SIGMA,BETA,SAT,AIN,FONM)
  310          CONTINUE
             END IF
             RAME = (RAME*KN+RA)/(KN+KOSTE)
             KN   = KN+KOSTE
             IDUM = 0
             CALL CATDWR(TIDCAT,KN,NSF,IDUM,RAMA,RAME,IXR,IYR,IDX,IDY,
     2                          A1,A2,A3,A4,IXW,IYW)
          END IF
       END IF
C
C *** close files
       CALL TBTCLO(TIDSKY,ISTAT)
       CALL TBSINI(TIDCAT,ISTAT)
       CALL TBTCLO(TIDCAT,ISTAT)
       CALL STSEPI
       END