C @(#)neczebpar.for	19.1 (ES0-DMD) 02/25/03 14:20:26
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
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1991 Special Astrophisical Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, 26 Aug 1991
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: A.Yu.Kniazev, V.Shergin
C
C.IDENTIFICATION
C                  ZEB_PAR.FOR
C
C.KEYWORDS
C  echelle (ZEBRA), parameters
C
C.PURPOSE
C Determination of the parameters for echelle images (ZEBRA)
C in single observation set.
C
C
C.ALGORITHM
C
C.INPUT/OUTPUT
C  input keywords are:
C  ZEBRAC(1:8)   INPUT   - IMAGE NAME
C  ZEBRAC(20)    INPUT   - PLOT KEY
C  ZEBRAI(1)     INPUT   - the number of column for begin work
C  ZEBRAI(2)     INPUT   - the half of zone for begin work
C
C Program output parameters in keywords:
C   ZEBRAR(1)  -    the mean slope of echelle spectra
C   ZEBRAR(2)  -    the mean width of orders
C
C  COMMAND :
C  DETERMIN/ZEB INPUT
C
C-----------------------------------------------------------------
      PROGRAM ZEBPAR

      IMPLICIT NONE
      INTEGER   ISLOST, LWORD
      REAL      ACCUR
      PARAMETER (ISLOST = 12)
      PARAMETER (ACCUR  = 0.005)
      INTEGER MADRID(1)
      INTEGER ACTVAL,I,K,ISTAT
      INTEGER NAXIS,NDIM
      INTEGER IVEC1, IVEC2
      INTEGER PNTRA,STATUS,WINDOW
      INTEGER NPIX(3), KUN, KNUL, IMNO, IJ
      DOUBLE PRECISION START(3),STEP(3)
      DOUBLE PRECISION ST
      CHARACTER*80 CA
      CHARACTER*72 IDENA
      CHARACTER*80 LINE
      CHARACTER*64 CUNITA
      CHARACTER*8  FRMIN
      CHARACTER*1  PLOT
C
      REAL INCL, XINCL  !  orders inclination coefficient
      REAL SINCL        !  step of change INCL
      INTEGER ORDWID    !  average width of ORDER
      INTEGER COLUMN(2) !  number of column for begin work
      INTEGER IDIR, FLAG
      REAL    AM        !  amplitude
      REAL    XC        !
      REAL    HWHI      !  average width of ORDER (It's main)
C      REAL    CUTS(2)   !  minimum, maximum
C      INTEGER PIX(2)    !  index of minimum and maximum pixel
      REAL    TMP(-1:1) !  output value width of order
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      CA = 'COUNTS'
      ST = 1.D0
      NDIM   = 3
C
C ... initialize system and read parameters
C
      CALL STSPRO('ZEBPAR')
      CALL STKRDC('IN_A',1,1,60,ACTVAL,FRMIN,KUN,KNUL,STATUS)
      CALL STKRDI('INPUTI',1,2,ACTVAL,COLUMN,KUN,KNUL,STATUS)
      CALL STKRDC('INPUTC',1,1,1,ACTVAL,PLOT,KUN,KNUL,STATUS)
C
      CALL STTPUT(' ZEBRA PARAMETERS DETERMINATION',ISTAT)
      CALL STTPUT(' ------------------------------',ISTAT)
      CALL STTPUT(' INPUT IMAGE  : '//FRMIN,ISTAT)
      WRITE (LINE,9040) COLUMN(1)
      CALL STTPUT(LINE ,ISTAT)
      WRITE (LINE,9050) COLUMN(2)
      CALL STTPUT(LINE ,ISTAT)
      CALL STTPUT(' PLOT FLAG  : '//PLOT,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,IMNO,STATUS)
C
C ... generate working space
C
      CALL STIPUT('ZEBRA.TMP',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,1,
     .            NPIX(2)*4,ST,ST,CA,CA,IVEC1,IJ,STATUS)
      IVEC2 = IVEC1 + NPIX(2)
C
C ... start step-by-step fitting process to find optimal inclination
C ... coefficient corresponding to minimum order width
      WINDOW  = NPIX(2)/5
      FLAG    = 0
      IDIR    = 0
      INCL    = 0.
      XINCL   = 0.
      SINCL   = 0.1
C
C ... next step  (or 3 steps when IDIR = 0)
 10   DO I = -1, 1, 1
        IF (IDIR .EQ. I .OR. IDIR .EQ. 0) THEN
        IF (I .NE. 0 .OR. FLAG .EQ. 0) THEN
           INCL = XINCL + SINCL*I
C ... extract and plot "averaging" column
C23456789012345678901234567890123456789012345678901234567890123456789012345678
           CALL EXTRAV(MADRID(PNTRA),NPIX(1),NPIX(2),COLUMN(1),
     .                 COLUMN(2),INCL,MADRID(IVEC1))
         IF (PLOT(1:1) .EQ. 'Y' .OR. PLOT(1:1) .EQ. 'y')
     .    CALL PLTARR(MADRID(IVEC1),NPIX(2),100,(NPIX(2)-150),FRMIN,0)
C ...   calculate  "autocorrelation"
          CALL CONVL(MADRID(IVEC1),1,NPIX(2),1,1,MADRID(IVEC2),WINDOW)
           IF (PLOT(1:1) .EQ. 'Y' .OR. PLOT(1:1) .EQ. 'y') THEN
               K = NPIX(2)/2-WINDOW
               CALL PLTARR(MADRID(IVEC2),NPIX(2),K,K+WINDOW*2,FRMIN,1)
           ENDIF
C ... determine order width
           ORDWID = LWORD(MADRID(IVEC2),NPIX(2),NPIX(2)/2,AM,XC,HWHI)
           TMP(I) = HWHI
        ENDIF
        ENDIF
      ENDDO
C
      IF ((TMP(-1) .LT. TMP(0)) .AND. (TMP(-1) .LT. TMP(1))) THEN
        IDIR   = -1
        TMP(1) = TMP (0)
        TMP(0) = TMP (-1)
        XINCL  = INCL - 2*SINCL
      ELSE IF ((TMP(1) .LT. TMP(0)) .AND. (TMP(1) .LT. TMP(-1))) THEN
        IDIR   =  1
        TMP(-1)= TMP (0)
        TMP(0) = TMP (1)
        XINCL  = INCL
      ELSE
        IDIR   = 0
        FLAG   = 1
        SINCL  = SINCL/2
        WRITE (*,*) 'New step =',SINCL
      ENDIF
C
C      WRITE (LINE,9000) TMP(0)*2
C      CALL STTPUT(LINE,ISTAT)
      WRITE (LINE,9010) XINCL
      CALL STTPUT(LINE,ISTAT)
      CALL STTPUT(' ------------------------------',ISTAT)
      IF (SINCL .GT. ACCUR) GOTO 10
C
      CALL STKWRR('OUTPUTR',XINCL,1,1,KUN,STATUS)
      CALL STKWRR('OUTPUTR',TMP(0)*2,2,1,KUN,STATUS)
      CALL STTPUT('--------------------',ISTAT)
      WRITE (LINE,9020) TMP(0)*2
      CALL STTPUT(LINE,ISTAT)
      WRITE (LINE,9030) XINCL
      CALL STTPUT(LINE,ISTAT)
C ... that's all ...
C ... free data
      CALL STSEPI
 9000 FORMAT ('OrderWidth: ',F8.4)
 9010 FORMAT ('Inclination_ of_Order: ',F8.4)
 9020 FORMAT ('Min_Order_Width: ',F8.4)
 9030 FORMAT ('Optimum_Inclination_of_Order: ',F8.4)
 9040 FORMAT ('Center_Column: ',I6)
 9050 FORMAT ('Half_of_width: ',I6)

      END

C*****************************************************************
      SUBROUTINE EXTRAV(X,N1,N2,I,WIDTH,INCL,OUT)
C
C       Extract (with averaging) 2*WIDHT+1 vertical columns
C       from image array X(N1,N2) centered at I-column
C
      INTEGER N1            ! IN
      INTEGER N2            ! IN
      REAL    X(N1, N2)     ! IN - input image
      INTEGER I             ! IN - position of center column
      INTEGER WIDTH         ! IN - +- from I-column
      REAL    INCL          ! IN - orders inclination coefficient
      REAL    OUT(N2)       ! OUT- output trace
C
      INTEGER J, J1, K, K1
C
      DO 10 J = 1, N2
         OUT(J) = 0.
         DO 20 K = I-WIDTH,I+WIDTH
            J1 = MAX(1, MIN(N2, INT(J + (K-I)*INCL + 0.5)))
            IF (K .LT. 1)  THEN
                K1 = -(K-1)
            ELSE IF (K .GT. N1) THEN
                K1 = N1-(K-N1)
            ELSE
                K1 = K
            ENDIF
            OUT(J) = OUT(J) + X(K1,J1)
20       CONTINUE
         OUT(J) = OUT(J)/(2*WIDTH+1)
10    CONTINUE
      RETURN
      END


C*****************************************************************
      SUBROUTINE CONVL(X,N1,N2,I1,I2,OUT,NSTEP)
C
C       Perform "convolution" of two vertical columns: I1, I2
C       of image array X(N1,N2)
C
      INTEGER N1                ! IN
      INTEGER N2                ! IN
      REAL    X(N1, N2)         ! IN
      INTEGER I1, I2            ! IN - NUMERICS OF COLUMN
      REAL    OUT(N2)           ! OUT
      INTEGER NSTEP             ! IN
C
      INTEGER J, K, N, N22
      REAL S, S1
C
      N22 = N2/2
      N = MIN0(NSTEP,N22-1)
      S1 = 0.
      DO 10 K = 1, N2
         OUT(K) = 0.
         S1 = S1 + X(I1,K)
10    CONTINUE
      DO 30 J = -N, N
         S = 0.
         DO 20 K = 1, N2
            S = S + X(I1,K)*X(I2,MOD(N2+K+J-1,N2)+1)
20       CONTINUE
         OUT(N22+J) = S/S1
30    CONTINUE
      RETURN
      END

C*****************************************************************
      INTEGER FUNCTION LWORD(INP,NP,CENTER,AM,XC,HWHI)
C
C      Determine new order width
C
      INTEGER NP                   ! IN arr.length
      REAL    INP(NP)              ! IN input array
      INTEGER CENTER               ! IN center pixel of contour
      REAL    AM, XC, HWHI         ! OUT output params
C
      INTEGER I, K, WND, IC
      REAL    A, B, D, A11, A12, A22, B1, B2
      REAL    X, Y, Y1, Y2
C
C ... at first - rough estimation of HWHI
      K = NP
      DO I = CENTER+1,NP-1
        IF (INP(I) .LE. INP(CENTER)/2.) THEN
            K = I
            GO TO 10
        ENDIF
        IF ((INP(I-1).GE.INP(I)) .AND. (INP(I+1).GE.INP(I))) THEN
            K = I-2
            GO TO 10
        ENDIF
      ENDDO
C
10    WND = K - CENTER
      LWORD = WND
C
C ... and then try to determine parameters with more precision
      AM = INP(CENTER)
      IC = CENTER
C ... make sure that we stay at maximum
      DO I = CENTER-WND+1, CENTER+WND
        IF (INP(I) .GT. AM) THEN
            AM = INP(I)
            IC = I
        ENDIF
      ENDDO
C ... redefine window of interest
      Y1 = AM
      Y  = (INP(IC+1) + INP(IC-1))/2.
      DO I = 1, WND
        Y2  = (INP(IC+I+1) + INP(IC-I-1))/2.
        IF (Y1-Y .GE. Y-Y2) THEN
            WND = I+1
            GO TO 20
        ENDIF
        Y1 = Y
        Y  = Y2
      ENDDO
C ... perform very simple fit
20    A11 = 0.
      A12 = 0.
      A22 = WND*2.
      B1 = 0.
      B2 = 0.
      DO I = IC-WND+1, IC+WND
        AM = MAX(AM, INP(I))
        X = I - 0.5
        Y = INP(I) - INP(I-1)
        A11 = A11 + X**2
        A12 = A12 + X
        B1  = B1  + X*Y
        B2  = B2  + Y
      ENDDO
      D = (A11*A22 - A12*A12)
      A = (B1*A22  - B2*A12)/D
      B = (B2*A11  - B1*A12)/D
      XC = -B/A
      HWHI = SQRT( MAX(0., -AM/A))
      WRITE (*,9000) AM, XC, HWHI*2.
C
      RETURN
 9000 FORMAT ('Am=',F7.1,'  Xc=',F6.1,'  FWHI=',F6.3)
      END

C*****************************************************************

      SUBROUTINE PLTARR(OUT,N2,NB,NE,NAME,OVR)
C
C       Plot elemets from NB to NE of array OUT(N2) with identifier NAME
C       (for debuging)
C
      INTEGER N2                    ! total length of data array
      REAL    OUT(N2)               ! array with data to plot
      INTEGER NB, NE                ! edges of plotting part of array
      INTEGER OVR                   ! overplot?
      CHARACTER*(*) NAME            ! label
C
      INTEGER   ISTAT, IAC, KUN, KNUL
      INTEGER   FMODE, IMAGE(4), NPL
      REAL      XS, SX, YOFF
      REAL      CUTS(2)
      REAL      FRAME(8)
      REAL      MUDAKI(2)
      CHARACTER LABEL1*80, LABEL2*80, LABEL3*80, LABEL4*80
      CHARACTER XFRAME*4, YFRAME*4
C
      INCLUDE 'MID_INCLUDE:PLTDEC.INC/NOLIST'
      INCLUDE 'MID_INCLUDE:PLTDAT.INC/NOLIST'
C
      DATA PMODE/-1/
      DATA MUDAKI/0.,0./
C
      CALL STKWRR('PLRSTAT',MUDAKI,19,2,KUN,ISTAT)
C
      CALL STKRDI('PLISTAT',1,1,IAC,FMODE,KUN,KNUL,ISTAT)
      FMODE = MIN0(FMODE,1)       ! we need not too complicated legenda
      CALL STKRDC('PLCSTAT',1,5,4,IAC,XFRAME,KUN,KNUL,ISTAT)
      CALL STKRDC('PLCSTAT',1,9,4,IAC,YFRAME,KUN,KNUL,ISTAT)
      IF (PMODE.LT.0 .OR. OVR.EQ.0) THEN
         PMODE = 0
      ELSE
         PMODE = 2
      ENDIF
C
      IMAGE(1) = MAX( 1, NB)
      IMAGE(2) = MIN(N2, NE)
      IMAGE(3) = 1
      IMAGE(4) = 1
      NPL = MAX( 2,  IMAGE(2) - IMAGE(1) + 1 )
      XS = 1.
      SX = 1.
C
C *** only for plot command
      IF (PMODE.EQ.0) THEN
         YOFF = 0.0
C
C *** find min and max of the data values
C         CALL MINMAX(OUT,N2,CUTS(1),CUTS(2))
C
C*** get the frame parameters
         IF (XFRAME(1:4).EQ.'MANU') THEN                        
            CALL STKRDR('PLRSTAT',11,4,IAC,FRAME(1),KUN,KNUL,ISTAT)
C            CALL BOXWTP(FRAME(1),NPL,XS,SX,IMAGE(1))
         ELSE
            XFRAME   = 'AUTO'
            FRAME(1) = FLOAT(IMAGE(1))
            FRAME(2) = FLOAT(IMAGE(2))
         END IF
C
         IF (YFRAME(1:4).EQ.'MANU') THEN
             CALL STKRDR('PLRSTAT',15,4,IAC,FRAME(5),KUN,KNUL,ISTAT)
         ELSE
             YFRAME   = 'AUTO'
             FRAME(5) = CUTS(1)
             FRAME(6) = CUTS(2)
         END IF

C         CALL GETFRM(FRAME(1),FRAME(2),FRAME(3),FRAME(4),XFRAME)
C         CALL GETFRM(FRAME(5),FRAME(6),FRAME(7),FRAME(8),YFRAME)
C
C        WRITE(*,*) 'XFRAME=', XFRAME
C        WRITE(*,*) 'YFRAME=', YFRAME
C        WRITE(*,*) 'FRAME=', FRAME(1), FRAME(2), FRAME(3), FRAME(4)
C        WRITE(*,*) 'FRAME=', FRAME(5), FRAME(6), FRAME(7), FRAME(8)
C
         CALL STKWRR('PLRSTAT',FRAME,11,8,KUN,ISTAT)
C
C     ELSE                                                       ! overplot mode
C        YOFF = YOFF + ABS(CUTS(1) - CUTS(2))/10.
C        CALL STKRDR('INPUTR',1,1,IAC,YOFF,KUN,KNUL,ISTAT)
      ENDIF
C
C *** do the plot setup
C      CALL PLOPN(NAME,PMODE,FMODE)
C
C *** do the work
C      CALL PLKEY(OUT(IMAGE(1)),NPL,IMAGE,XS,SX,YOFF)
C
C *** plot the legenda
      IF (FMODE.EQ.1 .AND. PMODE.EQ.0) THEN
C         CALL PLFRAM(FRAME(1),FRAME(2),FRAME(3),FRAME(4),
C     2               FRAME(5),FRAME(6),FRAME(7),FRAME(8))
         LABEL1 = 'Sequence Number'
         LABEL2 = 'Value'
C
         IF (FMODE.EQ.1) THEN
            LABEL3 = 'Array: '//NAME
            LABEL4 = ' '
C            CALL PLLABL(LABEL1,LABEL2,LABEL3,LABEL4)
         END IF
      END IF
C
C
C *** good bye and return
C      CALL PLCLS
C
      RETURN
      END
