C @(#)rfotselect.for	19.1 (ES0-DMD) 02/25/03 13:30:15
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 SELECT
C++
C.IDENTIFICATION: RFOTSELECT
C.PURPOSE: Select objects from a displayed frame and stored them in an
C          intermediate table
C.AUTHOR:  R. Buonanno, G. Buscema, C. Corsi, I. Ferraro, G. Iannicola
C          Osservatorio Astronomico di Roma
C.VERSION: 870417  First running version at ESO (outside MIDAS)   R. Buonanno
C.         881031  New version                                    R. Buonanno
C.         890220  Major changes for portable MIDAS               R.H. Warmels
C.                 load of frame done by LOAD/IMAGE
C.                 cursor interaction with GETCURS
C----
      IMPLICIT         NONE
      INCLUDE          'MID_REL_INCL:RFOTDECL.INC/NOLIST'
C
C     INTEGER          AREA(4)
      INTEGER          CCOUNT,COOFF
      INTEGER          EC, ED, EL
      INTEGER          IWND(2)    
      INTEGER          IVX(9)
      INTEGER*8        IPNTR
      INTEGER          IMF
      INTEGER          I, IT, INX
      INTEGER          IAC, ICOM, IROW, IDUM
      INTEGER          ISTAT,ISTAT1,ISTAT2
      INTEGER          J
      INTEGER          KUN, KNUL, KL
      INTEGER          KDX, KDY
      INTEGER          LX1, LY1
      INTEGER          LX, LY
      INTEGER          MADRID(1)
      INTEGER          NPIX(3)
      INTEGER          NC, NLIV
      INTEGER          NCP,NHL, NCH
      INTEGER          NAXIS, NPL,NL
      INTEGER          XY1,XY2
C
      REAL             AL
      REAL             BB, B1
      REAL             DIMX, DIMY
      REAL             FPIX1(2),FPIX2(2)
      REAL             FG, FO1
      REAL             H1
      REAL             PAR(4)
      REAL             PARS(4)
      REAL             RINPUT(4)
      REAL             RX(100),RY(100)
      REAL             RINTD1, RINTD2, RINTD3, RINTD4
      REAL             RINTD5, RINTD6, RINTD7
      REAL             RME
      REAL             SQR, SQM
      REAL             SL
      REAL             X1, Y1
      REAL             VET(100)
      REAL             VALUE1,VALUE2
      REAL             VME
      REAL             UU, U1
      REAL             VV
      REAL             WC01(2),WC02(2)
C
      DOUBLE PRECISION BEGIN(3),STEP(3),DDUM
C
      LOGICAL          NUL
     
      CHARACTER*60     FRAME
      CHARACTER*72     CUNIT
      CHARACTER*80     IDENT
      CHARACTER*60     INTFIL
      CHARACTER*3      ACTION
      CHARACTER*80     STRING
      CHARACTER*60     CINPUT
C
C *** this part into the common blocks
      INTEGER          NCINT, NRINT, NSINT, NACINT, NARINT
      INTEGER          TIDINT
      INTEGER          NGRP, NOBJ, NINT
C 
      INCLUDE          'MID_INCLUDE:IDIMEM.INC'
      INCLUDE          'MID_INCLUDE:IDIDEV.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/
      DATA             NCP/1/
      DATA             NHL/0/
C
 9001 FORMAT('*** INFO: Intermediate table contains', I5,
     2       ' groups with ',I5,' components')
 9004 FORMAT('          table will be appended')
 9003 FORMAT('Enter identification, mag., col1 and col2 [',I6,
     2       ',0.0,0.0,0.0]: ')
C
C *** start the code
      CALL STSPRO('SELECT')
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,INTFIL,KUN,KNUL,ISTAT)   ! intermediate file
C
C *** open or create the intermediate table
      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 = '*** INFO: Intermediate table not present; '//
     2            'will create a new one ...'
         IDNGRP = 0
         NGRP   = 0
         NOBJ   = 0
         IROW   = 0
         CALL STTPUT(STRING,ISTAT)
         CALL INTINI(INTFIL,TIDINT)

      ELSE
         CALL TBIGET(TIDINT,NCINT,NRINT,NSINT,NACINT,NARINT,ISTAT)
         IF (ISTAT.NE.0) THEN
            STRING = '*** FATAL: Problems with getting info for '//
     2               'intermediate table ...'
           CALL STTPUT(STRING,ISTAT)
           CALL STSEPI
         ENDIF
C
         IF (NRINT.EQ.0) THEN
            STRING = '*** INFO: No points in the intermediate table'
            CALL STTPUT(STRING,ISTAT)
         ENDIF
C
C *** get the intermediate table parameters
         CALL INTDRD(TIDINT,NGRP,NOBJ,NINT,RINTD1,RINTD2,RINTD3,
     2                      RINTD4,RINTD5,RINTD6,RINTD7)
         IROW = NRINT                                   ! number of rows written
         WRITE(STRING,9001) NGRP,NOBJ
         CALL STTPUT(STRING,ISTAT)
         WRITE(STRING,9004)
         CALL STTPUT(STRING,ISTAT)
         CALL TBERDI(TIDINT,NRINT,1,IDNGRP,NUL,ISTAT)
      ENDIF
      CALL STECNT('PUT',EC,ED,EL)
C
C *** open the display device
      CALL DTOPEN(1,ISTAT)
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 *** get the window size
      CALL STKRDI('INPUTI',1,2,IAC,IWND,KUN,KNUL,ISTAT)
      DIMX = FLOAT(IWND(1))
      DIMY = FLOAT(IWND(2))
      IF (DIMX.GT.55) THEN
         DIMX = 55
      ENDIF
      IF (DIMY.GT.55) THEN
         DIMY = 55
      ENDIF
      KDX = DIMX/2+1
      KDY = DIMY/2+1
C
C *** infinite loop for cursor interaction
      COOFF     = 0
      CCOUNT    = 0
      ACTION    = 'YYN'
      FRAME     = ' '
      FLGCMP(1) = 1
      CALL STTPUT('*** INFO: Use cursor control panel to '//
     2                       'select object',ISTAT)
C
 1000 CONTINUE
      IDNGRP = IDNGRP + 1
 1001 CONTINUE
      CALL GETCUR(ACTION,FRAME,XY1,FPIX1,WC01,VALUE1,ISTAT1,
     2                         XY2,FPIX2,WC02,VALUE2,ISTAT2)
      IF (ISTAT1.EQ.0 .AND. ISTAT2. EQ. 0) THEN
         IF (CCOUNT.EQ.0) THEN
            IF (COOFF.EQ.1) THEN
               CALL DTCLOS(QDSPNO)
               CALL STSEPI
            ELSE
      	       FRAME     = ' '
               CALL STTPUT('*** WARNING: Switch cursor on ...'//
     2                     ' next time we exit',ISTAT)
            ENDIF
            COOFF = 1
            GO TO 1001

         ELSE
            NINT   = 0
            RINTD1 = 15000.0
            RINTD2 = 1.
            RINTD3 = 3.
            RINTD4 = 4.
            RINTD5 = 0.0
            RINTD6 = 0.0
            RINTD7 = 0.0
            CALL INTDWR(TIDINT,NGRP,NOBJ,NINT,RINTD1,RINTD2,RINTD3,
     2                         RINTD4,RINTD5,RINTD6,RINTD7)
            CALL TBSINI(TIDINT,ISTAT)
            CALL TBTCLO(TIDINT,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) IDNGRP
         NCH = INDEX(STRING,':')+1               
         CINPUT = ' '
         CALL STKPRC(STRING(1:NCH),'INPUTC',1,1,60,IAC,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
            ICOM=4
            CALL GENCNV(CINPUT,2,ICOM,IDUM,RINPUT,DDUM,ISTAT)
C           CALL USRINP(RINPUT,4,'R',CINPUT)
C           ICOM = NEL(RINPUT,4)
            IF (ICOM.GT.0 .AND. ICOM.LT.4) THEN
               DO 3003 ICOM = ICOM+1, 4
                  RINPUT(ICOM) = 0.0
 3003          CONTINUE
            ENDIF
            IDNGRP= INT(RINPUT(1))
            IF (IDNGRP.GT.99999) THEN
               CALL STTPUT('*** WARNING: Group identification should'//
     2                     ' be less than 100000; try again ...',ISTAT)
               GO TO 3002

            ELSE
               VV     = RINPUT(2)
               BB     = RINPUT(3)
               UU     = RINPUT(4)
            ENDIF 
         ENDIF
C
C *** prepare the arrays for writting
         IDNGRP     = IDNGRP
         IDNCMP(1)  = 101
         IROW       = IROW + 1
         B1         = VV+BB
         U1         = B1+UU
         LX1        = BEGIN(1)+X1-KDX
         LY1        = BEGIN(2)+Y1-KDY
         PARINT(1)  = FLOAT(LX1)
         PARINT(2)  = FLOAT(LY1)
         PARINT(3)  = VV
         PARINT(4)  = B1
         PARINT(5)  = U1
         PARINT(6)  = DIMX
         PARINT(7)  = DIMY
         PARINT(8)  = 0.0
         PARINT(9)  = 0.0
         PARINT(10) = FO1
         PARINT(11) = 4.0
         PARINT(12) = 0.0
         PARINT(13) = 0.0
         PARINT(14) = 0.0
         PARINT(15) = NCP
         PARINT(16) = NHL
C
         FITCMP(1)  = H1
         FITCMP(2)  = FLOAT(KDX)
         FITCMP(3)  = FLOAT(KDY)
         FITCMP(4)  = 3.0
         FITCMP(5)  = 0.0
         FITCMP(6)  = 0.0
C           
         CALL INTWWR(TIDINT,IROW,NCP,NHL)
         NOBJ       = NOBJ + 1
         NGRP       = NGRP + 1
         GO TO 1000
      ENDIF

      CALL STSEPI
      END
