C @(#)neciden.for	17.1.1.1 (ESO-IPG) 01/25/02 17:51:29
C @(#)neciden.for	13.2 (ES0-DMD) 03/22/99 10:25:15
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  18:08 - 6 DEC 1987
C          1.1  read LINTAB to get the name of the line table
C               instead using fixed line.tbl. (SW)
C          1.2  Store the mean RMS in AVRMS descriptor (SW,1999-10-01)
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: P. BALLESTER
C
C
C.IDENTIFICATION
C
C  program ECHIDENEW.FOR
C
C.PURPOSE
C
C  Execute the command
C  CALIBRATE/ECHELLE CALIB-table error [degree]
C
C.KEYWORDS
C
C  echelle, line identification
C
C.ALGORITHM
C
C  SOURCE:
C  Adaptation of ECHIDEN2.FOR to involve the "echelle relation" in the
C  order by order polynomial fitting.
C  INPUT:
C  user- table contains detected lines (LINE.tbl)
C  Wavelengths have been assigned by the procedure ECHIDENEW
C  Coefficients of echelle relation are written in descriptor REGRD in LINE.tbl
C  Degree of the echelle relation fitting is written in keyword ECHI(3)
C  ALGORITHM:
C   This sample is used to find a high order approximation 
C  to the dispersion coeffs.
C
C
C.INPUT/OUTPUT
C
C  P1-P2 contain input parameters
C
C-----------------------------------------------------------
C
      PROGRAM ECHIDNEW
      IMPLICIT         NONE
C
C ... parameter definition
C
      INTEGER          STATUS
      INTEGER          IC(2),TID
      INTEGER          MADRID
      INTEGER          IAV, IDEG, IORD, ISTAT, KUN, KNUL
      INTEGER          NCOL1, NROW1, NSC1, NACOL, NAROW
      INTEGER          ICOLX, ICOLW, ICOLO, IC1, IR, IS
      INTEGER*8        IM1, IX, IDNT, K1, IWAV1, IRES
      INTEGER          NN1
      INTEGER          ISS, NACT, NO, IAC
      INTEGER*8        IBUF
      INTEGER          ISEL
      INTEGER*8        SELECT
      CHARACTER*60     TABLE
      CHARACTER*80     LINE
      CHARACTER*16     MSG
      REAL             W(2),ERROR

C******************************************************************************
C**   COMMON BLOCK. Don't forget to update all copies of this block if modified
      INTEGER          DEGMAX,ORDMAX
      PARAMETER        (DEGMAX=7,ORDMAX=500) ! Don't change 7
C     DEGMAX: Max. number of poly. coefficients. ORDMAX: Max. number of orders
      DOUBLE PRECISION A(DEGMAX,ORDMAX)   ! Array of order by order coefficients
      DOUBLE PRECISION A1(0:DEGMAX,0:DEGMAX)! Global 2D solution
C Initial coefficients (global echelle relation)
      INTEGER          IA(ORDMAX)  ! Number of identified values in each order
      INTEGER          ABSORD(ORDMAX) ! Absolute order number
      COMMON   /POLY/  A,A1,IA,ABSORD
C******************************************************************************
C
      DOUBLE PRECISION     AINIT(DEGMAX*DEGMAX)
      DOUBLE PRECISION     AVRMS(2)
      INTEGER              LOOP1,LOOP2,IPOS,MODEG

      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      COMMON/VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      DATA MSG/'ERR:ECHIDENxxxx'/
C
C ... get into MIDAS
C
      CALL STSPRO('ECHIDENEW')
C
C ... Checks pre-defined arrays dimensions
C
      CALL STKRDI('INPUTI',1,1,IAV,IDEG,KUN,KNUL,ISTAT)  ! Degree of polynomial
      CALL STKRDI('INPUTI',2,1,IAV,IORD,KUN,KNUL,ISTAT)  ! Number of orders
      CALL STKRDI('INPUTI',4,1,IAV,MODEG,KUN,KNUL,ISTAT) ! Number of orders
      CALL STKRDR('INPUTR',1,1,IAV,ERROR,KUN,KNUL,ISTAT) ! Allowed residual
      IF (IDEG.GE.DEGMAX) THEN
          WRITE (LINE,9020) IDEG,DEGMAX
          CALL  STETER(9,LINE)
      ENDIF
      IF (IORD.GT.ORDMAX) THEN
          WRITE (LINE,9020) IORD,ORDMAX
          CALL  STETER(9,LINE)
      ENDIF
C
C ... Displays input summary
C
      TABLE  = 'line'
      CALL STKRDC('LINTAB',1,1,60,IAC,TABLE,KUN,KNUL,ISTAT)
      CALL STTPUT(' COMPUTE DISPERSION COEFFICIENTS',ISTAT)
      CALL STTPUT(' -------------------------------',ISTAT)
      CALL STTPUT(' INPUT TABLE : '//TABLE,ISTAT)
      WRITE (LINE,9000) IDEG
      CALL STTPUT(LINE,ISTAT)
C
C ... Init input table
C
      CALL TBTOPN(TABLE,F_U_MODE,TID,ISTAT)
C
      CALL STDRDD(TID,'REGRD',1,50,IAC,AINIT,KUN,KNUL,ISTAT)
      DO 100 LOOP2 = 0,IDEG
      DO 100 LOOP1 = 0,IDEG
         IPOS = LOOP1*(IDEG+1) + LOOP2 + 1
         A1(LOOP1,LOOP2) = AINIT(IPOS)
100   CONTINUE
C
      CALL TBIGET(TID,NCOL1,NROW1,NSC1,NACOL,NAROW,ISTAT)
      CALL TBCSER(TID,':X ',ICOLX,ISTAT)
      CALL TBCSER(TID,':IDENT',ICOLW,ISTAT)
      CALL TBCSER(TID,':ORDER ',ICOLO,ISTAT)
      CALL TBCSER(TID,':WAVEC ',IC1,ISTAT)
      CALL TBCSER(TID,':RESIDUAL ',IR,ISTAT)
      CALL TBCSER(TID,':SELECT ',ISEL,ISTAT)
      IF (ICOLX.EQ.-1 .OR. ICOLW.EQ.-1 .OR. ICOLO.EQ.-1) THEN
          CALL STTPUT(' Column not present ',IS)
          GO TO 10
      END IF
C
      CALL TBCMAP(TID,0,IM1,ISTAT)         ! Select. col.  MADRID(IM1)
      IF (ISTAT.NE.0) GO TO 10
      CALL TBCMAP(TID,ICOLX,IX,ISTAT)      ! :X            MADRID(IX)
      IF (ISTAT.NE.0) GO TO 10
      CALL TBCMAP(TID,ICOLW,IDNT,ISTAT)    ! :IDENT        MADRID(IDNT)
      IF (ISTAT.NE.0) GO TO 10
      CALL TBCMAP(TID,ICOLO,K1,ISTAT)      ! :ORDER        MADRID(K1)    
      IF (ISTAT.NE.0) GO TO 10
      CALL TBCMAP(TID,IC1,IWAV1,ISTAT)     ! :WAVEC        MADRID(IWAV1)
      IF (ISTAT.NE.0) GO TO 10
      CALL TBCMAP(TID,IR,IRES,ISTAT)       ! :RESIDUAL     MADRID(IRES)
      IF (ISTAT.NE.0) GO TO 10
      CALL TBCMAP(TID,ISEL,SELECT,ISTAT)   ! :SELECT       MADRID(SELECT)
      IF (ISTAT.NE.0) GO TO 10
C
      NN1 = 6*8*NROW1
      CALL TDMGET(NN1,IBUF,ISS)            ! WORK SPACE    MADRID(IBUF)

C
C ... Loads values in work space
C
      CALL SPCOP3(NROW1,MADRID(IM1),MADRID(K1),MADRID(IX),
     +             MADRID(IDNT),MADRID(IBUF),NACT)
C
C ... Computes coeffs for each order
C
      CALL CCOEF (NROW1,MADRID(IBUF),NACT,IDEG,NO,W(1),W(2),ERROR,MODEG
     $     ,AVRMS(1),ISTAT)

      AVRMS(2) = ISTAT          !save number of bad RMS (2D-fit)
C
C ... Improves solution for undefined orders
C
C      CALL IMPCOE (IORD,IDEG)
C
C ... Computes residuals for each wavelength and unselect those which
C      are bigger than ERROR.
C
      CALL CRES(NROW1,MADRID(IM1),MADRID(K1),MADRID(IX),MADRID(IDNT),
     +    MADRID(IWAV1),MADRID(IRES),IDEG,NO,NROW1,MADRID(IBUF),NACT,
     +    MADRID(SELECT),ERROR)

C
C ... Write results
C
      IC(1)  = IDEG
      IC(2)  = NO
      CALL STKWRR('OUTPUTR',W,1,2,KUN,ISTAT)
      CALL STDWRI(TID,'COEFS',IC,1,2,KUN,ISTAT)
      CALL STDWRI(TID,'COEFI',IA,1,NO,KUN,ISTAT)
      NO     = NO*7
      CALL STDWRD(TID,'COEFD',A,1,NO,KUN,ISTAT)
      CALL STDWRD(TID,'AVRMS',AVRMS,1,2,KUN,ISTAT)
C
C ... End
C
      CALL TBTCLO(TID,ISTAT)
   10 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9010) ISTAT
          CALL TDERRR(ISTAT,MSG,STATUS)
      END IF
      CALL STSEPI
 9000 FORMAT (' POLYNOMIAL DEGREE : ',I3)
 9010 FORMAT (I4)
 9020 FORMAT (' OVERFLOW IN MODULE ECHIDENEW.FOR - VALUE : ',I5,
     +        '>> MAXI : ',I5)
      END

C******************************************************************************
      SUBROUTINE SPCOP3(N,M,Y,X,W,BUF,NACT)
C******************************************************************************
C
C PREPARE WORK SPACE
C
C******************************************************************************

      IMPLICIT NONE
      INTEGER          N, I
      INTEGER          TINULL, NACT, M(N), INTNB
      REAL             Y(N), X(N), TRNULL, YMIN, YMAX, YCURR
      DOUBLE PRECISION TDNULL, BUF(N,6), W(N)
C
      CALL TBMNUL(TINULL,TRNULL,TDNULL)
C 
C New storage procedure
C
      NACT   = 0
      YMIN   = MIN(Y(1),Y(N))
      YMAX   = MAX(Y(1),Y(N))

      DO 20 YCURR  =  YMAX,YMIN,-1
      INTNB  =  0

      DO 10 I = 1, N
          IF (M(I).GT.0.and.Y(I).EQ.YCURR) THEN ! If the row is selected then
              INTNB       = 1
              NACT        = NACT + 1
              BUF(NACT,1) = Y(I)           ! Load ORDER number
              BUF(NACT,2) = X(I)           ! Load X position
              IF (W(I).EQ.TDNULL) THEN
                  BUF(NACT,3) = -1.D0      ! Load Wavelength
              ELSE
                  BUF(NACT,3) = 0.0001*W(I)
              END IF
              BUF(NACT,4) = I              ! Original position
              BUF(NACT,5) = 1.D0           ! Selection flag
              BUF(NACT,6) = 0.D0           ! Residual
          END IF
   10 CONTINUE

      IF (INTNB .EQ. 0) THEN
         NACT = NACT + 1
         BUF(NACT,1)  =  YCURR
         BUF(NACT,2)  =  0.D0
         BUF(NACT,3)  =  0.D0
         BUF(NACT,4)  =  0.D0
         BUF(NACT,5)  =  0.D0           ! Selection flag
         BUF(NACT,6)  =  0.D0           ! Residual
      ENDIF

   20 CONTINUE

      RETURN
      END
C******************************************************************************
      SUBROUTINE CCOEF(NROW1,BUF,NACT,IDEG,NO,WST,WEN,ERROR,MODEG,AVRMS
     $     ,IBAD)
C******************************************************************************
C
C IDENTIFY DETECTED LINES
C       INPUT STORED IN BUF WITH THE FOLLOWING STRUCTURE
C BUF(I,1) - ORDER NUMBER
C BUF(I,2) - SAMPLE VALUE
C BUF(I,3) - WAVELENGTH
C
C******************************************************************************

      IMPLICIT NONE
      INTEGER          NROW1, IDEG, NO, ISEQ, IBAD, IK, IDEG1
      INTEGER          I, IP, IK1, ND, ISTAT, NACT, NB, MODEG
      REAL             W1, W2, WST, WEN, ERROR
      DOUBLE PRECISION BUF(NROW1,6),TOL,RMS,AVRMS
      CHARACTER*80     LINE
      CHARACTER*80     LINE1
      CHARACTER*80     LINE2
      CHARACTER*80     LINE3,LINE4
C******************************************************************************
C**   COMMON BLOCK. Don't forget to update all copies of this block if modified
      INTEGER          DEGMAX,ORDMAX
      PARAMETER        (DEGMAX=7,ORDMAX=500) ! Don't change 7
C     DEGMAX: Max. number of poly. coefficients. ORDMAX: Max. number of orders
      DOUBLE PRECISION A(DEGMAX,ORDMAX)   ! Array of order by order coefficients
      DOUBLE PRECISION A1(0:DEGMAX,0:DEGMAX)
C Initial coefficients (global echelle relation)
      INTEGER          IA(ORDMAX)  ! Number of identified values in each order
      INTEGER          ABSORD(ORDMAX) ! Absolute order number
      COMMON   /POLY/  A,A1,IA,ABSORD
C******************************************************************************
      INTEGER          IFIRST(ORDMAX), INDATA(ORDMAX), NBLINE
C
      DATA LINE/
     +      ' SEQ.NO  SPECTRAL  NO.LINES  WL START    WL END  STD. DEV.'
     +             /
      DATA LINE1/
     +      '          ORDER                                  ANGSTROEM'
     +             /
      DATA LINE2/
     +      ' ------  --------  --------  --------    ------  ---------'
     +             /
      DATA LINE3/
     +      ' ---------------------------------------------------------'
     +             /
C
      AVRMS  = 0.0
      TOL    = ERROR/10000.D0
      NBLINE = 0
      NO     = 0
      ISEQ   = 0
      IBAD   = 0
      WST    = 0.
      WEN    = 0.
      IDEG1  = IDEG + 1
      IK     = -999
      CALL STTPUT(LINE,ISTAT)
      CALL STTPUT(LINE1,ISTAT)
      CALL STTPUT(LINE2,ISTAT)
      DO 10 I = 1,NACT
          IK1    = BUF(I,1)
          IF (IK.NE.IK1) THEN
              IK     = IK1
              IF (NO.NE.0) THEN
                  IP     = IFIRST(NO)
                  ND     = INDATA(NO)
                  ISEQ   = ISEQ + 1
                  CALL FITPOL(ND,BUF(IP,2),BUF(IP,3),IDEG1,A(1,NO),
     +                 IA(NO),BUF(IP,1),ISEQ,W1,W2,BUF(IP,5),NB,
     +                 TOL,BUF(IP,6),MODEG,RMS)
                  ABSORD(NO) = BUF(IP,1)
                  NBLINE = NBLINE + NB
                  IF (WST.EQ.0)    WST    = W1
                  WEN    = W2
                  IF (RMS .LT. 99.99) THEN
                     AVRMS = AVRMS + RMS
                  ELSE
                     IBAD = IBAD + 1
                  ENDIF
              END IF
              NO     = NO + 1
              IFIRST(NO) = I
              INDATA(NO) = 1
          ELSE
              INDATA(NO) = INDATA(NO) + 1
          END IF
   10 CONTINUE
      ISEQ   = ISEQ + 1
      IP     = IFIRST(NO)
      ND     = INDATA(NO)
      CALL FITPOL(ND,BUF(IP,2),BUF(IP,3),IDEG1,A(1,NO),IA(NO),
     +            BUF(IP,1),ISEQ,W1,W2,BUF(IP,5),NB,TOL,
     +            BUF(IP,6),MODEG,RMS)
      IF (IBAD.NE.ISEQ) THEN
         AVRMS = (AVRMS + RMS)/(ISEQ-IBAD)
      ELSE
         WRITE (LINE4,'(A)') 'BAD RESULT !'
         AVRMS = 99.9999
      ENDIF
      ABSORD(NO) = BUF(IP,1)
      NBLINE = NBLINE + NB
      IF (WST.EQ.0)    WST    = W1
      WEN    = W2
      CALL STTPUT(LINE3,ISTAT)
      WRITE (LINE4,'(39X,A,F9.5)') 'MEAN RMS: ',AVRMS
      CALL STTPUT(LINE4,ISTAT)
      WRITE  (LINE4, 1000) NBLINE
1000  FORMAT ('** TOTAL NUMBER OF LINES : ',I6,' **')
      CALL STTPUT(LINE4,ISTAT)
      RETURN
      END

C******************************************************************************
      SUBROUTINE CRES(N,M,K,X,W,WT,RES,IDEG,NO,NROW1,BUF,NACT,
     +                SELECT,ERROR)
C******************************************************************************
C
C COMPUTE RESIDUALS
C
C******************************************************************************

      IMPLICIT NONE
      INTEGER          N, IDEG, NO, I
      INTEGER          TINULL, M(N),SELECT(N)
      DOUBLE PRECISION TDNULL
      DOUBLE PRECISION W(N), WT(N), RES(N)
      REAL             K(N),X(N)
      REAL             TRNULL, ERROR
C
      INTEGER          NROW1,NACT,POS
      DOUBLE PRECISION BUF(NROW1,6)
C
C******************************************************************************
C**   COMMON BLOCK. Don't forget to update all copies of this block if modified
      INTEGER          DEGMAX,ORDMAX 
      PARAMETER        (DEGMAX=7,ORDMAX=500)  ! Don't change 7
C     DEGMAX: Max. number of poly. coefficients. ORDMAX: Max. number of orders
      DOUBLE PRECISION A(DEGMAX,ORDMAX)   ! Array of order by order coefficients
      DOUBLE PRECISION A1(0:DEGMAX,0:DEGMAX)
C Initial coefficients (global echelle relation)
      INTEGER          IA(ORDMAX)  ! Number of identified values in each order
      INTEGER          ABSORD(ORDMAX) ! Absolute order number
      COMMON   /POLY/  A,A1,IA,ABSORD
C******************************************************************************
C
      CALL TBMNUL(TINULL,TRNULL,TDNULL)
C
      DO 15 I=1,N
         SELECT(I) = 0
   15 CONTINUE
C
      DO 20  I = 1,NACT

         POS = NINT(BUF(I,4))

         RES(POS) = BUF(I,6)

         IF (M(POS).GT.0.AND.BUF(I,5).GT.0.D0.AND.
     +       ABS(BUF(I,6)).LT.ERROR)  THEN
             SELECT(POS) = 1
C            if select(i) = 1 or 2, the value is unchanged.
         ENDIF

   20 CONTINUE

      RETURN
      END

C******************************************************************************
      SUBROUTINE FITPOL(N,X,W,IDEG1,COEFF,NOMB,KK,ISEQ,WST,
     +                  WEN,SEL,NUMBER,TOL,RES,MODEG,REF)
C******************************************************************************
C
C COMPUTE COEFFS
C
C******************************************************************************

      IMPLICIT NONE
      INTEGER          N, IDEG1, ISEQ, I, J, NR, NUMBER, ORDER, OUT
      INTEGER          ISTAT,NP,IAV,KUN,KNUL,NOMB,MODEG
      REAL             WST, WEN
      DOUBLE PRECISION X(N), W(N), KK, SEL(N)
      DOUBLE PRECISION TOL,REF
      DOUBLE PRECISION XX, POLEV, WAVE1, WAVE2, STDERR
      CHARACTER*80     LINE
C******************************************************************************
C**   COMMON BLOCK. Don't forget to update all copies of this block if modified
      INTEGER          DEGMAX,ORDMAX
      PARAMETER        (DEGMAX=7,ORDMAX=500) ! Don't change 7
C     DEGMAX: Max. number of poly. coefficients. ORDMAX: Max. number of orders
      DOUBLE PRECISION A(DEGMAX,ORDMAX)   ! Array of order by order coefficients
      DOUBLE PRECISION A1(0:DEGMAX,0:DEGMAX)
C Initial coefficients (global echelle relation)
      INTEGER          IA(ORDMAX)  ! Number of identified values in each order
      INTEGER          ABSORD(ORDMAX) ! Absolute order number
      COMMON   /POLY/  A,A1,IA,ABSORD
C******************************************************************************
      DOUBLE PRECISION COEFF(DEGMAX),X1(ORDMAX),W1(ORDMAX)
      INTEGER          SEL1(ORDMAX)
      DOUBLE PRECISION START, STEP, RES(ORDMAX), IX_IR8
C
      NR     = 0
      ORDER  = KK
      CALL STKRDI('INPUTI',3,1,IAV,NP,KUN,KNUL,ISTAT) ! Number of col.
      CALL STKRDD('INPUTD',1,1,IAV,START,KUN,KNUL,ISTAT) 
      CALL STKRDD('INPUTD',2,1,IAV,STEP,KUN,KNUL,ISTAT) 

      DO 10 I = 1,N
          IF (W(I).GT.0.D0) THEN
              NR       = NR + 1
              X1(NR)   = X(I)
              W1(NR)   = W(I)
              SEL1(NR) = I
          END IF
   10 CONTINUE 

      NOMB   = MIN(NR,IDEG1)
      NUMBER = NR
      IF (NR.GT.IDEG1.AND.MODEG.EQ.1) THEN
          CALL LSOLVE(NR,X1,W1,COEFF,IDEG1,REF)
          OUT = 1
   50     CONTINUE
          IF (NUMBER.GT.(IDEG1+1).AND.OUT.NE.0) THEN

              REF=STDERR(NUMBER,X1,W1,DEGMAX,COEFF,OUT,TOL,RES)
              IF (OUT.NE.0) THEN
                  SEL(SEL1(OUT)) = 0
              ENDIF
              NUMBER = 0
              DO 60 I = 1,N
                 IF (W(I).GT.0.D0.AND.SEL(I).GT.0.5D0) THEN
                     NUMBER     = NUMBER + 1
                     X1(NUMBER) = X(I)
                     W1(NUMBER) = W(I)
                     SEL1(NUMBER) = I
                 END IF
   60         CONTINUE 
              CALL LSOLVE (NUMBER,X1,W1,COEFF,IDEG1,REF)
              GOTO 50
          ENDIF

          REF    = REF*10000.D0
          XX     = IX_IR8(1,START,STEP)
          WAVE1  = POLEV(DEGMAX,COEFF,XX)
          WAVE1  = 10000.D0*WAVE1
          XX     = IX_IR8(NP,START,STEP)
          WAVE2  = POLEV(DEGMAX,COEFF,XX)
          WAVE2  = 10000.D0*WAVE2
          WRITE (LINE,9000) ISEQ,ORDER,NUMBER,WAVE1,WAVE2,REF
          CALL STTPUT(LINE,ISTAT)
          WST    = WAVE1
          WEN    = WAVE2
      ELSE
          DO 25  I = 1,IDEG1
              COEFF(I)   = DBLE(0.)
              DO 20 J=1,IDEG1
              COEFF(I)   = COEFF(I) + A1(J-1,I-1)*(FLOAT(ORDER)**
     +                     FLOAT(J-2))
20        CONTINUE
          COEFF(I) = COEFF(I)/DBLE(10000.)
25        CONTINUE
          XX     = IX_IR8(1,START,STEP)
          WAVE1  = POLEV(DEGMAX,COEFF,XX)
          WAVE1  = 10000.D0*WAVE1
          XX     = IX_IR8(NP,START,STEP)
          WAVE2  = POLEV(DEGMAX,COEFF,XX)
          WAVE2  = 10000.D0*WAVE2
          IF (NR.GT.0) THEN
              REF = STDERR(NR,X1,W1,DEGMAX,COEFF,OUT,TOL,RES)
              OUT    = 0
              REF    = REF*10000.D0
          ELSE
              REF    = 99.99999
          ENDIF

          IF (MODEG.EQ.1) THEN
            WRITE (LINE,9010) ISEQ,ORDER,NR,WAVE1,WAVE2,REF
          ELSE
            WRITE (LINE,9020) ISEQ,ORDER,NR,WAVE1,WAVE2,REF
          ENDIF

          CALL STTPUT(LINE,ISTAT)
          WST    = WAVE1
          WEN    = WAVE2
      END IF
      RETURN
 9000 FORMAT (1X,I6,2X,I8,2X,I8,2X,F8.2,2X,F8.2,2X,F9.5)
 9010 FORMAT (1X,I6,2X,I8,2X,I8,2X,F8.2,2X,F8.2,2X,F9.5,
     +                                           ' *NOT ENOUGH LINES*')
 9020 FORMAT (1X,I6,2X,I8,2X,I8,2X,F8.2,2X,F8.2,2X,F9.5,
     +                                           ' *FROM 2D SOLUTION*')
      END

C******************************************************************************
      DOUBLE PRECISION FUNCTION POLEV(DEGMAX,A,XX)
C******************************************************************************
C
C COMPUTE FUNCTION VALUE
C
C******************************************************************************

      IMPLICIT         NONE
      INTEGER          DEGMAX, I
      DOUBLE PRECISION A(1),XX,YY
C
      YY = DBLE(0.)
      DO 10 I = DEGMAX, 1, -1
         YY = A(I) + XX * YY
10    CONTINUE
C
      POLEV  = YY
      RETURN
      END

C******************************************************************************
      SUBROUTINE LSOLVE(NR,X1,W1,A,IA,REF)
C******************************************************************************
C
C L - S POLYNOMIAL
C
C******************************************************************************

      IMPLICIT NONE
C
      INTEGER NR,IA
      DOUBLE PRECISION X1(NR),W1(NR),A(IA),REF
      DOUBLE PRECISION G(20,20),X(20)
      INTEGER N,L,N1,I,J,M,KMIN,K
      COMMON /LS/G,X,N
C
      L      = 0
      N      = IA
      N1     = N + 1
      DO 30 I = 1,NR
          J      = L + 1
          G(J,1) = 1.D0
          DO 10 M = 2,N
              G(J,M) = G(J,M-1)*X1(I)
   10     CONTINUE
          G(J,N1) = W1(I)
          IF (L.NE.0) THEN
              KMIN   = MIN(L,N1)
              DO 20 K = 1,KMIN
                  CALL HT(K,J)
   20         CONTINUE
          END IF

          L      = MIN(J,N1)
   30 CONTINUE
      CALL SOLVE
      DO 40 I = 1,IA
          A(I)   = X(I)
   40 CONTINUE
C      REF    = DABS(G(N1,N1)) ! Result non normalized. 19.12.90 PB
      REF    = DSQRT(G(N1,N1)*G(N1,N1)/NR)
      RETURN
      END

C******************************************************************************
      SUBROUTINE HT(K,L1)
C******************************************************************************
C     FUNCTION : Apply a Householder transformation to the
C                two rows K and L1 of G.
C
C     INPUT :    * matrix G(1:L1,1:n+1)
C                * numbers L1 and K of the two rows
C                * n : rank of the matrix
C
C     OUTPUT :   modified rows K and L1 of the matrix G.
C******************************************************************************

      IMPLICIT NONE
      INTEGER  K,L1
      DOUBLE PRECISION G,X,S,H,B
      INTEGER N,KH,J
      COMMON /LS/G(20,20),X(20),N
C
C     Construct the transformation
C
      S      = DSQRT(G(K,K)**2+G(L1,K)**2)
      IF (G(K,K).GE.0) THEN
          S      = -S
      END IF
      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)) 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
      END

C******************************************************************************
      SUBROUTINE SOLVE
C******************************************************************************
C     FUNCTION : 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     INPUT :    * matrix G(1:n,1:n+1)
C                * n : rank of the matrix
C
C     OUTPUT :   solution vector x(1:n)
C******************************************************************************

      IMPLICIT NONE
      DOUBLE PRECISION G,X,S
      INTEGER N,I,II,J
      COMMON /LS/G(20,20),X(20),N
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

      END
C****************************************************************
C
      DOUBLE PRECISION FUNCTION STDERR (NVAL, X, Y, DIM, 
     +                   COEFF, OUTLIE, TOL,RES)
C
      IMPLICIT NONE
C
      INTEGER          NVAL,DIM
      DOUBLE PRECISION X(1), Y(1), RES(1), COEFF(DIM), STDMAX
      DOUBLE PRECISION YESTIM, SX2, RESIDU, XX, POLEV
      INTEGER          I,OUTLIE,FLAG
      DOUBLE PRECISION TOL
C
      FLAG = 0
      SX2  = DBLE(0.)
C
      DO 10 I = 1, NVAL

           XX = X(I)
           YESTIM = POLEV(DIM,COEFF,XX)
           RES(I) = (Y(I)-YESTIM)*10000.D0
           RESIDU = ABS(Y(I)-YESTIM)

           SX2    = SX2 + RESIDU * RESIDU

           IF (RESIDU.GT.TOL) FLAG = 1

           IF (I.EQ.1) THEN
               STDMAX = RESIDU
               OUTLIE = 1
           ELSE IF (RESIDU.GT.STDMAX) THEN
               STDMAX = RESIDU
               OUTLIE = I
           ENDIF

10    CONTINUE

      IF (FLAG.EQ.0) OUTLIE = 0
      STDERR = DSQRT ( SX2 / NVAL )

      RETURN
      END

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

      SUBROUTINE IMPCOE (NORD,DEG)

      IMPLICIT NONE

      INTEGER  NORD,DEG

      INTEGER    ORDER,POSINF,POSSUP,STAT,DEGREE,NCOE,DEG1
      CHARACTER  LINE*80

C******************************************************************************
C**   COMMON BLOCK. Don't forget to update all copies of this block if modified
      INTEGER          DEGMAX,ORDMAX
      PARAMETER        (DEGMAX=7,ORDMAX=500) ! Don't change 7
C     DEGMAX: Max. number of poly. coefficients. ORDMAX: Max. number of orders
      DOUBLE PRECISION A(DEGMAX,ORDMAX)   ! Array of order by order coefficients
      DOUBLE PRECISION A1(0:DEGMAX,0:DEGMAX)
C Initial coefficients (global echelle relation)
      INTEGER          IA(ORDMAX)  ! Number of identified values in each order
      INTEGER          ABSORD(ORDMAX) ! Absolute order number
      COMMON   /POLY/  A,A1,IA,ABSORD
C******************************************************************************


      DEG1 = DEG + 1

      DO 10 ORDER = 1,NORD

         IF (IA(ORDER).LT.DEG1) THEN

C ...        Looking for the inferior and superior correctly estimated order
             POSINF = ORDER
             POSSUP = ORDER
20           POSINF = POSINF - 1
             IF (POSINF.GT.0.AND.IA(POSINF).LT.DEG1) GOTO 20
30           POSSUP = POSSUP + 1
             IF (POSSUP.LE.NORD.AND.IA(POSSUP).LT.DEG1) GOTO 30

C ...        Computes coefficients from POSINF and from POSSUP (if possible)
             IF (POSINF.GT.0.OR.POSSUP.LE.NORD) THEN

                 NCOE = 0

                 DO 40 DEGREE = 1,DEG1
                       A(DEGREE,ORDER) = 0.D0
40               CONTINUE

                 IF (POSINF.GT.0) THEN
                     NCOE = NCOE + 1
                     DO 50 DEGREE = 1,DEG1
              A(DEGREE,ORDER) = A(DEGREE,POSINF)*ABSORD(POSINF)
     +                                               /ABSORD(ORDER)
50                   CONTINUE
                 ENDIF

                 IF (POSSUP.LE.NORD) THEN
                     NCOE = NCOE + 1
                     DO 60 DEGREE = 1,DEG1
                       A(DEGREE,ORDER) = A(DEGREE,ORDER) +
     +       A(DEGREE,POSSUP)*ABSORD(POSSUP)/ABSORD(ORDER)
60                   CONTINUE
                 ENDIF

                 IF (NCOE.GT.1) THEN
                     DO 70 DEGREE = 1,DEG1
                        A(DEGREE,ORDER) = A(DEGREE,ORDER)/NCOE
70                   CONTINUE
                 ENDIF

                 WRITE (LINE,1000) ABSORD(ORDER)
1000             FORMAT ('Replaced coefficients order : ',I6)
                 CALL   STTPUT(LINE,STAT)

           ENDIF
        ENDIF
10    CONTINUE

      RETURN
      END