C @(#)rfotfind.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 FIND
C+++
C.IDENTIFICATION: RFOTFIND
C.PURPOSE: Select objects from a displayed frame and put them into a
C.         catalogue table
C.AUTHOR:  R. Buonanno, G. Buscema, C. Corsi, I. Ferraro, G. Iannicola
C          Osservatorio Astronomico di Roma
C.VERSION: 870417 RXB First running version at ESO (outside MIDAS)o
C.         881031 RXB New version                                  
C.         890228 RHW Major changes for portable MIDAS   
C.                    load of frame done by LOAD/IMAGE
C.                    cursor interaction with GETCURS
C.VERSION: 900122 RHW IMPLICIT NONR added; all variables defined
C----
      IMPLICIT         NONE
C
C     INTEGER          AREA(4)
      INTEGER          CCOUNT,COOFF
      INTEGER          ED, EC, EL
      INTEGER          IVX(9)
      INTEGER          ICMP, IROW, IDUM
      INTEGER          INX, I, IT
      INTEGER          IAV, ICOM, IAC
      INTEGER          ISTAT,ISTAT1,ISTAT2
      INTEGER*8        IPNTR
      INTEGER          IMF
      INTEGER          ICATD2, ICATD3, ICATD6, ICATD7
      INTEGER          ICATD8, ICATD9, ICATD14, ICATD15
      INTEGER          J
      INTEGER          LX, LY
      INTEGER          LX1, LY1
      INTEGER          KL
      INTEGER          KUN, KNUL
      INTEGER          MADRID(1)
      INTEGER          NPIX(3)
      INTEGER          NAXIS, NPL,NL, NC
      INTEGER          NCCAT, NRCAT, NSCAT
      INTEGER          NARCAT, NACCAT
      INTEGER          NOBJ, NLIV
      INTEGER          NCH
      INTEGER          TIDCAT
      INTEGER          XY1, XY2
C
      DOUBLE PRECISION BEGIN(3),STEP(3), DDUM
C
      REAL             AL
      REAL             BB, B1
      REAL             FPIX1(2),FPIX2(2)
      REAL             FG, FO1
      REAL             H1
      REAL             PAR(4)
      REAL             PARS(4)
      REAL             RINPUT(4)
      REAL             RCATD4, RCATD5, RCATD10
      REAL             RCATD11, RCATD12, RCATD13
      REAL             RX(100),RY(100)
      REAL             RME, VME
      REAL             SQM, SQR, SL
      REAL             UU, U1
      REAL             VALUE1,VALUE2
      REAL             VET(100)
      REAL             VV
      REAL             WCO1(2),WCO2(2)
      REAL             X1, Y1

      CHARACTER*60     FRAME
      CHARACTER*72     CUNIT
      CHARACTER*80     IDENT
      CHARACTER*60     CATFIL
      CHARACTER*3      ACTION
      CHARACTER*80     STRING
      CHARACTER*60     CINPUT
      LOGICAL          NUL
C
      INCLUDE          'MID_INCLUDE:IDIDEV.INC'
      INCLUDE          'MID_INCLUDE:IDIMEM.INC'
      INCLUDE          'MID_INCLUDE:ST_DEF.INC'
      COMMON           /VMR/MADRID
      INCLUDE          'MID_INCLUDE:ST_DAT.INC'
C
      DATA             IVX/1,2,3,4,5,6,7,8,9/
C
 9001 FORMAT('*** INFO: Catalogue table already contains', I5,
     2       ' objects')
 9003 FORMAT('Enter identification., mag., col1 and col2: [', I6,
     2       ',0.0,0.0,0.0]: ')
C
C *** start the code
      CALL STSPRO('FIND')
      CALL DTOPEN(1,ISTAT)
C
C *** input parameters
      CALL STKRDC('IN_A',1,1,60,IAC,FRAME,KUN,KNUL,ISTAT)      ! name of image
      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)
C
C *** get the intemediate file name
      CALL STKRDC('OUT_A',1,1,60,IAC,CATFIL,KUN,KNUL,ISTAT)   ! catalogue file
C
C *** open or create the calogue file 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 = '*** INFO: Catalogue table not present; '//
     2            'will create a new one ...'
         CALL STTPUT(STRING,ISTAT)
         CALL CATINI(CATFIL,TIDCAT)
         ICMP    = 0
         IROW    = 0
         ICATD2  = 0
         ICATD3  = 0
         RCATD4  = 0.0
         RCATD5  = 0.0
         ICATD6  = 0
         ICATD7  = 0
         ICATD8  = 0
         ICATD9  = 0
         RCATD10 = 0.0
         RCATD11 = 0.0
         RCATD12 = 0.0
         RCATD13 = 0.0
         ICATD14 = 0
         ICATD15 = 0

      ELSE
         CALL TBIGET(TIDCAT,NCCAT,NRCAT,NSCAT,NACCAT,NARCAT,ISTAT)
         IF (ISTAT.NE.0) THEN
            STRING = '*** FATAL: Problems with getting info for '//
     2               'catalogue table ...'
            CALL STTPUT(STRING,ISTAT)
            CALL STSEPI
         ENDIF
C
         IF (NRCAT.EQ.0) THEN
            STRING = '*** INFO: No points in the catalogue table'
            CALL STTPUT(STRING,ISTAT)
         ENDIF
C
C *** get the catalogue table parameters
         CALL CATDRD(TIDCAT,NOBJ,ICATD2,ICATD3,RCATD4,RCATD5,ICATD6,
     2                          ICATD7,ICATD8,ICATD9,RCATD10,RCATD11,
     3                          RCATD12,RCATD13,ICATD14,ICATD15)
         IROW = NRCAT
         WRITE(STRING,9001) NOBJ
         CALL STTPUT(STRING,ISTAT)
         CALL TBERDI(TIDCAT,NRCAT,1,ICMP,NUL,ISTAT)
      ENDIF
      CALL STECNT('PUT',EC,ED,EL)
C
C *** get the area
C     CALL STKRDC('INPUTC',1,1,40,IAC,STRING,KUN,KNUL,ISTAT)        ! user area
C     CALL EXTRCO(STRING,NAXIS,NPIX,BEGIN,STEP,NDUM,AREA(1),
C   2                   AREA(3),ISTAT)
C     IX0  = AREA(1)
C     IY0  = AREA(2)
C     NX   = AREA(3) - AREA(1) + 1
C     NY   = AREA(4) - AREA(2) + 1
C
C     IF (IX0.LT.1) IX0 = 1
C     IF (IY0.LT.1) IY0 = 1
C     IF ((IX0+NX-1).GT.NPL) NX = NPL-IX0+1
C     IF ((IY0+NY-1).GT.NL)  NY = NL-IY0+1
C     WRITE (STRING,9001) AREA(1),AREA(2),AREA(3),AREA(4)
C     CALL STTPUT(STRING,ISTAT)
C
C
C *** infinite loop for cursor interaction
      COOFF  = 0
      CCOUNT = 0
      ACTION = 'YYN'
      FRAME  = ' '
      CALL STTPUT('*** INFO: Use cursor control panel to '//
     2            'select object',ISTAT)

 1000 CONTINUE
      ICMP = ICMP + 1
 1001 CONTINUE
      CALL GETCUR(ACTION,FRAME,XY1,FPIX1,WCO1,VALUE1,ISTAT1,
     2                         XY2,FPIX2,WCO2,VALUE2,ISTAT2)
      IF (ISTAT1.EQ.0) THEN
         IF (CCOUNT.EQ.0) THEN
            IF (COOFF.EQ.1) THEN
               CALL DTCLOS(QDSPNO)
               CALL STSEPI
            ELSE
               CALL STTPUT('*** WARNING: Switch cursor on ...'//
     2                     ' next time we exit',ISTAT)
            ENDIF
            COOFF = 1
            GO TO 1001

         ELSE
            CALL CATDWR(TIDCAT,NOBJ,ICATD2,ICATD3,RCATD4,RCATD5,ICATD6,
     2                             ICATD7,ICATD8,ICATD9,RCATD10,RCATD11,
     3                             RCATD12,RCATD13,ICATD14,ICATD15)
            CALL TBSINI(TIDCAT,ISTAT)
            CALL TBTCLO(TIDCAT,ISTAT)
            CALL DTCLOS(QDSPNO)
            CALL STSEPI
         ENDIF

      ELSE
         CCOUNT = 1
         LX     = INT(FPIX1(1))
         LY     = INT(FPIX1(2))
C
C  *** fit data
         AL     = 4*ALOG(2.)
         NC     = 0
         DO 110 J = 1,9
            RX(J) = 0
            RY(J) = 0
  110    CONTINUE
C
         H1 = -10.E15
         DO 120 I = LY-4,LY+4
            NC  = NC+1
            INX = LX-4
            CALL REALIN(NPL,NL,I,INX,9,MADRID(IPNTR),VET)
            RME = 0
C
            DO 121 J=1,9
               H1  = AMAX1(H1,VET(J))
               RME = RME+VET(J)
               RX(J) = RX(J)+VET(J)
  121       CONTINUE
            RY(NC) = RME/9.
  120    CONTINUE
C
         PAR(1) = 0.
         PAR(4) = 10.E15
         DO 130 J = 1,9
            RX(J) = RX(J)/9.
            IF (RX(J).GT.PAR(1)) THEN
               PAR(1) = RX(J)
               PAR(2)=J
            END IF
            IF (RX(J).LT.PAR(4)) PAR(4)=RY(J)
  130    CONTINUE
C
         PAR(1) = PAR(1)-PAR(4)
         KL     = PAR(2)
         VME    = (RX(KL-1)+RX(KL+1))/2.
         IF (VME.GE.PAR(1)) VME=PAR(1)/2.
         PAR(3) = SQRT(AL/ALOG(PAR(1)/VME))
         SQR    = 10.**15
         SL     = SQR
         IT     = 0
         NLIV   = 9
C
         DO 140 J = 1,4
            PARS(J)=PAR(J)
  140    CONTINUE

 2001    CONTINUE
            IF (SL.LE.0.0001.OR.IT.GT.20) GO TO 2000
            IT = IT+1
            CALL MONO4(IVX,RX,NLIV,PAR,.7)
            IF (PAR(1).LT..1 .OR. ABS(PAR(2)-PARS(2)).GT.4 .OR.
     2         PAR(3).LE.0. .OR. PAR(3).GT.20.) THEN
               DO 150 J=1,4
                  PAR(J) = PARS(J)
  150          CONTINUE
               GO TO 2000
            ENDIF
C
            SQM = 0.
            DO 160  I = 1,NLIV
               FG  = PAR(1)*EXP(-4*ALOG(2.)*((I-PAR(2))/PAR(3))**2)
               SQM = SQM+(RX(I)-FG-PAR(4))**2
  160       CONTINUE
            SQM    = SQRT(SQM/NLIV)
            SL     = ABS(SQR-SQM)/SQR
            SQR    = SQM
           GO TO 2001
 2000    CONTINUE

         X1     = PAR(2)+LX-5
         FO1    = PAR(4)
C 
         PAR(1) = 0
         PAR(4) = 10.E15
         DO 170  J=1,9
            IF (RY(J).GT.PAR(1)) THEN
               PAR(1) = RY(J)
               PAR(2) = J
            END IF
            IF (RY(J).LT.PAR(4)) THEN
               PAR(4) = RY(J)
            ENDIF
  170    CONTINUE
C
         PAR(1) = PAR(1)-PAR(4)
         KL     = PAR(2)
         VME    = (RY(KL-1)+RY(KL+1))/2.
         IF (VME.GE.PAR(1)) VME = PAR(1)/2.
         PAR(3) = SQRT(AL/ALOG(PAR(1)/VME))
         SQR    = 10.**15
         SL     = SQR
         IT     = 0
         NLIV   = 9
C
         DO 180 J = 1,4
            PARS(J) = PAR(J)
  180    CONTINUE
C
 3001    CONTINUE
            IF (SL.LE.0.0001.OR.IT.GT.20) GO TO 3000
            IT = IT+1
            CALL MONO4(IVX,RY,NLIV,PAR,.7)
            IF (PAR(1).LT..1 .OR. ABS(PAR(2)-PARS(2)).GT.4 .OR. 
     2          PAR(3).LE.0. .OR. PAR(3).GT.20.) THEN
               DO 190 J = 1,4
                  PAR(J) = PARS(J)
  190          CONTINUE
               GO TO 3000
            END IF
C
            SQM = 0.
            DO 200 I = 1,NLIV
               FG = PAR(1)*EXP(-4*ALOG(2.)*((I-PAR(2))/PAR(3))**2)
               SQM = SQM+(RY(I)-FG-PAR(4))**2
  200       CONTINUE
C
            SQM = SQRT(SQM/NLIV)
            SL  = ABS(SQR-SQM)/SQR
            SQR = SQM
            GO TO 3001
 3000    CONTINUE
C
         Y1  = PAR(2)+LY-5
         FO1 = (FO1+PAR(4))/2.
         H1  = H1-FO1      
C
C *** legge magnitudine v , colori b-v e u-b e nome della
C     stella prescelta e li scrive sul file f
3002     CONTINUE
         WRITE (STRING,9003) ICMP
         NCH = INDEX(STRING,':')+1               
         CINPUT = ' '
         CALL STKPRC(STRING(1:NCH),'INPUTC',1,1,60,IAV,CINPUT,
     2               KUN,KNUL,ISTAT)
         ICOM   = MIN(INDEX(CINPUT,',')-1,6)
         IF (ICOM.LE.0) THEN
            VV    = 0.0
            BB    = 0.0
            UU    = 0.0
         ELSE
            ICMP = 4
            CALL GENCNV(CINPUT,2,ICMP,IDUM,RINPUT,DDUM,ISTAT)
C           CALL USRINP(RINPUT,4,'R',CINPUT)
C           ICMP = NEL(RINPUT,4)
            IF (ICMP.GT.0 .AND. ICMP.LT.4) THEN
               DO 3003 ICMP = ICMP+1, 4
                  RINPUT(ICMP) = 0.0
 3003          CONTINUE
            ENDIF
            ICMP    = INT(RINPUT(1))
            IF (ICMP.GT.99999) THEN
               CALL STTPUT('*** WARNING: Identification '//
     2                     'should be less than 100000; '//
     3                     'try again ...',ISTAT)
               GO TO 3002
            ELSE
               VV     = RINPUT(1)
               BB     = RINPUT(2)
               UU     = RINPUT(3)
            ENDIF
         ENDIF 
C
C *** prepare the arrays for writting
         IROW      = IROW + 1
         B1        = VV+BB
         U1        = B1+UU
         LX1       = BEGIN(1)+X1-1
         LY1       = BEGIN(2)+Y1-1
         NOBJ      = NOBJ+1
         CALL CATTWR(TIDCAT, IROW, ICMP, FLOAT(LX1), FLOAT(LY1),
     2                       H1, VV, B1, U1, 0.0, 0.0, 0.0, 0.0, 0.0,
     3                       0.0, 0.0)
         GO TO 1000
      ENDIF

      CALL STSEPI
      END
