C @(#)necsord.for	19.1 (ESO-DMD) 02/25/03 14:20:25
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 ECHSRD
C
C+
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  23:20 - 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
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(1)              W IS THE ESTIMATED WIDTH OF THE ORDER ACCROSS THE
C                         DISPERSION DIRECTION IN USER COORDINATES
C  ECHR(2)              T THRESHOLD IN THE DETECTION
C  ECHR(4)              S SLOPE OF THE ORDER IN THE MIDDLE OF THE
C                         IMAGE TAG(ALFA)
C  ECHC(17:20)  METH   - STD STANDARD STAR METHOD (DEFAULT)
C      COM COMPLEMENT
C
C  010621		last modif
C 
C-
      IMPLICIT NONE
C 
      INTEGER MADRID
      INTEGER ACTVAL,I,IADD,IDELTA,IPOINT
      INTEGER IST,ISTAT,NACOL,NAUX
      INTEGER NAXIS,NDIM,NORDER,NTOT
      INTEGER NWIN1,NWIN2,PNTRA,STATUS,TID
C 
      INTEGER*8  IXCOL,IYAUX
C
      INTEGER NPIX(3),ICOL(7),IAUX(300),KUN,KNUL,IMNI
      INTEGER NSTART(300),NDATA(300),NPOINT(300),IHW(300)
C 
      REAL RPAR(4),LEVEL,SLOPE
      REAL ORDPOS(300),XAUX(300)
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
C
      CALL STSPRO('ECHSRD')
      CALL STKRDC('ECHC',1,1,8,ACTVAL,FRMIN,KUN,KNUL,STATUS)
      CALL STKRDC('ECHC',1,17,4,ACTVAL,METH,KUN,KNUL,STATUS)
C CALL SXKGET('ECHC', 'C', 21, 8, ACTVAL, TABLE, STATUS)
      TABLE  = 'ORDER'
      CALL STKRDR('ECHR',1,4,ACTVAL,RPAR,KUN,KNUL,STATUS)
      LEVEL  = RPAR(2)
      SLOPE  = RPAR(4)
      IF (METH.EQ.'CO') THEN
          LABCOL(3) = 'YBKG'
          LABCOL(4) = 'Y'
      END IF
      CALL STTPUT(' ECHELLE DEFINITION',ISTAT)
      CALL STTPUT(' ------------------',ISTAT)
      CALL STTPUT(' INPUT IMAGE  : '//FRMIN,ISTAT)
      CALL STTPUT(' OUTPUT TABLE : '//TABLE,ISTAT)
      CALL STTPUT(' PARAMETERS  ',ISTAT)
      WRITE (LINE,9000) RPAR(1)
      CALL STTPUT(LINE,ISTAT)
      WRITE (LINE,9010) RPAR(2)
      CALL STTPUT(LINE,ISTAT)
      WRITE (LINE,9020) SLOPE
      CALL STTPUT(LINE,ISTAT)
      IF (METH.EQ.'CO') THEN
          LINE   = '       METHOD : COMPLEMENT'
      ELSE
          LINE   = '       METHOD : STANDARD'
      END IF
      CALL STTPUT(LINE,ISTAT)
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)
      NWIN2  = RPAR(1)
      NWIN1  = NWIN2/3 + 1
C
C  extract column
C
      CALL TDMGET(4*NPIX(2),IXCOL,IST)
      CALL TDMGET(4*NPIX(2),IYAUX,IST)
      IPOINT = NPIX(1)/2
      CALL EXTR3(MADRID(PNTRA),NPIX(1),NPIX(2),IPOINT,1,NPIX(2),
     +           MADRID(IXCOL))
      IF (METH.EQ.'CO') CALL VCOM(NPIX(2),MADRID(IXCOL),MADRID(IXCOL))
C
C  find order positions in the middle of the image
C
      NAUX   = 200
      CALL FINDM1(MADRID(IXCOL),NPIX(2),NWIN1,NWIN2,LEVEL,NSTART,NDATA,
     +            NORDER,MADRID(IYAUX),XAUX,IAUX,NAUX)
      CALL TDMFRE(4*NPIX(2),IXCOL,IST)
      CALL TDMFRE(4*NPIX(2),IYAUX,IST)
      DO 10 I = 1,NORDER
          ORDPOS(I) = FLOAT(NSTART(I)) + FLOAT(NDATA(I)-1)*0.5
   10 CONTINUE
      WRITE (LINE,9030) NORDER
      CALL STTPUT(LINE,ISTAT)
C
C  get order positions
C
      IDELTA = (NPIX(1)-IPOINT)/10
      IF (NORDER.EQ.1) THEN
          IHW(1) = MIN0(NPIX(2),2*NDATA(1))
      ELSE
          DO 20 I = 1,NORDER - 1
              IHW(I) = (ORDPOS(I+1)-ORDPOS(I))
   20     CONTINUE
          IHW(NORDER) = IHW(NORDER-1)
      END IF
C
C  allocate table
C
      NTOT   = NORDER*30
      CALL TBTINI(TABLE,F_TRANS,F_O_MODE,12,NTOT,TID,STATUS)
      DO 30 I = 1,NACOL
          CALL TBCINI(TID,D_R4_FORMAT,1,FORM(I),UNIT(I),
     .                LABCOL(I),ICOL(I),STATUS)
   30 CONTINUE
      CALL TBCMAP(TID,0,IADD,STATUS)
C
C  follow maxima
C
      CALL FOLLOW(MADRID(PNTRA),NPIX(1),NPIX(2),IPOINT,ORDPOS,NORDER,
     +            IDELTA,IHW,SLOPE,NWIN1,NWIN2,LEVEL,MADRID(IADD),NTOT,
     +            NACOL,NPOINT,METH)
C
C  update number of elements
C
      NTOT   = 0
      DO 40 I = 1,NORDER
          NTOT   = NTOT + NPOINT(I)
   40 CONTINUE
      CALL STKWRI('ECHI',NORDER,4,1,KUN,STATUS)
      CALL STKWRI('ECHI',NPIX,7,2,KUN,STATUS)
      CALL TBIPUT(TID,NACOL,NTOT,STATUS)
      IDENA = 'ORDER POSITION'
      CALL STDWRC(TID,'IDENT',1,IDENA,1,72,KUN,STATUS)
C      CALL DSCUPX(TABLEF,TABLEF,' ',STATUS)
      CALL TBSINI(TID,STATUS)
      CALL TBTCLO(TID,STATUS)
C
C  free data
C
      CALL STSEPI
 9000 FORMAT ('  ORDER WIDTH : ',F6.1)
 9010 FORMAT ('    THRESHOLD : ',F6.1)
 9020 FORMAT ('        SLOPE : ',F6.3)
 9030 FORMAT (' NUMBER OF DETECTED ORDERS : ',I4)
      END

      SUBROUTINE VCOM(N,X,Y)
C
C COMPUTE COMPLEMENT IN VECTOR FORM
C
      IMPLICIT NONE
      INTEGER    N,I
      REAL       X(N),Y(N)

      DO 10 I = 1,N
          Y(I)   = -X(I)
   10 CONTINUE
      RETURN

      END


      SUBROUTINE FOLLOW(X,N1,N2,IXST,ORDPOS,NORDER,IDELTA,IHW,SLOPE,
     +                  NWI1,NWI2,THRES,TAB,NTOT,NACOL,NP,METH)
      IMPLICIT NONE
C
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  SUBROUTINE  FOLLOW      VERSION 1   28 FEB 1983
C  J.D.PONZ
C
C.KEYWORDS:
C  SPECTROSCOPY, CASPEC
C
C.PURPOSE:
C  DEFINE ORDER POSTIONS BY FOLLOWING STRUCTURES ABOVE A GIVEN THRESHOLD
C
C.ALGORITHM:
C
C  ORDER POSITIONS ARE CALCULATED BY THE ROUTINE FINDLT
C
C.INPUT/OUTPUT
C  INPUT ARGUMENTS
C
C  X(N1, N2) REAL*4 INPUT ARRAY
C  N1, N2       INTG*4 DIMENSIONS OF X
C  IXST         INTG*4  STARTING X POSITION IN PIXELS
C  ORDPOS(NORDER)REAL*4 Y POSITION OF THE ORDERS IN PIXELS
C  NORDER INTG*4 NUMBER OF ORDERS
C  IDELTA INTG*4 X INTERVAL IN PIXELS
C  IHW(NORDER) INTG*4 HALFWIDTH OF THE Y-INTERVAL TO BE SEARCHED
C  SLOPE REAL INITIAL SLOPE OF THE ORDERS
C  NWI INTG*4 NUMBER OF PIXELS IN THE 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  NPOS  INTG NUMBER OF LOCATED POSITIONS
C
C-------------------------------------------------------------
C
C
      INTEGER   N1,N2,NLIM,I,IXST,NORDER,IDELTA,NWI1,NWI2
      INTEGER   NTOT,NACOL,NP,ISTAT,IORD,NPF,IDELT1,NL
      INTEGER   NPIX,IX0,IY01,IY0,KZ6907,NPIX1,NMAX,I1,I2
      INTEGER   IB1,IB2,NR,KZ6910
      REAL      SLOPE,THRES,SLOPE1,SLMAX,SLMIN,BACK,AVE,Y,A
      REAL      SLOPE2,F

      REAL X(N1,N2),ORDPOS(NORDER),XPIX(300)
      REAL XL(50),YL(50),XAUX(300),YAUX(300)
      INTEGER IHW(NORDER),IAUX(300),NSTART(10),NDATA(10)
      REAL TAB(NTOT,0:NACOL),TRNULL
      DOUBLE PRECISION       TDNULL
      INTEGER                TINULL
      CHARACTER*2 METH
      CHARACTER*80 LINE
C
      CALL TBMNUL(TINULL,TRNULL,TDNULL)
C
      NLIM   = N1 - 2
      NP     = 0
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
      DO 90 IORD = 1,NORDER
C
C ... INCLUDE START + END POINTS
C
          NPF = 0
          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
C
C ... LEFT PART OF THE IMAGE
C
          IDELT1 = IDELTA
          NL     = 1
          SLOPE1 = SLOPE
          SLMAX  = SLOPE
          SLMIN  = SLOPE
          NPIX   = 2*IHW(IORD)
          IX0    = IXST
          IY01   = ORDPOS(IORD) - IHW(IORD) + 0.5
          IY0    = MAX0(1,IY01)
          XL(1)  = IX0 + IDELT1
          YL(1)  = ORDPOS(IORD) + IDELT1*SLOPE1
          NPIX   = MIN0(NPIX,N2-IY0)
          DO 20 KZ6907 = 1,2147483647
              IF ( .NOT. (IX0.GE.2.AND.NPIX.GT.NWI2)) GO TO 30
C
C ... EXTRACT COLUMN
              NPIX1  = NPIX - (IY0-IY01)
              CALL EXTR3(X,N1,N2,IX0,IY0,NPIX1,XPIX)
              IF (METH.EQ.'CO') CALL VCOM(NPIX1,XPIX,XPIX)
C
C
C ... FIND MAXIMUM
              IF (METH.EQ.'ST' .OR. METH.EQ.'CO') THEN
                  CALL FINDM1(XPIX,NPIX1,NWI1,NWI2,THRES,NSTART,
     +                        NDATA,NMAX,YAUX,XAUX,IAUX,300)

              ELSE
                  CALL FINDM5(XPIX,NPIX1,NWI2,THRES,NSTART,NDATA,
     +                        NMAX)
              END IF

              IF (NMAX.EQ.1) THEN
                  NL     = NL + 1
                  XL(NL) = IX0
                  YL(NL) = FLOAT(NSTART(1)) + FLOAT(NDATA(1)-1)*0.5 +
     +                     IY0 - 1.
                  IDELT1 = IDELTA
C
C ... COMPUTE MOMENTA
                  NP     = NP + 1
                  I1     = NSTART(1)
                  I2     = NSTART(1) + NDATA(1) - 1
                  BACK   = (XPIX(I1-1)+XPIX(I1-2)+XPIX(I1-3)+XPIX(I2+1)+
     +                     XPIX(I2+2)+XPIX(I2+3))/6.
                  AVE    = 0.
                  Y      = 0.
                  DO 10 I = I1,I2
                      F      = XPIX(I) - BACK
                      AVE    = AVE + F
                      Y      = Y + F*I
   10             CONTINUE
                  TAB(NP,1) = IORD
                  TAB(NP,2) = IX0
                  A      = Y/AVE
                  NPF = NPF + 1
                  IF (IORD.NE.1 .OR. IX0.EQ.IXST) THEN
                      TAB(NP,3) = A + IY0 - 1

                  ELSE
                      TAB(NP,3) = TRNULL
                  END IF

                  IB1    = NINT(A-IHW(IORD)/2.)
                  IB2    = NINT(A+IHW(IORD)/2.)
                  IF (IB2.LE.NPIX .AND. IB2.GE.1) THEN
                      IF (METH.EQ.'CO') THEN
                          TAB(NP,5) = -XPIX(NINT(A))

                      ELSE
                          TAB(NP,5) = XPIX(IB2)
                      END IF

                      TAB(NP,4) = IB2 + IY0 - 1

                  ELSE
                      TAB(NP,5) = TRNULL
                      TAB(NP,4) = TRNULL
                  END IF
C       ENDIF
              ELSE
                  NP     = NP + 1
                  TAB(NP,1) = IORD
                  TAB(NP,2) = IX0
                  TAB(NP,3) = TRNULL
                  TAB(NP,4) = TRNULL
                  TAB(NP,5) = TRNULL
                  IDELT1 = IDELT1 + IDELTA
              END IF

              IX0    = XL(NL) - IDELT1
              IF (NL.GT.2) THEN
                  SLOPE2 = (YL(NL)-YL(NL-1))/ (XL(NL)-XL(NL-1))
C       IF (ABS(SLOPE2-SLOPE1).GT.ABS(SLOPE1)*1.0) THEN
C         TYPE *,SLOPE1,SLOPE2
C         NP = NP - 1
C         GO TO 10
C       ENDIF
                  SLOPE1 = SLOPE2
                  SLMAX  = MAX(SLMAX,SLOPE1)
                  SLMIN  = MIN(SLMIN,SLOPE1)
              END IF

              IY01   = YL(NL) - SLOPE1*IDELT1 - IHW(IORD) + 0.5
              IY0    = MAX0(1,IY01)
              NPIX   = 2*IHW(IORD)
              NPIX   = MIN0(NPIX,N2-IY0)
   20     CONTINUE
   30     CONTINUE
C
C ... RIGTH PART OF THE IMAGE
C
          IDELT1 = IDELTA
          NR     = 1
          SLOPE1 = SLOPE
          XL(1)  = IXST
          YL(1)  = ORDPOS(IORD)
          NPIX   = 2*IHW(IORD)
          IX0    = IXST + IDELT1
          IY01   = ORDPOS(IORD) + IDELT1*SLOPE1 - IHW(IORD) + 0.5
          IY0    = MAX0(1,IY01)
          NPIX   = MIN0(NPIX,N2-IY0)
          DO 60 KZ6910 = 1,2147483647
              IF ( .NOT. (IX0.LT.NLIM.AND.NPIX.GT.NWI2)) GO TO 70
C
C ... EXTRACT COLUMN
              NPIX1  = NPIX - (IY0-IY01)
              CALL EXTR3(X,N1,N2,IX0,IY0,NPIX1,XPIX)
              IF (METH.EQ.'CO') CALL VCOM(NPIX1,XPIX,XPIX)
C
C ... FIND MAXIMUM
              IF (METH.EQ.'ST' .OR. METH.EQ.'CO') THEN
                  CALL FINDM1(XPIX,NPIX1,NWI1,NWI2,THRES,NSTART,
     +                        NDATA,NMAX,YAUX,XAUX,IAUX,300)

              ELSE
                  CALL FINDM5(XPIX,NPIX1,NWI2,THRES,NSTART,NDATA,
     +                        NMAX)
              END IF

              IF (NMAX.EQ.1) THEN
                  NR     = NR + 1
                  XL(NR) = IX0
                  YL(NR) = FLOAT(NSTART(1)) + FLOAT(NDATA(1)-1)*0.5 +
     +                     IY0 - 1.
                  IDELT1 = IDELTA
C
C ... COMPUTE MOMENTA
                  NP     = NP + 1
                  I1     = NSTART(1)
                  I2     = NSTART(1) + NDATA(1) - 1
                  BACK   = (XPIX(I1-1)+XPIX(I1-2)+XPIX(I1-3)+XPIX(I2+1)+
     +                     XPIX(I2+2)+XPIX(I2+3))/6.
                  AVE    = 0.
                  Y      = 0.
                  DO 50 I = I1,I2
                      F      = XPIX(I) - BACK
                      AVE    = AVE + F
                      Y      = Y + F*I
   50             CONTINUE
                  TAB(NP,1) = IORD
                  TAB(NP,2) = IX0
                  A      = Y/AVE
                  NPF = NPF + 1
                  IF (IORD.NE.NORDER) THEN
                      TAB(NP,3) = A + IY0 - 1

                  ELSE
                      TAB(NP,3) = TRNULL
                  END IF

                  IB1    = NINT(A-IHW(IORD)/2.)
                  IB2    = NINT(A+IHW(IORD)/2.)
                  IF (IB2.LE.NPIX .AND. IB2.GE.1) THEN
                      IF (METH.EQ.'CO') THEN
                          TAB(NP,5) = -XPIX(NINT(A))

                      ELSE
                          TAB(NP,5) = XPIX(IB2)
                      END IF

                      TAB(NP,4) = IB2 + IY0 - 1

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

              ELSE
                  NP     = NP + 1
                  TAB(NP,1) = IORD
                  TAB(NP,2) = IX0
                  TAB(NP,3) = TRNULL
                  TAB(NP,4) = TRNULL
                  TAB(NP,5) = TRNULL
                  IDELT1 = IDELT1 + IDELTA
              END IF

              IX0    = XL(NR) + IDELT1
              IF (NR.GT.2) THEN
                  SLOPE2 = (YL(NR)-YL(NR-1))/ (XL(NR)-XL(NR-1))
C       IF (ABS(SLOPE2-SLOPE1).GT.ABS(SLOPE1)*1.0) THEN
C         TYPE *,SLOPE1,SLOPE2
C         NP = NP - 1
C         GO TO 20
C       ENDIF
                  SLOPE1 = SLOPE2
                  SLMAX  = MAX(SLMAX,SLOPE1)
                  SLMIN  = MIN(SLMIN,SLOPE1)
              END IF

              IY01   = YL(NR) + SLOPE1*IDELT1 - IHW(IORD) + 0.5
              IY0    = MAX0(1,IY01)
              NPIX   = 2*IHW(IORD)
              NPIX   = MIN0(NPIX,N2-IY0)
   60     CONTINUE
   70     CONTINUE
          WRITE (LINE,9010) IORD,IXST,NINT(ORDPOS(IORD)),IHW(IORD),
     +      NPF,SLMAX,SLMIN
          CALL STTPUT(LINE,ISTAT)
   90 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,2X,F8.2,2X,F8.2)
 9020 FORMAT (' SEQ.NO.  X CENTER  Y CENTER  INTERORDER  TEMPLA',' MAX',
     +       '.SLOPE MIN.SLOPE')
 9030 FORMAT (' -------  --------  --------  ----------  --------',' -',
     +       '-------- ---------')
 9040 FORMAT (' -------------------------------------------------','--',
     +       '------------------')
      END



      SUBROUTINE EXTR3(X,N1,N2,IXS,IYS,N,XCOL)
C+
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  23:45 - 3 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C.KEYWORDS
C
C
C.PURPOSE
C
C  EXTRACTS A COLUMN IN A 2D ARRAY AS AVERAGE OF 3 CONTIGUOUS COLUMNS
C
C.ALGORITHM
C
C  OUTPUT COLUMN IS THE AVERAGE OF THREEE CONTIGUOUS COLUMNS
C
C.INPUT/OUTPUT
C
C  INPUT ARGUMENTS ARE
C  X(N1,N2) REAL*4 INPUT 2D FRAME
C  N1,N2 INTG*4 X-DIMENSION
C  IXS  INTG*4 CENTRAL COLUMN NUMBER
C  IYS  INTG*4 STARTING ROW NUMBER
C  N  INTG*4 NUMBER OF PIXELS IN THE COLUMN
C
C  OUTPUT ARGUMENTS ARE
C  XCOL(N) REAL*4 EXTRACTED COLUMN
C
C-
C
C  EXTRACT COLUMN AT IPOINT
C
C      IMPLICIT NONE
      INTEGER N1,N2,N,I1,IXS,IYS,I,J,I2,J1
      REAL X(N1,N2),XCOL(N)

      I1     = IXS - 1
      I2     = IXS + 1
      I      = IXS
      DO 10 J = 1,N
          J1     = J + IYS - 1
          XCOL(J) = (X(I1,J1)+X(I,J1)+X(I2,J1))/3.0
   10 CONTINUE
      RETURN
      END


C
      SUBROUTINE FINDM1(X,N,NW1,NW2,THRES,IST,NDAT,NPOS,YAUX,XAUX,IAUX,
     +                  NAUX)
C+
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  23:50 - 3 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C.KEYWORDS
C
C  CASPEC
C
C.PURPOSE
C
C  FIND VALUES ABOVE A GIVEN THRESHOLD
C
C.ALGORITHM
C
C
C.INPUT/OUTPUT
C
C  INPUT ARGUMENTS
C  X(N)  REAL INPUT ARRAY
C  N  INTG ARRAY DIMENSION
C  NW1          INTG    SMALL WINDOW
C  NW2          INTG    LARGE WINDOW
C  THRES REAL THRESHOLD
C
C  OUTPUT ARGUMENTS
C  IST  INTG STARTING POSITION OF THE REGION
C  NDAT         INTG    NUMBER OF PIXELS IN THE REGION
C  NPOS  NPOS NUMBER OF REGIONS ABOVE THRESHOLD
C
C  AUXILIARY SPACE
C  YAUX(N)      REAL
C  XAUX(NAUX)   REAL
C  IAUX(NAUX)   INTG
C
C-
C
C      IMPLICIT NONE
      REAL X(N),YAUX(N),XAUX(NAUX),THRES
      INTEGER IST(1),NDAT(1),IAUX(NAUX),N,I,NFIRST
      INTEGER NW1,NW2,NPIX,NPOS,NAUX,MWIND,NO
      LOGICAL FIRST
C
C ... median filtering using threshold
C
      MWIND = 2*NW2 + 1
      CALL MEDNTH(X,N,MWIND,THRES,YAUX,XAUX,IAUX,NAUX)
C
C ... original - filtered
C
      FIRST  = .TRUE.
      NO     = 0
      DO 10 I = NW2,N - NW2 + 1
          YAUX(I) = X(I) - YAUX(I)
          IF (YAUX(I).GT.0.) THEN
C ... above the threshold
              IF (FIRST) THEN
C ... first value above the threshold
                  FIRST  = .FALSE.
                  NFIRST = I
                  NPIX   = 1

              ELSE
                  NPIX   = NPIX + 1
              END IF

          ELSE
C ... below the threshold
              IF ( .NOT. FIRST) THEN
C ... first value below the threshold
                  FIRST  = .TRUE.
                  IF (NPIX.GE.NW1) THEN
                      NO     = NO + 1
                      IST(NO) = NFIRST
                      NDAT(NO) = NPIX
                  END IF

              END IF

          END IF

   10 CONTINUE
      NPOS   = NO
      RETURN

      END


      SUBROUTINE FINDM5(X,N,NW,T,IST,NDAT,NACT)
C+
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  23:56 - 3 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C  SUBROUTINE FINDM5  VERSION 1.0 11 MAR 1983
C  J.D.PONZ   ESO - GARCHING
C
C.KEYWORDS
C
C  FIND SPECTRAL FEATURES, SPECTROSCOPY
C
C.PURPOSE
C
C  FIND FLAT FIELD PROFILES
C
C.ALGORITHM
C
C  SIMPLE THRESHOLDING METHOD
C
C.INPUT/OUTPUT
C
C  INPUT ARGUMENTS
C  X(N)  REAL INPUT ARRAY
C  N  INTG ARRAY DIMENSION
C  NW           INTG    WINDOW SIZE IN PIXELS
C  T  REAL THRESHOLD
C
C  OUTPUT ARGUMENTS
C  IST  INTG POINTER OF THE STARTING POSITION
C  NDAT  INTG NUMBER OF PIXELS ABOVE THE THRESHOLD
C  NACT         INTG    NUMBER OF DETECTED FEATURES
C
C-
C      IMPLICIT NONE
      REAL X(N),T
      INTEGER IST(1),NDAT(1),N,NW,NACT
      INTEGER I,NPIX,NO
      LOGICAL BACK
C
      BACK   = .FALSE.
      NO     = 0
      DO 10 I = 1,N
          IF (X(I).GT.T) THEN
              IF (BACK) THEN
                  BACK   = .FALSE.
                  NO     = NO + 1
                  IST(NO) = I
                  NPIX   = 1

              ELSE
                  NPIX   = NPIX + 1
              END IF

          ELSE
              IF ( .NOT. BACK) THEN
                  BACK   = .TRUE.
                  IF (NO.GT.0) NDAT(NO) = NPIX
              END IF

          END IF

   10 CONTINUE
      NACT   = NO
      IF (NO.EQ.2) THEN
          NACT   = 1

      ELSE
          IF (NO.EQ.3) THEN
              NACT   = 1
              NDAT(1) = NDAT(2)
              IST(1) = IST(2)
          END IF

      END IF

      RETURN

      END


      SUBROUTINE MEDNTH(IA,LEN,WIDTH,THRES,OA,FA,LINK,NAUX)
C+
C
C.IDENTIFICATION
C  SUBROUTINE MEDIAN              VERSION 1.0     8 MAR 1983
C
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  0:34 - 4 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C.KEYWORDS
C
C  FILTERING, MEDIAN FILTERING, SPECTROSCOPY, THRESHOLDING
C
C.PURPOSE
C
C  PERFORM MEDIAN FILTERING AND THRESHOLDING ON THE INPUT ARRAY
C
C.ALGORITHM
C
C  DATA IN A RUNNING WINDOW IS SORTED TO FORM A LOCAL HISTOGRAM
C  THE VALUE IN THE MIDDLE OF THE LOCAL HISTOGRAM IS TAKEN AS
C  OUTPUT VALUE IF ABS(MEDIAN-VALUE) .GT. ABS(THRES)
C
C.INPUT/OUTPUT
C
C  INPUT ARGUMENTS :
C  IA      REAL*4     INPUT ARRAY
C  LEN     INTG*4     NUMBER OF INPUT PIXELS
C  WIDTH   INTG*4     WINDOW WIDTH
C  THRES   REAL*4     THRESHOLD
C
C  OUTPUT ARGUMENTS :
C  OA      REAL*4     OUTPUT ARRAY
C
C  WORK SPACE :
C  FA(NAUX)   REAL*4     AUXILIARY SPACE
C  LINK(NAUX) INTG*4     AUXILIARY SPACE
C  NAUX       INTG*4     SPACE SIZE. NAUX .GT. WIDTH+1
C
C-
C      IMPLICIT NONE
      REAL FA(NAUX),IA(LEN),OA(LEN),THRES,T
      INTEGER LINK(NAUX),MIN,MAX,WIDTH,LEN,MIDPT,FRONT,AVAIL
      INTEGER RANGE,NEWNOD,I,J,K,NODE,ITMP,N,NAUX
C
C  SET UP BEGINNING (MIN) AND END (MAX) OF FILTER ARRAY
C   ALSO SET UP POINTER TO AVAIL LIST AND LIST HEADER (FA(1),LINK(1)
C
      T      = ABS(THRES)
      MIN    = 2
      MAX    = NAUX
      AVAIL  = MIN
      FA(1)  = -1.0
      LINK(1) = 0
C
C  SET UP FIRST FILTER WINDOW
C
      RANGE  = WIDTH - 1
      DO 10 I = 1,RANGE,1
          NODE   = NEWNOD(IA(I),AVAIL,MIN,MAX,FA,LINK,NAUX)
          CALL INSNOD(NODE,FA,LINK,NAUX)
   10 CONTINUE
C
C COPY BORDERS TO OUTPUT BUFFER (CANNOT PROCESS BORDER POINTS)
C
      RANGE  = WIDTH/2
      K      = LEN - RANGE
      DO 20 J = 1,RANGE,1
          OA(J)  = IA(J)
          K      = K + 1
C          K POINTS TO HIGH END OF ARRAY, J POINTS TO LOW
          OA(K)  = IA(K)
   20 CONTINUE
C
C  FOR EACH POINT FROM WIDTH/2 TO LEN-WIDTH/2+1
C   FIND ITS MEDIAN
C
      FRONT  = MIN
C                  FRONT IS THE OLDEST POINT IN THE FILTER ARRAY
      K      = RANGE
C                  K IS THE POINT BEING PROCESSED (LOC IN OA)
      MIDPT  = RANGE + 1
C                  MIDPT IS THE MIDDLE OF THE WINDOW (IN FA)
      DO 40 J = WIDTH,LEN,1
          NODE   = NEWNOD(IA(J),AVAIL,MIN,MAX,FA,LINK,NAUX)
          CALL INSNOD(NODE,FA,LINK,NAUX)
          K      = K + 1
C    FIND THE MIDPT - - IT IS THE MEDIAN VALUE
          ITMP   = LINK(1)
          DO 30 I = 2,MIDPT,1
              ITMP   = LINK(ITMP)
   30     CONTINUE
C        GET THE LINK WHICH POINTS TO THE VALUE OF THE MEDIAN
C        MODIFY ONLY IF GREATHER THAN THRESHOLD
          IF (ABS(FA(ITMP)-IA(K)).GT.T) THEN
              OA(K)  = FA(ITMP)

          ELSE
              OA(K)  = IA(K)
          END IF

          CALL DELNOD(FRONT,FA,LINK,NAUX)
          FRONT  = FRONT + 1
          IF (FRONT.GT.MAX) FRONT  = MIN
C                    DELETE OLDEST AND ADVANCE PTR, FRONT ALSO WRAPS
   40 CONTINUE
C
C USE SMOOTH VALUE AT THE ENDS
C
      N      = WIDTH/2 + 1
      DO 50 I = 1,N
          OA(I)  = OA(N)
          OA(LEN-I+1) = OA(LEN-N)
   50 CONTINUE
      RETURN

      END
