C @(#)neczebord.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_ord.for
C
C.KEYWORDS
C
C  echelle, order searching, maxima following
C
C.PURPOSE
C
C Search for echelle(ZEBRA) orders on any echelle image.
C Output is a table with the following columns :
C
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  Maxima following algorithm.
C
C.INPUT/OUTPUT
C
C  COMMAND :
C  SEARCH/ORDER INPUT W OUTPUT
C  input keywords are
C  ECHC(1:8)   INPUT   - IMAGE NAME
C  ECHC(21:28)  OUTPUT  - TABLE NAME
C         W - PARAMETER WHERE :
C  ECHR(1)              average widht of orders
C
C  OUTPUT PARAMETRES:
C
C-----------------------------------------------------------------
      PROGRAM ZEBORD

      IMPLICIT NONE
      INTEGER MAXORD
      PARAMETER (MAXORD = 100)
      INTEGER MADRID
      INTEGER ACTVAL,I,ISTAT,IXSTR,COLUMN(3)
      INTEGER NACOL,NAXIS,NDIM,NORDER
C     INTEGER NPOINT
      INTEGER IVEC1, IVEC2, IVEC3, IVEC4, ITBL, NTBL
      INTEGER NTOT,PNTRA,STATUS,WINDOW,IACOL
      INTEGER NPIX(3),ICOL(7), KUN, KNUL, IMNO, IJ
      INTEGER UPPER(MAXORD),LOWER(MAXORD)
      REAL    FORD(MAXORD), RPAR(2), CORD(MAXORD)
      REAL    INCLIN
      INTEGER ISTEP, IS, ISTART, ISTOP, IDIR, TID
      INTEGER UP0(MAXORD),LOW0(MAXORD)
      REAL    CORD0(MAXORD)
      DOUBLE PRECISION START(3),STEP(3)
      DOUBLE PRECISION ST
      CHARACTER*80 CA
      CHARACTER*60 TABLE
      CHARACTER*72 IDENA
      CHARACTER*80 LINE
      CHARACTER*64 CUNITA
      CHARACTER*16 LABCOL(5),UNIT(5)
      CHARACTER*6  FORM(5)
      CHARACTER*8  FRMIN
      CHARACTER*1  PLOT
      LOGICAL      DBGPLT
C                   WIDTH  - HALF 0F WIDTH FOR AVERAGE
      INTEGER WIDTH, NTRACE, NTR, KX, KY
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/
      DATA PLOT/'N'/

      CA = 'COUNTS'
      ST = 1.D0
      NDIM   = 3
C
C ... initialize system and read parameters
C
      CALL STSPRO('ZEBORD')
      CALL STKRDC('IN_B',1,1,60,ACTVAL,TABLE,KUN,KNUL,STATUS)
      CALL STKRDC('IN_A',1,1,60,ACTVAL,FRMIN,KUN,KNUL,STATUS)
      CALL STKRDC('INPUTC',1,1,1,ACTVAL,PLOT,KUN,KNUL,STATUS)
      CALL STKRDI('INPUTI',1,3,ACTVAL,COLUMN,KUN,KNUL,STATUS)
      CALL STKRDR('INPUTR',1,2,ACTVAL,RPAR,KUN,KNUL,STATUS)
      IXSTR   = COLUMN(1)
      WIDTH   = COLUMN(2)
      NTRACE  = COLUMN(3)
      INCLIN  = RPAR(1)
      WINDOW  = RPAR(2) + 0.5
      DBGPLT  = (PLOT.EQ.'Y') .OR. (PLOT.EQ.'y')
C
      CALL STTPUT(' Echelle(ZEBRA) definition',ISTAT)
      CALL STTPUT(' ------------------',ISTAT)
      CALL STTPUT(' InputImage  : '//FRMIN,ISTAT)
      CALL STTPUT(' OutputTable : '//TABLE,ISTAT)
      CALL STTPUT(' Parameters  ',ISTAT)
      WRITE (LINE,9040) IXSTR
      CALL STTPUT(LINE ,ISTAT)
      WRITE (LINE,9050) WIDTH
      CALL STTPUT(LINE ,ISTAT)
      WRITE (LINE,9020) WINDOW
      CALL STTPUT(LINE,ISTAT)
      WRITE (LINE,9030) INCLIN
      CALL STTPUT(LINE,ISTAT)
      WRITE (LINE,9010) NTRACE
      CALL STTPUT(LINE,ISTAT)
      CALL STTPUT(' PLOT FLAG  : '//PLOT,ISTAT)
 9010 FORMAT ('Follow_Steps_Number: ',I6)
 9020 FORMAT ('Order_Width: ',I6)
 9030 FORMAT ('Order_Inclination: ',F8.4)
 9040 FORMAT ('Center_Column: ',I6)
 9050 FORMAT ('Averagin_Width: ',I6)
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('ZEB_ORD.WRK',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,1,
     .            NPIX(2)*4 + MAXORD*NTRACE*3,
     .            ST,ST,CA,CA,IVEC1,IJ,STATUS)
      IVEC2 = IVEC1 + NPIX(2)
      IVEC3 = IVEC2 + NPIX(2)
      IVEC4 = IVEC3 + NPIX(2)
      ITBL  = IVEC4 + NPIX(2)
C
C ... extract and plot "averaging" column
      CALL EXTRAV(MADRID(PNTRA),NPIX(1),NPIX(2),IXSTR,WIDTH,INCLIN,
     .           MADRID(IVEC1))
      IF ( DBGPLT ) THEN
        CALL PLTARR(MADRID(IVEC1),NPIX(2),100,NPIX(2)-150,FRMIN,0)
      ENDIF
C
C ... find order positions in the middle of the image by BGD-thresholding
C
      CALL ORDTHR(MADRID(IVEC1),MADRID(IVEC2),MADRID(IVEC3),
     .           MADRID(IVEC4),NPIX(2),WINDOW,UPPER,LOWER,FORD,
     .           CORD,NORDER,FRMIN,DBGPLT)
C
CD... simulate order miss ( for ORDCHK() debuging )
CD      NORDER = NORDER - 1
CD      DO I = 2, NORDER
CD          UPPER(I) = UPPER(I+1)
CD          LOWER(I) = LOWER(I+1)
CD          CORD(I)  = CORD(I+1)
CD      ENDDO
C
C ... try to find missing orders
      CALL ORDCHK(NORDER,UPPER,LOWER,CORD,MADRID(ITBL),MADRID(IVEC3)
     .           ,NPIX(2),FRMIN,DBGPLT)
C
C ... allocate table
C
      NTOT   = (NORDER+1)*NTRACE
      IACOL  = 12
      CALL TBTINI(TABLE,F_TRANS,F_O_MODE,IACOL,NTOT,TID,STATUS)
      DO 10 I = 1,NACOL
          CALL TBCINI(TID,D_R4_FORMAT,1,FORM(I),
     .    UNIT(I),LABCOL(I),ICOL(I),STATUS)
 10   CONTINUE
C
C ... now we need to follow orders
      ISTEP = NPIX(1)/NTRACE
      NTBL = 1
      NTR  = 0
      KY   = MAX(0, MIN(3, (NORDER+3)/5))
      DO IDIR = -1,1,2
        DO I = 1,NORDER
            UP0(I) = UPPER(I)
            LOW0(I) = LOWER(I)
            CORD0(I) = CORD(I)
        ENDDO
        IF (IDIR.GT.0) THEN
            CALL REVTBL(MADRID(ITBL),NTBL)
            ISTART = ((IXSTR/ISTEP) + 1) * ISTEP + 1
            ISTOP  = NPIX(1) - ISTEP
        ELSE
            ISTART = (IXSTR/ISTEP) * ISTEP + 1
            ISTOP  = ISTEP + 1
        ENDIF
        IS = ISTART - IXSTR
        DO I = ISTART, ISTOP, ISTEP*IDIR
            WRITE(*,*) '========================================='
            WRITE(*,*) 'AvrColunm=',I
            NTR = NTR + 1
            KX = MAX(0, MIN(2, NTBL/8/(KY+1)-1 ))
C
C23456789012345678901234567890123456789012345678901234567890123456789012345
C
            CALL ORDFIT(MADRID(ITBL),NTBL,KX,KY,I,UP0,LOW0,CORD0,NORDER
     .          ,INCLIN,IS,MADRID(IVEC3),NPIX(1),NPIX(2),FRMIN,.FALSE.)
            IS = ISTEP*IDIR
            CALL EXTRAV(MADRID(PNTRA),NPIX(1),NPIX(2),I,ISTEP,INCLIN,
     .                  MADRID(IVEC1))
CD            READ(*,*) III
            IF ( DBGPLT ) THEN
             CALL PLTARR(MADRID(IVEC1),NPIX(2),100,NPIX(2)-150,FRMIN,0)
            ENDIF
            CALL ORDCOR(MADRID(IVEC1),MADRID(IVEC2),MADRID(IVEC3),
     .                 MADRID(IVEC4),NPIX(2),WINDOW,UP0,LOW0,FORD,
     .                 CORD0,NORDER,FRMIN,DBGPLT)
            CALL ORDTBL(I,FORD,CORD0,NORDER,MADRID(ITBL),NTBL)
        ENDDO
      ENDDO
      CALL WRTTBL(TID,NORDER,MADRID(ITBL),NTBL,NPIX(1))
C
C ... update number of elements
C      NTOT   = NPOINT
      CALL STKWRI('ECHORD',NORDER,1,1,KUN,KNUL,STATUS)
      IDENA = 'ORDER POSITION'
C      CALL SXDPUT(TABLEFULL,'IDENT','C',IDENA,1,72,STATUS)
C      CALL DSCUPD(TABLEFULL,TABLEFULL,' ',STATUS)
CD    CALL TBCSRT(TID,1,1,1,STATUS)
      CALL TBSINI(TID,STATUS)
      CALL TBTCLO(TID,STATUS)
C
C ... that's all ...
C ... free data
   20 CALL STSEPI

      END

C*****************************************************************
      SUBROUTINE ORDTHR(VEC1,VEC2,VEC3,IVEC4,N2,WINDOW,
     .                  UP,LOW,FORD,CORD,NORD,FRMIN,DBGPLT)
C
C       Search for orders in the image trace (BGD-threshholding)
C
      INTEGER MAXORD
      REAL    SIGMA3
      PARAMETER (MAXORD = 100)
      INTEGER N2
      REAL    VEC1(N2)                   ! input trace
      REAL    VEC2(N2),VEC3(N2)          ! working space
      INTEGER IVEC4(N2)                  ! working space
      INTEGER WINDOW                     ! order width (IN/OUT!)
      INTEGER UP(MAXORD)                 ! orders beg. (OUT)
      INTEGER LOW(MAXORD)                ! orders end  (OUT)
      REAL    FORD(MAXORD)               ! avr.flux    (OUT)
      REAL    CORD(MAXORD)               ! ord. center (OUT)
      INTEGER NORD                       ! # of orders (OUT)
      CHARACTER*(*) FRMIN                ! label
      LOGICAL DBGPLT                     ! debug plotting flag
C
      INTEGER I, WND, ISTAT, IB, IE, K, L
C      REAL    AM, XC, HWHI
      REAL    S, S0, S1
      REAL    NOISE0, NOISE1, NOISE2
      REAL    MAXR, MAXL, MINC, Q
      CHARACTER*80 LINE
C
CD... culculate  "autocorrelation"
CD    CALL CONVL(VEC1, 1, N2, 1, 1, VEC2, WINDOW*5)
CD    K = N2/2-WINDOW*5
CD    L = WINDOW*10
CD    CALL PLTARR(VEC2, N2, K, K+L, FRMIN, 1)
C
CD... determine order width
CD    WINDOW = LWORD(VEC2, N2, N2/2, AM, XC, HWHI)
CD    WRITE (LINE,9000) WINDOW*2
CD    CALL STTPUT(LINE,ISTAT)
C
C ... estimate noise level
      NOISE0 = SIGMA3(VEC1, N2)
C
C ... extract "background" (variable threshold level)
      WND = WINDOW * 0.75 + 0.5
      I = 0
      CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE0, 4, I)
CD    CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1)
      NOISE1 = SIGMA3(VEC2, N2)
      CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE1, 4, I)
CD    CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1)
      NOISE2 = SIGMA3(VEC2, N2)
CD    CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE2, 4, I)
      IF ( DBGPLT ) THEN
        CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1)
CD        READ(*,*) III
      END IF
      WRITE (LINE,9010)  NOISE0, NOISE1, NOISE2
      CALL STTPUT(LINE,ISTAT)
C
C ... use threshold level to form orders mask
      IVEC4(1) = 0
      IVEC4(N2) = 0
      DO I = 2, N2-1
        VEC2(I) = VEC1(I)-VEC2(I)
        IF ( VEC1(I) .GT. VEC3(I)+NOISE1) THEN
            IVEC4(I) = 1
        ELSE
            IVEC4(I) = 0
        ENDIF
        VEC3(I) = IVEC4(I)
      ENDDO
      IF ( DBGPLT ) THEN
        CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1)
CD        READ(*,*) III
      ENDIF
C
C ... try to find orders
      WRITE (LINE, *) '_No._Begin__End__AvrFlux_Center_______'
      CALL STTPUT(LINE,ISTAT)
      NORD = 0
      DO I = 1, N2-1
        VEC3(I) = 0.
        IF (IVEC4(I).EQ.0 .AND. IVEC4(I+1).NE.0) THEN
            IB = I+1
            S  = 0.
            S0 = 0.
            S1 = 0.
        ELSEIF (IVEC4(I).NE.0 .AND. IVEC4(I+1).NE.0) THEN
            S  = S  + VEC2(I)
            S0 = S0 + VEC1(I)
            S1 = S1 + VEC1(I)*I
        ELSEIF (IVEC4(I).NE.0 .AND. IVEC4(I+1).EQ.0) THEN
            IE = I
            IF (IE.GT.IB) THEN
                S1 = S1/S0
                S0 = S0/(IE-IB+1)
            ENDIF
            IF (IE-IB+1 .GE. WINDOW/2) THEN
                Q = (IE-IB)/3.0             !    |  Q  |     |  Q  |
                K = IB+Q+0.5                !    |     |     |     |
                L = IE-Q+0.5                ! ..IB.....K.....L.....IE..
                CALL MINMAX(VEC2(IB), K-IB,  Q, MAXL)
                CALL MINMAX(VEC2(K),  L-K+1, MINC, Q)
                CALL MINMAX(VEC2(L+1),IE-L,  Q, MAXR)
                IF ((IE-IB+1.0 .GE. WINDOW*1.4) .AND.
     .              (MAXR.GT.MINC .AND. MAXL.GT.MINC)) THEN
                    IVEC4((IE+IB)/2) = 0
                    I = IB-2
CD                  WRITE (*,*) IB,K,L,IE
CD                  WRITE (*,*) MAXL,MINC,MAXR
                    WRITE (LINE,9040) IB, IE, S0, S1
                    CALL STTPUT(LINE,ISTAT)
C
C23456789012345678901234567890123456789012345678901234567890123456789012345
C
                    WRITE (LINE, *) 
     .              '-- Warning! Problems with order detection! -- '
                    CALL STTPUT(LINE,ISTAT)
                ELSEIF (S .GT. WINDOW*NOISE1*2.) THEN
                    IF (NORD.LT.MAXORD) THEN
                        NORD = NORD + 1
                    ELSE
                        GO TO 100
                    ENDIF
                    UP(NORD) = IB
                    LOW(NORD) = IE
                    FORD(NORD) = S0
                    CORD(NORD) = S1
                    WRITE (LINE,9020) NORD, IB, IE, S0, S1
                    CALL STTPUT(LINE,ISTAT)
                    VEC3( INT(S1+0.5) ) = FORD(NORD)
                ELSE
                    GO TO 100
                ENDIF
            ELSE
100             CONTINUE
                WRITE (*,9030) IB, IE, S0, S1
            ENDIF
        ENDIF
      ENDDO
      WRITE(*,*) 'NumOrders=',NORD
CD    READ(*,*) I
      IF ( DBGPLT ) THEN
        CALL PLTARR(VEC1, N2, 100, N2-150, FRMIN, 0)
CD        READ(*,*) III
        CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1)
CD        READ(*,*) III
      ENDIF
C
      RETURN
 9000 FORMAT ('OrderWidth: ',I6)
 9010 FORMAT ('NoiseLevels: ', 3F8.3)
 9020 FORMAT (I4,2I6,2F8.1)
 9030 FORMAT (' ---',2I6,2F8.1)
 9040 FORMAT (' ***',2I6,2F8.1)
      END

C*****************************************************************
      SUBROUTINE ORDCHK(NORD,UP,LOW,CORD,TBL,VEC,N2,FRMIN,DBGPLT)
C
C       Check orders positions (to find missing orders)
C
      INTEGER MAXORD
      REAL    EVAPOL
      PARAMETER (MAXORD = 100)
      INTEGER NORD                       ! # of orders (IN/OUT)
      INTEGER UP(MAXORD)                 ! orders beg. (IN/OUT)
      INTEGER LOW(MAXORD)                ! orders end  (IN/OUT)
      REAL    CORD(MAXORD)               ! ord. center (IN/OUT)
      REAL    TBL(3,MAXORD)              ! wrk.space
      INTEGER N2
      REAL    VEC(N2)                    ! wrk.space
      CHARACTER*(*) FRMIN                ! label
      LOGICAL DBGPLT                     ! debug plotting flag
C
      INTEGER I, ISTAT, J, K, L, LM, NEW
      REAL    COEF(30), ERR, ERR0, ERR1, ERR2, CL, C1, C2
      CHARACTER*80 LINE
      LOGICAL OK
C
C ... perform polynomial fitting of order positions
      OK = .FALSE.
  10  DO I = 1,NORD
        TBL(1,I) = CORD(I)
        TBL(2,I) = I
        TBL(3,I) = 0.
      ENDDO
C     K   = MAX(0, MIN(3, (NORDER+3)/5))
      IF (NORD .LE. 4) THEN
        RETURN
      ELSE IF (NORD .LE. 6) THEN
        K = 1
      ELSE IF (NORD .LE. 11) THEN
        K = 2
      ELSE
        K = 3
      ENDIF
      CALL POLFIT(TBL, NORD, K, 0, COEF, ERR)
      ERR0 = SQRT(ERR**2/NORD)
C
      IF ( OK ) GO TO 100
C
C ... for debuging plot:
CD    WRITE (*, 9000) K, 0
CD    WRITE (*, 9010) (COEF(I),I=1,K+1)
CD    WRITE (*, 9020) ERR0
CD    READ (*,*) I
      IF ( DBGPLT ) THEN
        J = N2/NORD
        DO I = 1,N2
          VEC(I) = CORD( MIN(NORD, (I - 1)/J + 1) )
        ENDDO
        CALL PLTARR(VEC, N2, 1, N2, FRMIN, 0)
        DO I = 1,N2
          VEC(I) = EVAPOL(COEF, K, 0, FLOAT(I-1)/J+0.5, 0.)
        ENDDO
        CALL PLTARR(VEC, N2, 1, N2, FRMIN, 1)
      ENDIF
C
C ... try to insert new orders between every pair of old orders
      NEW = 0
      DO L = 2,NORD
C ...   at first - one order
        DO I = 1,NORD+1
            IF (I .LT. L) THEN
                TBL(1,I) = CORD(I)
            ELSE IF (I .GT. L) THEN
                TBL(1,I) = CORD(I-1)
            END IF
            TBL(2,I) = I
            TBL(3,I) = 0.
        ENDDO
        TBL(1,L) = (CORD(L) + CORD(L-1)) / 2.0
        CALL POLFIT(TBL, NORD+1, K, 0, COEF, ERR)
        ERR1 = SQRT(ERR**2/(NORD+1))
        CL = EVAPOL(COEF, K, 0, FLOAT(L), 0.)
C ...   then - two (if posible)
        IF (NORD .GE. 7) THEN
            DO I = 1,NORD+2
                IF (I .LT. L) THEN
                    TBL(1,I) = CORD(I)
                ELSE IF (I .GT. L+1) THEN
                    TBL(1,I) = CORD(I-2)
                END IF
                TBL(2,I) = I
                TBL(3,I) = 0.
            ENDDO
            TBL(1,L) = CORD(L-1) + (CORD(L) - CORD(L-1)) * 0.333
            TBL(1,L+1) = CORD(L-1) + (CORD(L) - CORD(L-1)) * 0.667
            CALL POLFIT(TBL, NORD+2, K, 0, COEF, ERR)
            ERR2 = SQRT(ERR**2/(NORD+2))
        ELSE
            ERR2 = ERR1*2.
        END IF
CD      WRITE (*, 9030) L-1, L, ERR1, ERR2
C
        IF (ERR1.LT.ERR0 .OR. ERR2.LT.ERR0) THEN
            IF (ERR1 .LE. ERR2) THEN
                ERR0 = ERR1
                NEW = 1
                C1 = CL
            ELSE
                ERR0 = ERR2
                NEW = 2
                C1 = EVAPOL(COEF, K, 0, FLOAT(L), 0.)
                C2 = EVAPOL(COEF, K, 0, FLOAT(L+1), 0.)
            END IF
            LM = L
        END IF
      ENDDO
C
C ... if there are new orders, fix them
      IF (NEW .GT. 0) THEN
        WRITE (LINE, 9040) NEW, LM-1, LM
        CALL STTPUT(LINE,ISTAT)
        DO I = NORD+NEW, LM+NEW, -1
            CORD(I) = CORD(I-NEW)
            UP(I)   = UP(I-NEW)
            LOW(I)  = LOW(I-NEW)
        ENDDO
        CORD(LM) = C1
        L = (LOW(LM-1) - UP(LM-1) + LOW(LM) - UP(LM))/2 + 1
        LOW(LM) = MIN(INT(CORD(LM) + L/2.0 + 0.5), UP(LM)-1)
        UP(LM)  = MAX(INT(CORD(LM) - L/2.0 + 0.5), LOW(LM-1)+1)
        IF (NEW .GT. 1) THEN
            CORD(LM+1) = C2
            UP(LM+1)  = MAX(INT(CORD(LM+1) - L/2.0 + 0.5), LOW(LM)+1)
            LOW(LM+1) = MIN(INT(CORD(LM+1) + L/2.0 + 0.5), UP(LM+2)-1)
        END IF
        NORD = NORD+NEW
      ELSE
        OK = .TRUE.
      END IF
      GO TO 10
C
100   DO I = 1, NORD
        CORD(I) = EVAPOL(COEF, K, 0, FLOAT(I), 0.)
      ENDDO
C
      RETURN
CD9000 FORMAT ('Degrees: ',2I2)
CD9010 FORMAT ('Coef-s: ',3F8.3)
CD9020 FORMAT ('Error: ',F7.3)
CD9030 FORMAT ('Ord:',I2,'-',I2,' Err: ',2F6.3)
 9040 FORMAT ('-- Attention! ',I1,'-ord.miss.in:',I2,'-',I2)
      END

C*****************************************************************
      SUBROUTINE ORDCOR(VEC1,VEC2,VEC3,VEC4,N2,WINDOW,
     .                  UP,LOW,FORD,CORD,NORD,FRMIN,DBGPLT)
C
C       Correct orders positions (by correlation with model contour)
C
      INTEGER MXBEAM, MAXORD
      REAL    SIGMA3
      PARAMETER (MAXORD = 100)
      PARAMETER (MXBEAM = 100)
      INTEGER N2
      REAL    VEC1(N2)                   ! input trace
      REAL    VEC2(N2),VEC3(N2),VEC4(N2) ! working space
      INTEGER WINDOW                     ! order width (IN)
      INTEGER UP(MAXORD)                 ! orders beg. (IN/OUT)
      INTEGER LOW(MAXORD)                ! orders end  (IN/OUT)
      REAL    FORD(MAXORD)               ! avr.flux    (OUT)
      REAL    CORD(MAXORD)               ! ord. center (IN/OUT)
      INTEGER NORD                       ! # of orders (IN)
      CHARACTER*(*) FRMIN                ! label
      LOGICAL DBGPLT                     ! debug plotting flag
C
      INTEGER I, WND, J, NBEAM, K, KMAX, IMAX, L
      INTEGER IB, IE, PIX(2)
      REAL    S, SS, SX, SB, SMAX, CUTS(2)
      REAL    NOISE0, NOISE1, NOISE2
      REAL    BEAM(MXBEAM)
C
      DATA BEAM /MXBEAM*0./
C
C ... estimate noise level
      NOISE0 = SIGMA3(VEC1, N2)
C
C ... extract and subtract "background"
      WND = WINDOW * 0.75 + 0.5
      I = 0
      CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE0, 4, I)
      NOISE1 = SIGMA3(VEC2, N2)
      CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE1, 4, I)
      NOISE2 = SIGMA3(VEC2, N2)
C     CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE2, 4, I)
      IF ( DBGPLT ) THEN
        CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1)
      ENDIF
CD    WRITE (LINE,9010)  NOISE0, NOISE1, NOISE2
CD    CALL STTPUT(LINE,ISTAT)
      DO I = 1, N2
        VEC2(I) = VEC1(I)-VEC3(I)
      ENDDO
C
C ... construct sample order profile (we use gaussian)
      NBEAM =  MIN(WINDOW,MXBEAM)
      BEAM(1) = 1.0
      SB = 1.0
      DO I = 2, NBEAM
        BEAM(I) = EXP(-FLOAT(I-1)**2 / (WINDOW/2.0)**2)
        SB = SB + 2*BEAM(I)**2
      ENDDO
CD    WRITE(*,9030) (BEAM(I),I=1,NBEAM)
C
C ... for debuging plot:
      IF ( DBGPLT ) THEN
        DO I = 1,N2
            VEC4(I) = 0.
        ENDDO
      END IF
C
C ... for every known order...
      DO I = 1, NORD
C ...   search for max.correlation with model order within UP and LOW
        SMAX = -1.
        KMAX = CORD(I)
C ...   try to modify UP and LOW edges if need
        IF (I .GT. 1) THEN
            IB = LOW(I-1) + 1
        ELSE
            IB = 2
        ENDIF
        K = UP(I)
        DO J = UP(I), IB, -1
            IF (VEC2(J-1) .GE. VEC2(J)) GO TO 100
            K = J-1
        ENDDO
100     IB = K
        IF (I .LT. NORD) THEN
            IE = UP(I+1) - 1
        ELSE
            IE = N2-1
        ENDIF
        K = LOW(I)
        DO J = LOW(I), IE
            IF (VEC2(J+1) .GE. VEC2(J)) GO TO 200
            K = J+1
        ENDDO
200     IE = K
        SS = 0.
        SX = 0.
        DO K = IB, IE
C ...       calc.next val. of conv-n with BEAM
            S = VEC2(K)*BEAM(1)
            DO J = 2, NBEAM
                L = K - (J-1)
                IF (L .GE. IB) THEN
                    S = S + VEC2(L)*BEAM(J)
                ELSE
                    S = S + VEC2(IB)*BEAM(IB-L+1)*BEAM(J)
                ENDIF
                L = K + (J-1)
                IF (L .LE. IE) THEN
                    S = S + VEC2(L)*BEAM(J)
                ELSE
                    S = S + VEC2(IE)*BEAM(L-IE+1)*BEAM(J)
                ENDIF
            ENDDO
            S = S/SB
            SS = SS + S
            SX = SX + K*S
            IF (S .GT. SMAX) THEN
                SMAX = S
                KMAX = K
            ENDIF
            IF ( DBGPLT ) VEC4(K) = S
        ENDDO
C ... check is it significant or not?
        CALL MNMX(VEC2(IB), IE-IB+1, CUTS, PIX)
        IMAX = PIX(2)-1 + IB
CD      WRITE (*,*) IB,IE,IMAX,KMAX,SMAX
        IF ( (SMAX .GT. 5.*NOISE1) .AND.
     .       (IMAX .GT. IB) .AND. (IMAX .LT. IE) .AND.
     .       (KMAX .GT. IB) .AND. (KMAX .LT. IE) )  THEN
            CORD(I) = SX / SS
            FORD(I) = SMAX
        ELSE
            FORD(I) = -1.
        ENDIF
      ENDDO
      IF ( DBGPLT ) THEN
        CALL PLTARR(VEC4, N2, 100, N2-150, FRMIN, 1)
      ENDIF
C
C ... for debuging plot
      IF ( DBGPLT ) THEN
        DO I = 1,N2
            VEC4(I) = 0.
        ENDDO
        DO I = 1, NORD
            IF (FORD(I) .GT. 0.) THEN
                VEC4( INT(CORD(I)+0.5) ) = FORD(I)
            ENDIF
        ENDDO
        CALL PLTARR(VEC4, N2, 100, N2-150, FRMIN, 1)
      ENDIF
C
C ... as the last step correct edges location
CD    WRITE (LINE, *) '_No._Begin__End__AvrFlux_Center_______'
CD    CALL STTPUT(LINE,ISTAT)
      DO I = 1, NORD
        IF (I .GT. 1) THEN
            J = (CORD(I-1) + CORD(I))/2.0 + 1.5
        ELSE
            J = 2
        ENDIF
        IF (I .LT. NORD) THEN
            K = (CORD(I+1) + CORD(I))/2.0 + 0.5
        ELSE
            K = N2-1
        ENDIF
        L = LOW(I) - UP(I) + 1
        UP(I)  = MAX(INT(CORD(I)+0.5) - L/2, J)
        LOW(I) = MIN(INT(CORD(I)+0.5) + (L-1)/2, K)
CD      WRITE (LINE,9020) I, UP(I), LOW(I), FORD(I), CORD(I)
CD      CALL STTPUT(LINE,ISTAT)
      ENDDO
C
      RETURN
CD9010 FORMAT ('NoiseLevels: ', 3F8.3)
CD9020 FORMAT (I4,2I6,2F8.1)
CD9030 FORMAT (3F8.3)
      END

C*****************************************************************
      SUBROUTINE ORDFIT(TBL,NTBL,KX,KY,XPOS,UP,LOW,CORD,NORD,INCL
     .                 ,STEP,VEC,N1,N2,FRMIN,DBGPLT)
C
C       2D-poly.fit of current table
C       eval. edges&centers for new X-position
C
      INTEGER MAXORD
      REAL    EVAPOL
      PARAMETER (MAXORD = 100)
      INTEGER NTBL                       ! # of rows in table (IN)
      REAL    TBL(3,NTBL)                ! "in core" table (IN)
      INTEGER KX, KY                     ! polinomial degrees (IN)
      INTEGER XPOS                       ! X-posit.    (IN)
      INTEGER UP(MAXORD)                 ! orders beg. (IN/OUT)
      INTEGER LOW(MAXORD)                ! orders end  (IN/OUT)
      REAL    CORD(MAXORD)               ! ord.cent.   (IN/OUT)
      INTEGER NORD                       ! # of orders (IN)
      REAL    INCL                       ! inclination (IN/OUT)
      INTEGER STEP                       ! step of XPOS(IN)
      INTEGER N1, N2
      REAL    VEC(N1)                    ! wrk.space
      CHARACTER*(*) FRMIN                ! label
      LOGICAL DBGPLT                     ! debug plotting flag
C
      REAL    COEF(30), ERR, SI
      INTEGER OVR, I, J
C
      IF (KX.GT.0 .AND. KY.GT.0) THEN
        CALL POLFIT(TBL, NTBL, KX, KY, COEF, ERR)
        ERR = SQRT(ERR**2/NTBL)
        WRITE(*,*) 'TblLength:',NTBL
        WRITE(*,*) 'FitDegrees:',KX,KY
        WRITE(*,*) 'FitError:',ERR
        IF ( DBGPLT ) THEN
            OVR = 0
            VEC(1)    = N2-150.
            VEC(N1)   = N2-150.
            VEC(2)    = 100.
            VEC(N1-1) = 100.
        ENDIF
        SI = 0.
        DO I = 1,NORD
            CORD(I) = EVAPOL(COEF, KX, KY, FLOAT(XPOS), FLOAT(I))
            UP(I)   = EVAPOL(COEF, KX, KY, FLOAT(XPOS), I-0.4) + 0.5
            LOW(I)  = EVAPOL(COEF, KX, KY, FLOAT(XPOS), I+0.4) + 0.5
            SI = SI + (EVAPOL(COEF, KX, KY, XPOS+1.0, FLOAT(I)) -
     .                 EVAPOL(COEF, KX, KY, XPOS-1.0, FLOAT(I))) / 2.0
            IF ( DBGPLT ) THEN
                DO J = 2, N1-2
                    VEC(J) = EVAPOL(COEF, KX, KY, FLOAT(J), FLOAT(I))
                ENDDO
                CALL PLTARR(VEC, N1, 1, N1, FRMIN, OVR)
                OVR = 1
            END IF
        ENDDO
        INCL = SI / NORD
        WRITE (*,*) 'NewInclination:',INCL
      ELSE
        WRITE (*,*) 'OldInclination:',INCL
        SI = INCL * STEP
        DO I = 1,NORD
          CORD(I) = CORD(I) + SI
          UP(I)   = UP(I)   + SI + 0.5
          LOW(I)  = LOW(I)  + SI + 0.5
        ENDDO
      END IF
C
      RETURN
      END

C*****************************************************************
      SUBROUTINE ORDTBL(XPOS,FORD,CORD,NORD,TBL,NTBL)
C
C       Put next portion of rows to TBL-array
C
      INTEGER MAXORD
      PARAMETER (MAXORD = 100)
      INTEGER XPOS                       ! X-posit.(:X) (IN)
      REAL    FORD(MAXORD)               ! avr.flux     (IN)
      REAL    CORD(MAXORD)               ! ord.cent.(:Y)(IN)
      INTEGER NORD                       ! # of orders  (IN)
      REAL    TBL(3,1)                   ! "in core" table (IN/OUT)
      INTEGER NTBL                       ! # of rows in table (OUT)
C
      INTEGER NC                         ! curr.row number
      INTEGER I
C
      DATA      NC/1/
C
      DO I = 1, NORD
        IF (FORD(I) .GT. 0.) THEN
            TBL(1,NC) =  CORD(I)
            TBL(2,NC) =  XPOS
            TBL(3,NC) =  I
            NTBL = NC
            NC = NC + 1
        ENDIF
      ENDDO
C
      RETURN
      END

C*****************************************************************
      SUBROUTINE REVTBL(TBL,NTBL)
C
C       Reverse TBL-array
C
      INTEGER NTBL                       ! # of rows in TBL (IN)
      REAL    TBL(3,NTBL)                ! "in core" table (IN/OUT)
C
      REAL RVAL(3)
      INTEGER I, J, K
C
      IF (NTBL .GT. 1) THEN
        DO I = 1, NTBL/2
            K = NTBL-I+1
            DO J = 1, 3
                RVAL(J) = TBL(J,I)
            ENDDO
            DO J = 1, 3
                TBL(J,I) = TBL(J,K)
            ENDDO
            DO J = 1, 3
                TBL(J,K) = RVAL(J)
            ENDDO
        ENDDO
      END IF
C
      RETURN
      END

C*****************************************************************
      SUBROUTINE WRTTBL(TID,NORD,TBL,NTBL,N1)
C
C       Form the output table sorted by :ORDER
C
      INTEGER TID                        ! table id     (IN)
      INTEGER NORD                       ! # of orders  (IN)
      INTEGER NTBL                       ! # of rows in table (IN)
      REAL    TBL(3,NTBL)                ! "in core" table (IN)
      INTEGER N1                         ! image X-size       (IN)
C
      INTEGER NC                         ! curr.row number
      REAL RVAL(5)
      INTEGER I, K, ICOL(5), ISTAT
C
      DATA      ICOL/1,2,3,4,5/
      DATA      NC/1/
C
      DO I = 1, 5
        RVAL(I) = 0.
      ENDDO
C
      DO I = 1, NORD
        RVAL(1) =  I
        RVAL(2) =  1.0
        CALL TBRWRR(TID,NC,2,ICOL,RVAL,ISTAT)
        NC = NC + 1
        DO K = 1, NTBL
            IF ( INT(TBL(3,K)+0.1) .EQ. I) THEN
                RVAL(1) =  TBL(3,K)
                RVAL(2) =  TBL(2,K)
                RVAL(3) =  TBL(1,K)
                CALL TBRWRR(TID,NC,5,ICOL,RVAL,ISTAT)
                NC = NC + 1
            ENDIF
        ENDDO
        RVAL(1) =  I
        RVAL(2) =  N1
        CALL TBRWRR(TID,NC,2,ICOL,RVAL,ISTAT)
        NC = NC + 1
      ENDDO
C
      RETURN
      END

C*****************************************************************
      SUBROUTINE EXTRAV(X,N1,N2,I,WIDTH,INCL,OUT)
C
C       Extract and average 2*WIDHT+1 vertical columns
C       from image array X(N1,N2) centered at I-column
C       ( with order inclination account )
C
      INTEGER N1
      INTEGER N2
      REAL    X(N1, N2)     ! input image
      INTEGER I             ! position of center column
      INTEGER WIDTH         ! +- from I-column
      REAL    INCL          ! orders inclination coefficient
      REAL    OUT(N2)       ! 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

CC*****************************************************************
C      SUBROUTINE CONVL(X,N1,N2,I1,I2,OUT,NSTEP)
CC
CC       Perform "convolution" of two vertical columns: I1, I2
CC       of image array X(N1,N2)
CC
C      INTEGER N1
C      INTEGER N2
C      REAL    X(N1, N2)
C      INTEGER I1, I2
C      REAL    OUT(N2)
C      INTEGER NSTEP
CC
C      INTEGER J, K, N, N22
C      REAL S, S1
CC
C      N22 = N2/2
C      N = MIN0(NSTEP,N22-1)
C      S1 = 0.
C      DO 10 K = 1, N2
C         OUT(K) = 0.
C         S1 = S1 + X(I1,K)
C10    CONTINUE
C      DO 30 J = -N, N
C         S = 0.
C         DO 20 K = 1, N2
C            S = S + X(I1,K)*X(I2,MOD(N2+K+J-1,N2)+1)
C20       CONTINUE
C         OUT(N22+J) = S/S1
C30    CONTINUE
C      RETURN
C      END
C
CC*****************************************************************
C      INTEGER FUNCTION LWORD(INP,NP,CENTER,AM,XC,HWHI)
CC
CC      Determine new order width
CC
C      INTEGER NP                   ! IN arr.length
C      REAL    INP(NP)              ! IN input array
C      INTEGER CENTER               ! IN center pixel of contour
C      REAL    AM, XC, HWHI         ! OUT output params
CC
C      INTEGER I, J, K, WND, IC
C      REAL    A, B, D, A11, A12, A22, B1, B2
C      REAL    X, Y, Y1, Y2
CC
CC ... at first - rough estimation of HWHI
C      K = NP
C      DO I = CENTER+1,NP-1
C        IF (INP(I) .LE. INP(CENTER)/2.) THEN
C            K = I
C            GO TO 10
C        ENDIF
C        IF ((INP(I-1).GE.INP(I)) .AND. (INP(I+1).GE.INP(I))) THEN
C            K = I-2
C            GO TO 10
C        ENDIF
C      ENDDO
CC
C10    WND = K - CENTER
C      LWORD = WND
CC
CC ... and then try to determine parameters with more precision
C      AM = INP(CENTER)
C      IC = CENTER
CC ... make sure that we stay at maximum
C      DO I = CENTER-WND+1, CENTER+WND
C        IF (INP(I) .GT. AM) THEN
C            AM = INP(I)
C            IC = I
C        ENDIF
C      ENDDO
CC ... redefine window of interest
C      Y1 = AM
C      Y  = (INP(IC+1) + INP(IC-1))/2.
C      DO I = 1, WND
C        Y2  = (INP(IC+I+1) + INP(IC-I-1))/2.
C        IF (Y1-Y .GE. Y-Y2) THEN
C            WND = I+1
C            GO TO 20
C        ENDIF
C        Y1 = Y
C        Y  = Y2
C      ENDDO
CC ... perform very simple fit
C20    A11 = 0.
C      A12 = 0.
C      A22 = WND*2.
C      B1 = 0.
C      B2 = 0.
C      DO I = IC-WND+1, IC+WND
C        AM = MAX(AM, INP(I))
C        X = I - 0.5
C        Y = INP(I) - INP(I-1)
C        A11 = A11 + X**2
C        A12 = A12 + X
C        B1  = B1  + X*Y
C        B2  = B2  + Y
C      ENDDO
C      D = (A11*A22 - A12*A12)
C      A = (B1*A22  - B2*A12)/D
C      B = (B2*A11  - B1*A12)/D
C      XC = -B/A
C      HWHI = SQRT( MAX(0., -AM/A))
C      WRITE (*,9000) AM, XC, HWHI*2.
CC
C      RETURN
C 9000 FORMAT ('Am=',F7.1,'  Xc=',F6.1,'  FWHI=',F6.3)
C      END
C
C*****************************************************************
      REAL FUNCTION SIGMA3(INP,NP)
C
C       Rough noise estimation by a "3-point dispersion" algorithm
C
      INTEGER NP
      REAL    INP(NP)
C
      INTEGER K
      REAL S
C
      S = 0.
      DO 10 K = 2, NP-1
          S = S + (INP(K) - (INP(K-1)+INP(K+1))/2.)**2
10    CONTINUE
      SIGMA3 = SQRT(S/(NP-3))
      RETURN
      END

C*****************************************************************
      SUBROUTINE BGRND(INP,NP,WRK,OUT,PORT,NOISE,NIT,CIT)
C
C       Background extraction
C
      INTEGER NP                 ! vectors lehgth
      REAL    INP(NP)            ! input vector
      REAL    WRK(NP)            ! working array
      REAL    OUT(NP)            ! output vector
      INTEGER PORT               ! smoothing window = PORT*2+1
      REAL    NOISE              ! "cutting" level
      INTEGER NIT                ! # of steps to pass
      INTEGER CIT                ! current step (IN/OUT!!!)
C
      INTEGER IT, I, II, IO, N
      REAL S
C
C ...  copy input values to the working area
      IF (CIT .EQ. 0) THEN
           DO 5 I = 1,NP
5              WRK(I) =INP(I)
      ENDIF
C
C ...  execute NIT steps of iteration: "smoothing" + "cutting"
      DO 60 IT = 1,NIT
         CIT = CIT + 1
C
C ...    first - simlple smooth with a "sliding-window"
         S = 0.
         DO 10 I = 1,PORT
10           S = S + WRK(I)           ! fill smoothing window
         N = PORT
         DO 40 I = 1,NP
             II = I + PORT
             IO = I - PORT-1
             IF (II .LE. NP) THEN
                 S = S + WRK(II)      ! add input element to window
                 N = N + 1
             ENDIF
             IF (IO .GE. 1 ) THEN
                 S = S - WRK(IO)      ! del. output elem. from window
                 N = N - 1
             ENDIF
             OUT(I) = S / N
40      CONTINUE
C
C ...   and then supperss strong details with noise level account
        DO 50 I = 1,NP
             IF (INP(I) .GT. OUT(I)+NOISE) THEN
                 WRK(I) = OUT(I)
             ELSE
                 WRK(I) = INP(I)
             ENDIF
50      CONTINUE
C
60    CONTINUE
C
      RETURN
      END

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

      SUBROUTINE PLTARR(OUT,N2,NB,NE,NAME,OVR)
C
C       Plot elements 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
         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

C *********************************************************************
C   INSERTED POLYFIT.FOR
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  15:55 - 19 NOV 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.PURPOSE
C  Evaluation of 2D-polynomial
C
C.KEYWORDS
C  REGRESSION
C
C.ALGORITHM
C
C
C-----------------------------------------------------------
C
      REAL FUNCTION EVAPOL(PAR, K, L, VAL1, VAL2)
C
      INTEGER K, L                    ! polynomial degrees
      REAL    PAR(30)                 ! arr. of (K+1)*(L+1) coef-s
      REAL    VAL1, VAL2              ! arguments
C
      INTEGER NA, K1, L1, IP
      DOUBLE PRECISION WORK(30),RESULT,DC
C
      IP     = 0
      DC     = 1.D0
      RESULT = 0.D0
      NA = (K+1)*(L+1)
      IF ( NA.GE.1 .AND. NA.LE.30) THEN
          DO L1 = 0,L
              IP     = IP + 1
              WORK(IP) = DC
              RESULT = RESULT + WORK(IP)*PAR(IP)
              IF (K .GT. 0) THEN
                  DO K1 = 1,K
                      IP = IP + 1
                      WORK(IP) = WORK(IP-1)*VAL1
                      RESULT = RESULT + WORK(IP)*PAR(IP)
                  END DO
              END IF
              DC = DC*VAL2
          ENDDO
      END IF
C
      EVAPOL = RESULT
C
      RETURN
      END

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  16:02 - 19 NOV 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.PURPOSE
C
C  2D-polynomial fitting
C
C.KEYWORDS
C
C  POLYNOMIAL FIT
C
C.ALGORITHM
C
C  fitting algorithm using Householder transformation
C
C-----------------------------------------------------------
C
      SUBROUTINE POLFIT(VAL,NROW,IDEG1,IDEG2,COEF,ERROR)
      IMPLICIT NONE
C
      INTEGER NROW                       ! # data elements:
      REAL    VAL(3,NROW)                !    { Y, X1, X2 }
      INTEGER IDEG1,IDEG2                ! polynomial degrees
      REAL    COEF(30)                   ! polyn-l coefficients (OUT)
      REAL    ERROR                      ! fitting error        (OUT)
C
      INTEGER IPAR(13)
      INTEGER I, ISTAT, KUN, KMIN
      INTEGER NVAR, K, L, N, LL, N1
      INTEGER IROW, K1
      CHARACTER*20 CPAR
      REAL RPAR(13)
C
      DOUBLE PRECISION XMIN,XMAX,YMIN,YMAX,V1,V2,V3
      DOUBLE PRECISION G(30,30),X(30)
C
      ERROR = -1.
      DO 15 I = 1, 30
          X(I) = 0.D0
   15 CONTINUE
      K      = IDEG1
      L      = IDEG2
      IF (K.EQ.0 .AND. L.EQ.0) THEN
          GO TO 70
      ELSE
          N      = (K+1) * (L+1)
      END IF
      IF (N.GT.29) THEN
          GO TO 70
      END IF
      IF (K.EQ.0 .OR. L.EQ.0) THEN
          NVAR = 2
      ELSE
          NVAR = 3
      END IF
C
C
C ... SOLVE LEAST SQUARE PROBLEM
C
      XMIN   = 1.D30
      YMIN   = 1.D30
      XMAX   = -XMIN
      YMAX   = -YMIN
      LL     = 0
      N1     = 0
      DO 50 IROW = 1,NROW
          V1     = VAL(1,IROW)
          V2     = VAL(2,IROW)
          V3     = VAL(3,IROW)
          XMIN   = MIN(XMIN,V2)
          XMAX   = MAX(XMAX,V2)
          YMIN   = MIN(YMIN,V3)
          YMAX   = MAX(YMAX,V3)
          CALL TDSET2(LL+1,V2,V3,V1,K,L,G,X,N,30)
          IF (LL.NE.0) THEN
              KMIN   = MIN(LL,N+1)
              DO 40 K1 = 1,KMIN
                  CALL TDHHTR(K1,LL+1,G,X,N,30)
   40         CONTINUE
          END IF
          LL     = MIN(LL+1,N+1)
   50 CONTINUE
      CALL TDSOLV(G,X,N,30)
      ERROR  = ABS(G(N+1,N+1))
C
C ... COPY RESULT OF REGRESSION
C
      CPAR(9:16) = 'ORDER   '
      IPAR(1) = NROW
      IPAR(2) = NVAR - 1
      IPAR(3) = 3             ! :Y
      CPAR(17:20) = 'MULT'
      IPAR(4) = 2             ! :X
      IPAR(5) = 1             ! :ORDER
      IPAR(6) = IDEG1
      IPAR(7) = IDEG2
      RPAR(1) = XMIN
      RPAR(2) = XMAX
      RPAR(3) = YMIN
      RPAR(4) = YMAX
      RPAR(5) = SQRT((ERROR*ERROR)/IPAR(1))
      CALL STKWRC('OUTPUTC',1,CPAR,1,20,KUN,ISTAT)
      CALL STKWRI('OUTPUTI',IPAR,1,13,KUN,ISTAT)
      CALL STKWRR('OUTPUTR',RPAR,1,5,KUN,ISTAT)
      CALL STKWRD('OUTPUTD',X,1,N,KUN,ISTAT)
70    DO I = 1, N
          COEF(I) = X(I)
      END DO
C
C ... END
C
      RETURN
      END
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  18:24 - 11 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.KEYWORDS           TABLE, APPLICATIONS
C.ENVIRONMENT  MIDAS
C.PURPOSE
C LINEAR REGRESSION
C.ALGORITHM
C
C HOUSEHOLDER TRANSFORMATION
C REF.: CH.L.LAWSON AND R.J.HANSON
C       SOLVING LEAST SQUARES PROBLEMS
C       PRENTICE HALL INC. 1974
C
C
C------------------------------------------------------------------
C
      SUBROUTINE TDSET2(LL,X1,X2,Y,K,L,G,X,N,M)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.PURPOSE
C
C PREPARE DATA IN WORK SPACE FOR 1 OR 2D POLYNOMIAL
C
C------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER LL
      DOUBLE PRECISION X1
      DOUBLE PRECISION X2
      DOUBLE PRECISION Y
      INTEGER K
      INTEGER L
      INTEGER N
      INTEGER M
      DOUBLE PRECISION G(M,M)
      DOUBLE PRECISION X(M)
C
      INTEGER IP,L1,K1
      DOUBLE PRECISION DC

      IP     = 0
      DC     = 1.D0
      DO 20 L1 = 0,L
          IP     = IP + 1
          G(LL,IP) = DC
          DO 10 K1 = 1,K
              IP     = IP + 1
              G(LL,IP) = G(LL,IP-1)*X1
   10     CONTINUE
          DC     = DC*X2
   20 CONTINUE
      G(LL,N+1) = Y
      RETURN

      END

      SUBROUTINE TDHHTR(K,L1,G,X,N,M)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.PURPOSE
C
C     Apply a Householder transformation to the two rows K and L1 of G.
C
C------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER K
      INTEGER L1  !   IN : RANK OF THE MATRIX
      INTEGER N
      INTEGER M
      DOUBLE PRECISION G(M,M)
      DOUBLE PRECISION X(M)
C
      INTEGER KH,J
      DOUBLE PRECISION S,H,B
C
C ... CONSTRUCT THE TRANSFORMATION
C
      S      = DSQRT(G(K,K)**2+G(L1,K)**2)
      IF (G(K,K).GE.0.D0) S      = -S
      H      = G(K,K) - S
      G(K,K) = S
      B      = G(K,K)*H
C
C ... APPLY THE TRANSFORMATION TO THE ROWS K AND L1 OF G.
C
      KH     = N + 1 - K
      IF ((KH.NE.0) .AND. (B.NE.0.D0)) THEN
          DO 10 J = 1,KH
              S      = (G(K,K+J)*H+G(L1,K+J)*G(L1,K))/B
              G(K,K+J) = G(K,K+J) + S*H
              G(L1,K+J) = G(L1,K+J) + S*G(L1,K)
   10     CONTINUE
      END IF

      RETURN

      END

      SUBROUTINE TDSOLV(G,X,N,M)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.PURPOSE
C
C
C                Solve the set of linear equations which are
C                represented by the upper triangular part
C                of the matrix G(1:n,1:n+1).
C
C------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER N
      INTEGER M
      DOUBLE PRECISION G(M,M)  
      DOUBLE PRECISION X(M)
C
      INTEGER I,II,J
      DOUBLE PRECISION S
C
      DO 10 I = 1,N
          X(I)   = G(I,N+1)
   10 CONTINUE
      DO 30 II = 1,N
          I      = N - II + 1
          S      = 0.D0
          IF (I.NE.N) THEN
              DO 20 J = 2,N + 1 - I
                  S      = S + G(I,I+J-1)*X(I+J-1)
   20         CONTINUE
          END IF

          X(I)   = (X(I)-S)/G(I,I)
   30 CONTINUE
      RETURN

      END

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C Author: P. Ballester. Creation 17-AUG-1992
C
C Find minimum and maximum in a 1D array of type REAL
C
      SUBROUTINE MINMAX(TAB,NEL,MIN,MAX)

      INTEGER NEL
      REAL    TAB(NEL),MIN,MAX

      INTEGER CNT

      MIN = TAB(1)
      MAX = TAB(1)

      IF (NEL .GT. 1) THEN

          DO 1000 CNT = 2,NEL

             IF (TAB(CNT) .GT. MAX) MAX = TAB(CNT)
             IF (TAB(CNT) .LT. MIN) MIN = TAB(CNT)

1000      CONTINUE

      ENDIF

      RETURN
      END

      

