C @(#)necford.for	19.1 (ESO-DMD) 02/25/03 14:20:23
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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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 ECHFOR
C
C+
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  21:22 - 3 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C.IDENTIFICATION
C
C  program  ECHSORD  version 1.6 830310
C
C.MODIFICATIONS
C
C   910128   P. Ballester   Define and set variables for mmake -i
C
C.KEYWORDS
C
C  ECHELLE, CASPEC
C
C.PURPOSE
C
C SEARCH FOR ORDERS IN AN IMAGE. TYPE ECHELLE - FLAT FIELD
C OUTPUT IS A TABLE WITH THE FOLLOWING COLUMNS :
C COLUMN NO.  LABEL   UNIT      DESCRIPTION
C    1        :ORDER  -         ORDER NUMBER
C    2        :X      PIXEL     SAMPLE NUMBER
C    3        :Y      PIXEL     1. MOMENT AS SUM(F*Y)/SUM(F)
C    4        :YBKG   PIXEL     POSITION OF THE BACKGROUND
C    5        :BKG    DN        BACKGROUND LEVEL
C
C.ALGORITHM
C
C  MAXIMA FOLLOWING ALGORITHM IS USED
C
C  CONSTRAINS : MAXIMUM NUMBER OF ORDERS IN THE IMAGE IS 100
C               FOR EACH ORDER 19 POINTS ARE FOUND
C
C.INPUT/OUTPUT
C
C  COMMAND :
C  SEARCH/ORDER INPUT W,T,S OUTPUT METH
C  input keywords are
C  ECHC(1:8)   INPUT   - IMAGE NAME
C  ECHC(21:28)  OUTPUT  - TABLE NAME
C         W,T,S   - PARAMETERS WHERE :
C  ECHR(2)              T THRESHOLD IN THE DETECTION
C
C 010226		last modif
C 
C----------------------------------------------------------

      IMPLICIT  NONE
      INTEGER   MADRID
      INTEGER   ACTVAL,I
      INTEGER   IXSTR
      INTEGER   NACOL,NAXIS,NDIM,NORDER,NP,NPOINT
      INTEGER   NTOT
      INTEGER   STATUS,WINDOW,IACOL,ID
      INTEGER   NPIX(3),ICOL(7),KUN,KNUL,IMNI
      INTEGER   UPPER(100),LOWER(100),ORDPOS(100)
      INTEGER*8 IADD
      INTEGER*8 PNTRA
C 
      REAL    RPAR(4),LEVEL
C 
      DOUBLE PRECISION    START(3),STEP(3)
C 
      CHARACTER*8 FRMIN,TABLE
      CHARACTER*72 IDENA
      CHARACTER*80 LINE
      CHARACTER*64 CUNITA
      CHARACTER*2 METH
      CHARACTER*16 LABCOL(5),UNIT(5)
      CHARACTER*6 FORM(5)
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C 
      DATA LABCOL(1)/'ORDER'/,FORM(1)/'I6'/
      DATA LABCOL(2)/'X'/,FORM(2)/'F7.1'/
      DATA LABCOL(3)/'Y'/,FORM(3)/'F7.1'/
      DATA LABCOL(4)/'YBKG'/,FORM(4)/'F7.1'/
      DATA LABCOL(5)/'BKG'/,FORM(5)/'E12.3'/
      DATA UNIT(1)/'UNITLESS'/
      DATA UNIT(2)/'PIXEL'/
      DATA UNIT(3)/'PIXEL'/
      DATA UNIT(4)/'PIXEL'/
      DATA UNIT(5)/'DN'/
      DATA NACOL/5/
C
      NDIM   = 3
C
C ... initialize system and read parameters
C
      CALL STSPRO('ECHFOR')
      CALL STKRDC('ECHC',1,1,8,ACTVAL,FRMIN,KUN,KNUL,STATUS)
      CALL STKRDC('ECHC',1,17,4,ACTVAL,METH,KUN,KNUL,STATUS)
      TABLE  = 'ORDER'
      CALL STKRDR('ECHR',1,4,ACTVAL,RPAR,KUN,KNUL,STATUS)
      LEVEL  = RPAR(2)
      IF (METH.EQ.'CO') THEN
          LABCOL(3) = 'YBKG'
          LABCOL(4) = 'Y'
      END IF
      CALL STTPUT(' ECHELLE DEFINITION',STATUS)
      CALL STTPUT(' ------------------',STATUS)
      CALL STTPUT(' INPUT IMAGE  : '//FRMIN,STATUS)
      CALL STTPUT(' OUTPUT TABLE : '//TABLE,STATUS)
      CALL STTPUT(' PARAMETERS  ',STATUS)
      WRITE (LINE,9010) RPAR(2)
      CALL STTPUT(LINE,STATUS)
C
C ... read frame
C
      CALL STIGET(FRMIN,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .            NDIM,NAXIS,NPIX,START,STEP,IDENA,
     .            CUNITA,PNTRA,IMNI,STATUS)
C
C ... find order positions in the middle of the image
C
      NORDER = 0
      IXSTR  = NPIX(1)/2
      WINDOW = RPAR(1)
      NP     = NPIX(2) - 2*WINDOW
      CALL FINDM1(MADRID(PNTRA),NPIX(1),NPIX(2),IXSTR,WINDOW,NP,LEVEL,
     +            NORDER,ORDPOS,UPPER,LOWER)
      IF (NORDER.EQ.0) THEN
          CALL STTPUT('Order detection failed',STATUS)
          CALL STTPUT('Use a lower threshold',STATUS)
          GO TO 20
      END IF
      WRITE (LINE,9030) NORDER
      CALL STTPUT(LINE,STATUS)
C
C ... allocate table
C
      NTOT   = (NORDER+1)*NPIX(1)/10
      IACOL  = 12
      CALL TBTINI(TABLE,F_TRANS,F_O_MODE,IACOL,NTOT,ID,STATUS)
      DO 10 I = 1,NACOL
          CALL TBCINI(ID,D_R4_FORMAT,1,FORM(I),
     .                UNIT(I),LABCOL(I),ICOL(I),
     +                STATUS)
   10 CONTINUE
      CALL TBCMAP(ID,0,IADD,STATUS)
C
C ... follow orders
C
      CALL FOLLOW(MADRID(PNTRA),NPIX(1),NPIX(2),NORDER,ORDPOS,UPPER,
     +            LOWER,WINDOW,LEVEL,MADRID(IADD),NTOT,NACOL,NPOINT)
C
C ... update number of elements
C
      NTOT   = NPOINT
      CALL STKWRI('ECHI',NORDER,4,1,KUN,STATUS)
      CALL STKWRI('ECHI',NPIX,7,2,KUN,STATUS)
      CALL TBIPUT(ID,NACOL,NTOT,STATUS)
      IDENA = 'ORDER POSITION'
C      CALL SXDPUT(TABLEFULL,'IDENT','C',IDENA,1,72,STATUS)
C      CALL DSCUPD(TABLEFULL,TABLEFULL,' ',STATUS)
      CALL TBSINI(ID,STATUS)
      CALL TBTCLO(ID,STATUS)
C
C ... free data
C
   20 CALL STSEPI
C 9000 FORMAT ('  ORDER WIDTH : ',F6.1)
 9010 FORMAT ('    THRESHOLD : ',F6.1)
C 9020 FORMAT ('        SLOPE : ',F6.3)
 9030 FORMAT (' NUMBER OF DETECTED ORDERS : ',I4)
      END

      SUBROUTINE FOLLOW(X,N1,N2,NORDER,ORDPOS,UPPER,LOWER,WINDOW,THRES,
     +                  TAB,NTOT,NACOL,NP)
C
C  DEFINE ORDER POSTIONS BY FOLLOWING STRUCTURES ABOVE A GIVEN THRESHOLD
C      IMPLICIT NONE
C  INPUT ARGUMENTS
C
C  X(N1, N2) REAL*4 INPUT ARRAY
C  N1, N2       INTG*4 DIMENSIONS OF X
C  NORDER INTG*4 NUMBER OF ORDERS
C  ORDPOS INTG*4 Y POSITION OF THE ORDERS IN PIXELS
C  LOWER,UPPER  INTG*4 LOWER AND UPPER ORDER POSITIONS
C  WINDOW INTG*4 NUMBER OF PIXELS IN THE EXCLUDING WINDOW
C  THRES REAL THRESHOLD
C  NTOT  INTG*4 ALLOCATED SPACE
C
C  OUTPUT ARGUMENTS
C
C  TAB   REAL ARRAY WITH THE LOCATED POSITIONS IN FORMAT
C                   (NTOT,0:NACOL) 1 - ORDER NUMBER
C               2 - SAMPLE
C       3 - MOMENT ORDER 1
C                            4 - position of back.
C       5 - BACKGROUND
C  NP  INTG NUMBER OF LOCATED POSITIONS
C
C-------------------------------------------------------------
C
C
      IMPLICIT NONE

      INTEGER       N1,N2,NORDER,NTOT,NP,NACOL
      INTEGER       NSAMP,ISTAT,IORD,IRL,IL
      REAL          THRES
      REAL          X(N1,N2)
      INTEGER       ORDPOS(NORDER),LOWER(NORDER),UPPER(NORDER)
      INTEGER       ICPOS(20),ILPOS(20),IUPOS(20),WINDOW
      INTEGER       TINULL
      DOUBLE PRECISION TDNULL
      REAL          TAB(NTOT,0:NACOL),TRNULL
      CHARACTER*80  LINE
      LOGICAL       NEXT
      INTEGER       NF,IHW,IXSTR,IBACK,ISTEP,IX,IY,I1,IU
      INTEGER       KZ6902,NPTS,NFOUND
C
      CALL TBMNUL(TINULL,TRNULL,TDNULL)
      NSAMP  = 10
      IXSTR = N1/2
C
C ... DISPLAY RESULTS
C
      CALL STTPUT(' ',ISTAT)
      WRITE (LINE,9020)
      CALL STTPUT(LINE,ISTAT)
      WRITE (LINE,9030)
      CALL STTPUT(LINE,ISTAT)
C
C ... ITERATION ON ORDERS
C
      NP     = 0
      DO 40 IORD = 1,NORDER
C
C ... INCLUDE START + END POINTS
C
          NF     = 0
          IF (IORD.NE.NORDER) IHW    = (ORDPOS(IORD+1)-ORDPOS(IORD))/2
          NP     = NP + 1
          TAB(NP,1) = IORD
          TAB(NP,2) = 1.
          TAB(NP,3) = TRNULL
          TAB(NP,4) = TRNULL
          TAB(NP,5) = TRNULL
          NP     = NP + 1
          TAB(NP,1) = IORD
          TAB(NP,2) = N1
          TAB(NP,3) = TRNULL
          TAB(NP,4) = TRNULL
          TAB(NP,5) = TRNULL
          NP     = NP + 1
          TAB(NP,1) = IORD
          TAB(NP,2) = IXSTR
          TAB(NP,3) = ORDPOS(IORD)
          IBACK  = ORDPOS(IORD) + IHW
          IF (IBACK.LT.N2-WINDOW) THEN
              TAB(NP,4) = IBACK
              TAB(NP,5) = X(IXSTR,IBACK)

          ELSE
              TAB(NP,4) = TRNULL
              TAB(NP,5) = TRNULL
          END IF
C
C ... LEFT/RIGHT  PART OF THE IMAGE
C
          DO 30 IRL = 1,2
              IF (IRL.EQ.1) THEN
                  ISTEP  = 1

              ELSE
                  ISTEP  = -1
              END IF

              IX     = IXSTR
              IL     = LOWER(IORD)
              IU     = UPPER(IORD)
              NEXT   = .TRUE.
99            KZ6902 = 1
                  IF ( .NOT. (NEXT)) GO TO 20
                  IX     = IX + ISTEP
                  IY     = MAX(IL-3,WINDOW)
                  I1     = MIN(N2-WINDOW,IU+3)
                  NPTS   = I1 - IY + 1
                  CALL FINDM1(X,N1,N2,IX,IY,NPTS,THRES,NFOUND,ICPOS,
     +                        IUPOS,ILPOS)
C
C ... CHECK IF ONE ORDER FOUND
C
                  IF (NFOUND.EQ.1) THEN
C
C ... STORE INTO TABLE
C
                      IF (MOD(IX,NSAMP).EQ.0) THEN
                          NF     = NF + 1
                          NP     = NP + 1
                          TAB(NP,1) = IORD
                          TAB(NP,2) = IX
                          TAB(NP,3) = ICPOS(1)
                          IBACK  = ICPOS(1) + IHW
                          IF (IBACK.LT.N2-WINDOW) THEN
                              TAB(NP,4) = IBACK
                              TAB(NP,5) = X(IXSTR,IBACK)

                          ELSE
                              TAB(NP,4) = TRNULL
                              TAB(NP,5) = TRNULL
                          END IF

                      END IF

                      IL     = ILPOS(1)
                      IU     = IUPOS(1)
                      NEXT   = IL .GT. WINDOW .AND. IU .LT.
     +                         (N2-WINDOW) .AND. IX .GT. WINDOW .AND.
     +                         IX .LT. (N1-WINDOW)

                  ELSE
                      NEXT   = .FALSE.
                  END IF

              GOTO 99
   20         CONTINUE
   30     CONTINUE
          WRITE (LINE,9010) IORD,IXSTR,ORDPOS(IORD),IHW,NF
          CALL STTPUT(LINE,ISTAT)
   40 CONTINUE
      WRITE (LINE,9040)
      CALL STTPUT(LINE,ISTAT)
      RETURN

C 9000 FORMAT (I4)
 9010 FORMAT (1X,I7,2X,I8,2X,I8,2X,I10,2X,I8)
 9020 FORMAT (' SEQ.NO.  X CENTER  Y CENTER  INTERORDER  TEMPLA')
 9030 FORMAT (' -------  --------  --------  ----------  --------')
 9040 FORMAT (' -------------------------------------------------')
      END

      SUBROUTINE FINDM1(X,N1,N2,IX,IY,NP,THRES,NORDER,ORDPOS,UPPER,
     +                  LOWER)
C
      IMPLICIT NONE
C 
C  FIND CENTER OF THE ORDERS IN THE MIDDLE OF THE IMAGE
C
C  INPUT ARGUMENTS
C  X(N1,N2) REAL INPUT ARRAY
C  N1,N2 INTG ARRAY DIMENSION
C  IX, IY INTG    STARTING PIXEL POSITIONS
C  NP  INTG    NUMBER OF PIXELS TO SEARCH
C  THRES REAL THRESHOLD
C
C  OUTPUT ARGUMENTS
C  NORDER INTG NUMBER OD ORDERS
C  ORDPOS INTG    POSITION OF THE ORDER (UPPER-LOWER)/2
C  UPPER INTG UPPER ORDER LIMIT
C  LOWER INTG LOWER ORDER LIMIT
C
      INTEGER    N1,N2,IX,IY,NP,NORDER,ISTART,NO,I,NFIRST,NPIX
      INTEGER ORDPOS(1),UPPER(1),LOWER(1),WIDTH,WIDTH1,STATUS
C 
      REAL       THRES
      REAL X(N1,N2)
C 
      LOGICAL FIRST
C
C 
C 
      ISTART = IX
      WIDTH  = 0
C
C ... FIND UPPER LEVEL
C
      FIRST  = .TRUE.
      NO     = 0
      DO 10, I = IY,IY + NP - 1
          IF (X(ISTART,I).GT.THRES) THEN
C
C ... ABOVE THE THRESHOLD
C
              IF (FIRST) THEN
C
C ... FIRST VALUE ABOVE THE THRESHOLD
C
                  FIRST  = .FALSE.
                  NFIRST = I
                  NPIX   = 1
              END IF

          ELSE
C
C ... BELOW THE THRESHOLD
C
              IF ( .NOT. FIRST) THEN
C
C ... FIRST VALUE BELOW THE THRESHOLD
C
                  FIRST  = .TRUE.
                  NO     = NO + 1
                  UPPER(NO) = I - 1
                  LOWER(NO) = NFIRST
                  ORDPOS(NO) = LOWER(NO) + (UPPER(NO)-LOWER(NO))/2
                  IF (WIDTH.EQ.0) THEN
                      WIDTH  = UPPER(NO) - LOWER(NO)

                  ELSE
                      WIDTH1 = UPPER(NO) - LOWER(NO)
                      IF (ABS(WIDTH1-WIDTH).GT.0.1*WIDTH) 
     +                   CALL STTPUT(
     +                  'Warning: Order width changes',STATUS)
                      WIDTH  = WIDTH1
                  END IF

              END IF

          END IF

   10 CONTINUE
      NORDER = NO
      RETURN

      END
