C @(#)tdregt.for	19.1 (ESO-DMD) 02/25/03 14:11:20
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  16:14 - 19 NOV 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.IDENTIFICATION
C
C  program TDREGT
C
C.PURPOSE
C
C  POLYNOMIAL FIT OF ONE OR TWO VARIABLES BETWEEN TWO TABLES
C  THE PROGRAM CAN BE USED TO FIND CORRESPONDENCES BETWEEN
C  ENTRIES IN THE TWO TABLES
C
C  REGRESSION/TABLE TABLE X1[,X2] CATAL Y1[,Y2] DEGREE TOL [COEFFS]
C
C  THE METHOD IS NOT YET IMPLEMENTED FOR TWO VARIABLES
C
C.KEYWORDS
C
C  TABLES, REGRESSION, IDENTIFICATIONS
C
C.ALGORITHM
C
C  THERE ARE TWO MODES OF OPERATION :
C  1. INITIAL GUESS FOR THE CORRESPONDENCE IS DEFINED BY THE
C     REFERENCE COLUMNS OF BOTH TABLES
C     IN THIS CASE REFERNCE COLUMN IN SECOND OUT HAS TO BE SORTED
C  2. INITIAL GUESS IS DEFINED BY A SET OF COEFFS. WITH THE
C     GEOMETRIC TRANSFORMATION
C  THE RESULT IS
C  1. THE SET OF COEFFS. OF THE TRANSFORMATION
C     X1 <- Y1    OR
C     X1 <- Y1,Y2
C
C  2. THE CORRESPONDENCE BETWEEN THE TWO TABLES IS DEFINED BY
C     UPDATING THE REFERENCE COLUMN OF THE FIRST INPUT TABLE
C
C  3. If two columns of the problem table are specified, the
C     second column defines the iteration variable.
C
C.VERSION
C 020807	last modif
C
C-----------------------------------------------------------
C
      SUBROUTINE TDREGT
      IMPLICIT NONE
C
C ... PARAMETER DEFINITION
C
      CHARACTER    LINE1*64,KEY*64
      CHARACTER*17 COLX1,COLX2,COLY1,COLY2
      CHARACTER*8  FORM
      CHARACTER*8  FORM1,FORM2
      CHARACTER*16 LABEL1,LABEL2,UNIT1,UNIT2
      CHARACTER*64 TABLE,CATAL
      CHARACTER*80 OUTPUT1, OUTPUT2, OUTPUT3
C
      INTEGER     MADRID
      INTEGER DTYPE1, DTYPE2, DTYPE, DTYPR1, DTYPR2
      INTEGER IDEG(2),ICOL1(3),ICOL2(3),ICOL3(2)
      INTEGER NPIX(2),INDX1(400),INDX2(400)
      INTEGER IACT,STATUS,INDX,NINC,NINC1,NINC2
      INTEGER NOUTC, TID, TID1, LEN, NACT, NRCOL, IST, I, J
      INTEGER NSTART, NORD, NCOL2, NROW2, NSC2, NAC2, NAR2
      INTEGER*8  IBUF
      INTEGER NCOL1, NROW1, NSC1, NAC, NAR, NBYTES, IS
      INTEGER INDEX, KUN, KNUL
C
      REAL TOL
      REAL START(2)
      REAL STEP(2)
C
      DOUBLE PRECISION ERROR,COEFF(30),COEFFR(30)
C
       INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
       INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      DATA STEP/1.,1./
      DATA START/0.,0./
      DATA OUTPUT1
     ./'  scan         total  unid.  ident.  abs(res)    '/
      DATA OUTPUT2
     ./'  no.  iterat. lines  lines  lines   >3*std.dev.    stdev.'/
      DATA OUTPUT3
     ./'  ---- ------  -----  -----  -----   -----------    ------'/
C
C ... GET INPUT PARAMETERS + DEFAULT
C
      CALL STKRDC('P1',1,1,64,IACT,TABLE,KUN,KNUL,STATUS)
      CALL STKRDC('P2',1,1,64,IACT,LINE1,KUN,KNUL,STATUS)
      INDX   = INDEX(LINE1,',')
      IF (INDX.EQ.0) THEN
          NINC = 1
          COLX1  = LINE1
      ELSE
          NINC = 2
          COLX1  = LINE1(1:INDX-1)
          COLX2  = LINE1(INDX+1:)
      END IF
      CALL STKRDC('P3',1,1,64,IACT,CATAL,KUN,KNUL,STATUS)
      CALL STKRDC('P4',1,1,64,IACT,LINE1,KUN,KNUL,STATUS)
      INDX   = INDEX(LINE1,',')
      IF (INDX.EQ.0) THEN
          NOUTC = 1
          COLY1  = LINE1
      ELSE
          NOUTC = 2
          COLY1  = LINE1(1:INDX-1)
          COLY2  = LINE1(INDX+1:)
          CALL STTPUT(' Not yet implemented ',STATUS)
          GO TO 50
      END IF
      CALL STKRDI('INPUTI',1,2,IACT,IDEG,KUN,KNUL,STATUS)
      CALL STKRDR('INPUTR',1,1,IACT,TOL,KUN,KNUL,STATUS)
      CALL STKRDC('P7',1,1,64,IACT,KEY,KUN,KNUL,STATUS)
C
C ... INIT INPUT TABLE
C
      TID = -1
      CALL TBTOPN(TABLE,F_U_MODE,TID,STATUS)
      CALL TBIGET(TID,NCOL1,NROW1,NSC1,NAC,NAR,STATUS)
      CALL TBKGET(TID,ICOL1(1),STATUS)
      IF (ICOL1(1).LE.0) THEN
          CALL STTPUT(' There is no REFERENCE column in first table',
     +                STATUS)
          GO TO 50
      END IF
      CALL TBFGET(TID,ICOL1(1),FORM,LEN,DTYPR1,STATUS)
      CALL TBCSER(TID,COLX1,ICOL1(2),STATUS)
      CALL TBFGET(TID,ICOL1(2),FORM,LEN,DTYPE,STATUS)
      IF (DTYPE.EQ.D_C_FORMAT) THEN
          CALL STTPUT(' Illegal column type in first table',STATUS)
          GO TO 50
      END IF
      IF (NINC.EQ.2) THEN
          CALL TBCSER(TID,COLX2,ICOL1(3),STATUS)
          CALL TBFGET(TID,ICOL1(3),FORM,LEN,DTYPE,STATUS)
          IF (DTYPE.EQ.D_C_FORMAT) THEN
              CALL STTPUT(' Illegal column type in first table',STATUS)
              GO TO 50
          END IF
      END IF
      IF (ICOL1(2).EQ.-1 .OR. ICOL1(3).EQ.-1) THEN
          CALL STTPUT(' Column in first table not found',STATUS)
          GO TO 50
      END IF
C
C ... INIT REFERENCE TABLE
C
      TID1 = -1
      CALL TBTOPN(CATAL,F_I_MODE,TID1,STATUS)
      CALL TBIGET(TID1,NCOL2,NROW2,NSC2,NAC2,NAR2,STATUS)
      CALL TBKGET(TID1,ICOL2(1),STATUS)
      IF (ICOL2(1).LE.0) THEN
          CALL STTPUT(' There is no REFERENCE column in second table',
     +                STATUS)
          GO TO 50
      END IF
      IF (NSC2 .NE. ICOL2(1)) THEN
        CALL STTPUT(' Second table is not sorted by REFERENCE column',
     +                STATUS)
          GO TO 50
      END IF
      CALL TBFGET(TID1,ICOL2(1),FORM,LEN,DTYPR2,STATUS)
      IF (DTYPR1.NE.DTYPR2) THEN
          CALL STTPUT(' Warning: columns with different types',STATUS)
C          GO TO 50
      END IF
      CALL TBCSER(TID1,COLY1,ICOL2(2),STATUS)
      CALL TBFGET(TID1,ICOL2(2),FORM1,LEN,DTYPE1,STATUS)
      IF (DTYPE1.EQ.D_C_FORMAT) THEN
          CALL STTPUT(' Illegal column type in second table',STATUS)
          GO TO 50
      END IF
      IF (NOUTC.EQ.2) THEN
          CALL TBCSER(TID1,COLY2,ICOL2(3),STATUS)
          CALL TBFGET(TID1,ICOL2(3),FORM2,LEN,DTYPE2,STATUS)
          IF (DTYPE2.EQ.D_C_FORMAT) THEN
              CALL STTPUT(' Illegal column type in second table',STATUS)
              GO TO 50
          END IF
      END IF
      IF (ICOL2(2).EQ.-1 .OR. ICOL2(3).EQ.-1) THEN
          CALL STTPUT(' Column in second table not found',STATUS)
          GO TO 50
      END IF
C
C ... IDENTIFY FEATURES(iteration on equal order numbers )
C
      NINC1 = NINC + 1
      NINC2 = NOUTC + 1
      NPIX(1) = NROW1
      NPIX(2) = 6*2
C
C ... PREPARE WORK SPACE
C
      NBYTES = NPIX(1)*NPIX(2)*4
      CALL TDMGET(NBYTES,IBUF,IS)
      CALL SPCOP1(NROW1,TID,NINC1,ICOL1,TID1,NINC2,ICOL2,
     +             MADRID(IBUF),NACT,NORD,INDX1,INDX2,NSTART)
C
C ... HEADER
C
      CALL STTPUT(OUTPUT1, STATUS)
      CALL STTPUT(OUTPUT2, STATUS)
      CALL STTPUT(OUTPUT3, STATUS)
C
C ... IDENTIFICATION OF MAIN REFERENCE ORDER
C
      NRCOL  = NINC2 - 1
      ERROR  = TOL
      CALL IDENTB(NROW1,MADRID(IBUF),INDX1(NSTART),INDX2(NSTART),
     +     TID1,NROW2,NRCOL,ICOL2(2),ERROR,COEFF,IDEG,NSTART,IST)
      DO 10 I = 1,30
         COEFFR(I) = COEFF(I)
   10 CONTINUE
C
C ... IDENTIFICATION OF PREVIOUS ORDERS
C
      IF (NSTART.NE.1) THEN
          DO 20 I = NSTART - 1,1,-1
              DO 15 J = 1,30
                 COEFF(J) = COEFFR(J)
   15         CONTINUE
              CALL IDENTA(NROW1,MADRID(IBUF),INDX1(I),INDX2(I),
     +        TID1,NROW2,NRCOL,ICOL2(2),ERROR,COEFF,IDEG,I,IST)
   20     CONTINUE
      END IF
C
C ... IDENTIFICATION OF NEXT ORDERS
C
      IF (NSTART.NE.NORD) THEN
          DO 40 I = NSTART + 1,NORD
              DO 30 J = 1,30
                 COEFF(J) = COEFFR(J)
   30         CONTINUE
              CALL IDENTA(NROW1,MADRID(IBUF),INDX1(I),INDX2(I),
     +        TID1,NROW2,NRCOL,ICOL2(2),ERROR,COEFF,IDEG,I,IST)
   40     CONTINUE
      END IF
C
C ... COPY FINAL RESULTS WITH THE IDENTIFICATIONS IN OUT TABLE
C
      CALL TBCSER(TID,COLY1,ICOL3(1),STATUS)
      IF (ICOL3(1).LE.0) THEN
        CALL TBFGET(TID1,ICOL2(2),FORM1,LEN,DTYPE1,STATUS)
        CALL TBUGET(TID1,ICOL2(2),UNIT1,STATUS)
        CALL TBLGET(TID1,ICOL2(2),LABEL1,STATUS)
        CALL TBCINI(TID,DTYPE1,1,FORM1,UNIT1,LABEL1,ICOL3(1),STATUS)
      ENDIF
      IF (NOUTC.EQ.2) THEN
          CALL TBCSER(TID,COLY2,ICOL3(2),STATUS)
          IF (ICOL3(1).LE.0) CALL TBCINI(TID,DTYPE2,1,FORM2,UNIT2,
     +                                   LABEL2,ICOL3(2),STATUS)
      END IF
      CALL SPCOP2(NROW1,MADRID(IBUF),TID,NOUTC,ICOL3)
C
C ... end
C
   50 CALL TBTCLO(TID,STATUS)
      CALL TBTCLO(TID1,STATUS)
      CALL TDMFRE(NBYTES,IBUF,IS)
      RETURN
      END

      SUBROUTINE IDENTA(NROW1,BUF,I1,I2,TID,NROW2,NC,IC,ERROR,
     +                  COEFF,IDEG,INDX,IST)
C
C IDENTIFY DETECTED LINES
C INITIAL GUESSES DEFINED BY THE COEFFS
C       INPUT STORED IN BUF WITH THE FOLLOWING STRUCTURE
C BUF(I,1) - X1 (SAMPLE)
C BUF(I,2) - X2 (ORDER)
C BUF(I,3) - Y1 (IDENTIFIED WAVELENGTH)
C BUF(I,4) - Y2 (SECOND VARIABLE)
C BUF(I,5) - LINE WEIGHT AS    1 MANUAL IDENTIFICATION
C                              0 NOT IDENTIFIED
C                             -1 AUTOMATIC IDENTIFICATION
C BUF(I,6) - RESIDUALS
C
      IMPLICIT NONE
C
      INTEGER NROW1, I1, I2, NROW2, NC, IST, NSTEP, NV, NTOT
      INTEGER NTOTAL, NTOTU, NTOTI, NTOTR1,NTOTR2, I, TID, ISTAT
      INTEGER IC(2),IDEG(2),INDX
C
      DOUBLE PRECISION ERROR,BUF(NROW1,6),COEFF(30)
      DOUBLE PRECISION STDEV,STDEV1,STDEV3
C
      CHARACTER*80 OUT
C
      STDEV1 = -1.D20
      IST    = 0
C
C ITERATION
C
      DO 20 NSTEP = 1,10
          NV     = (IDEG(1)+1)* (IDEG(2)+1)
C
C   BUILD L - S MATRIX AND SOLVE SYSTEM
C
          IF (NSTEP.GT.1) CALL LSOLVE(NROW1,BUF,I1,I2,IDEG,NV,COEFF)
C
C   IDENTIFY FEATURES
C
          CALL LSIDEN(IDEG,NV,COEFF,NROW1,BUF,I1,I2,TID,NROW2,NC,
     +                IC,ERROR,NTOT,STDEV)
C
C   SELECT LINES WITH LOW RESIDUALS
C
          STDEV3 = STDEV*3.D0
          IF (NTOT.EQ.0) THEN
              OUT = ' '
              CALL STTPUT(' Wrong identifications in input table',ISTAT)
              IST    = 1
              RETURN
          END IF
C
C STATISTICS ON THE IDENTIFICATIONS
C  !   TOTAL NUMBER OF LINES
          NTOTAL = 0  !   TOTAL NUMBER OF UNIDENTIFIED LINES
          NTOTU  = 0  !   TOTAL NUMBER OF IDENTIFIED LINES
          NTOTI  = 0  !   TOTAL NUMBER OF LINES WITH RESIDUALS .LT. 0.1 PIXEL
          NTOTR1 = 0  !   TOTAL NUMBER OF LINES WITH RESIDUALS .GT.3*STDEV
          NTOTR2 = 0
          DO 10 I = I1,I2
              NTOTAL = NTOTAL + 1
              IF (DABS(BUF(I,5)).GT.0.5D0) THEN
                  NTOTI  = NTOTI + 1
                  IF (DABS(BUF(I,6)).GT.STDEV3) THEN
                      BUF(I,4) = 0.0D0
                      BUF(I,5) = 0.0D0
                      BUF(I,6) = 0.0D0
                      NTOTR2   = NTOTR2 + 1
                  END IF
              END IF
   10     CONTINUE
          NTOTU  = NTOTAL - NTOTI
          IF (NTOTI.LT.NV) THEN
              CALL STTPUT(' Error : not enough  identified entries',
     +                    ISTAT)
              IST    = 1
              RETURN
          END IF
          IF (DABS((STDEV-STDEV1)/STDEV).LE.0.005D0) GOTO 100
          STDEV1 = STDEV
   20 CONTINUE
      RETURN
100   WRITE(OUT,300) INDX,NSTEP,NTOTAL,NTOTU,NTOTI,NTOTR2,STDEV
      CALL STTPUT(OUT,ISTAT)
      RETURN
300   FORMAT(1X,I5,I6,2X,I5,2X,I5,2X,I6,2X,I8,5X,E10.3)
      END

      SUBROUTINE IDENTB(NROW1,BUF,I1,I2,TID,NROW2,NC,IC,ERROR,
     +                  COEFF,IDEG,INDX,IST)
C
C IDENTIFY DETECTED LINES
C INITIAL GUESSES DEFINED BY THE IDENTIFICATIONS IN BUF(I,3)
C       INPUT STORED IN BUF WITH THE FOLLOWING STRUCTURE
C BUF(I,1) - X1 (SAMPLE)
C BUF(I,2) - X2 (ORDER)
C BUF(I,3) - Y1 (IDENTIFIED WAVELENGTH)
C BUF(I,4) - Y2 (SECOND VARIABLE)
C BUF(I,5) - LINE WEIGHT AS    1 MANUAL IDENTIFICATION
C                              0 NOT IDENTIFIED
C                             -1 AUTOMATIC IDENTIFICATION
C BUF(I,6) - RESIDUALS
C
      IMPLICIT NONE
C
      INTEGER NROW1, I1, I2, NROW2, NC, IST, NSTEP, NV, NTOT
      INTEGER NTOTAL, NTOTU, NTOTI, NTOTR1,NTOTR2, I, TID
      INTEGER ISTAT,INDX
      INTEGER IC(2),IDEG(2)
C
      DOUBLE PRECISION ERROR,BUF(NROW1,6),COEFF(30)
      DOUBLE PRECISION STDEV,STDEV1,STDEV3
C
      CHARACTER*80 OUT
C
      STDEV1 = -1.D20
      IST    = 0
C
C ITERATION
C
      DO 20 NSTEP = 1,10
          NV     = (IDEG(1)+1)* (IDEG(2)+1)
C
C   BUILD L - S MATRIX AND SOLVE SYSTEM
C
          CALL LSOLVE(NROW1,BUF,I1,I2,IDEG,NV,COEFF)
C
C   IDENTIFY FEATURES
C
          CALL LSIDEN(IDEG,NV,COEFF,NROW1,BUF,I1,I2,TID,NROW2,NC,
     +                 IC,ERROR,NTOT,STDEV)
C
C   SELECT LINES WITH LOW RESIDUALS
C
          STDEV3 = STDEV*3.D0
          IF (NTOT.EQ.0) THEN
              OUT = ' '
              CALL STTPUT(' Wrong identifications in input table',ISTAT)
              IST    = 1
              RETURN
          END IF
C
C STATISTICS ON THE IDENTIFICATIONS
C  !   TOTAL NUMBER OF LINES
          NTOTAL = 0  !   TOTAL NUMBER OF UNIDENTIFIED LINES
          NTOTU  = 0  !   TOTAL NUMBER OF IDENTIFIED LINES
          NTOTI  = 0  !   TOTAL NUMBER OF LINES WITH RESIDUALS .LT. 0.1 PIXEL
          NTOTR1 = 0  !   TOTAL NUMBER OF LINES WITH RESIDUALS .GT.3*STDEV
          NTOTR2 = 0
          DO 10 I = I1,I2
              NTOTAL = NTOTAL + 1
              IF (DABS(BUF(I,5)).GT.0.5D0) THEN
                  NTOTI  = NTOTI + 1
                  IF (DABS(BUF(I,6)).GT.STDEV3) THEN
                      BUF(I,4) = 0.0D0
                      BUF(I,5) = 0.0D0
                      BUF(I,6) = 0.0D0
                      NTOTR2   = NTOTR2 + 1
                  END IF
              END IF
   10     CONTINUE
          NTOTU  = NTOTAL - NTOTI
          IF (NTOTI.LT.NV) THEN
              CALL STTPUT(' Error : not enough  identified entries',
     +                    ISTAT)
              IST    = 1
              RETURN
          END IF
          IF (DABS((STDEV-STDEV1)/STDEV).LE.0.005D0) GOTO 100
          STDEV1 = STDEV
   20 CONTINUE
      RETURN
100   WRITE(OUT,300) INDX,NSTEP,NTOTAL,NTOTU,NTOTI,NTOTR2,STDEV
      CALL STTPUT(OUT,ISTAT)
      RETURN
300   FORMAT(1X,I5,I6,2X,I5,2X,I5,2X,I6,2X,I8,5X,E10.3)
      END

      SUBROUTINE LSIDEN(IDEG,NV,COEFF,NROW1,BUF,I1,I2,TID,NROW2,
     +                   NC,IC,ERROR,NTOT,STDEV)
C
C IDENTIFY DETECTED LINES USING DISPERSION COEFFS
C
      IMPLICIT NONE
C
      INTEGER NTOT, NV, NROW1, I1, I2, TID, NROW2, NC
      INTEGER K, L, NACT, I, IPOINT, STATUS, ISTART
      INTEGER IC(NC),IDEG(2)
C
      LOGICAL NULL(2)
C
      DOUBLE PRECISION COEFF(NV),BUF(NROW1,6),ERROR,SAMPLE
      DOUBLE PRECISION DD(2),STDEV, POLEV
C
      DATA DD/0.D0,0.D0/
      DATA NULL/.FALSE.,.FALSE./
C
C LOOP ON THE LINES IN THE CATAL
C
      NTOT   = 0
      STDEV  = 0.D0
C      K      = IDEG(1) + 1
C      L      = IDEG(2) + 1
      K      = IDEG(1)
      L      = IDEG(2)
      NACT   = I2 - I1 + 1
      ISTART = 1
      DO 10 I = 1,NROW2
C
C ... USE ONLY SUITABLE RANGE !!
C
          CALL TBRRDD(TID,I,NC,IC,DD,NULL,STATUS)
          IF ( .NOT. NULL(1) .AND. .NOT. NULL(2)) THEN
              SAMPLE = POLEV(K,L,NV,COEFF,DD)
C
              CALL SEABAD(BUF(I1,1),NACT,SAMPLE,ERROR,ISTART,IPOINT)
C
C ... LINE IDENTIFIED
C
              IF (IPOINT.GT.0) THEN
                  IPOINT = IPOINT + I1 - 1
                  BUF(IPOINT,3) = DD(1)
                  BUF(IPOINT,4) = DD(1)
                  BUF(IPOINT,5) = -1.D0
                  BUF(IPOINT,6) = BUF(IPOINT,1) - SAMPLE
                  STDEV  = STDEV + BUF(IPOINT,6)*BUF(IPOINT,6)
                  NTOT   = NTOT + 1
              END IF
          END IF
   10 CONTINUE
      IF (NTOT.GT.0) STDEV  = DSQRT(STDEV/NTOT)
      RETURN
      END
C
      SUBROUTINE LSOLVE(NROW1,BUF,I1,I2,IDEG,NV,COEFF)
C
C SOLVE THE L-S FIT TO THE DATA
      IMPLICIT NONE
C
      INTEGER NROW1, I1, I2, NV, K, L, NTOTAL, LL, N, I
      INTEGER KMIN, N1, K1, ISTAT, LL1
      INTEGER IDEG(2)
C
      DOUBLE PRECISION BUF(NROW1,6),COEFF(NV)
      DOUBLE PRECISION G(20,20)
C
C BUILD INPUT ARRAYS
C
      K      = IDEG(1)
      L      = IDEG(2)
      NTOTAL = 0
      LL     = 0
      N      = NV
      N1     = N + 1
      DO 20 I = I1,I2
        IF (DABS(BUF(I,5)).GT.0.5D0) THEN
          NTOTAL = NTOTAL + 1
          LL1    = LL+1
          CALL TDSET2(LL1,BUF(I,3),BUF(I,4),BUF(I,1),K,L,G,COEFF,N,20)
            IF (LL.NE.0) THEN
              KMIN   = MIN(LL,N1)
               DO 10 K1 = 1,KMIN
                  CALL TDHHTR(K1,LL1,G,COEFF,N,20)
   10          CONTINUE
            END IF
          LL     = MIN(LL1,N1)
        END IF
   20 CONTINUE
C
C SOLVE THE LINEAR SYSTEM
C
      IF (NTOTAL.LT.NV) THEN
          CALL STTPUT(' Not enough identified features',ISTAT)
          RETURN
      END IF
      CALL TDSOLV(G,COEFF,N,20)
      RETURN
      END

      DOUBLE PRECISION FUNCTION POLEV(K,L,NV,DPAR,VAL)
C
C COMPUTE FUNCTION VALUE
C
      IMPLICIT NONE
C
      INTEGER K, L, NV, IP, L1, K1
C
      DOUBLE PRECISION DPAR(NV),VAL(2)
      DOUBLE PRECISION DC,RESULT,WORK(30)
C
      IP     = 0
      DC     = 1.D0
      RESULT = 0.D0
      DO 20 L1 = 0, L
          IP     = IP + 1
          WORK(IP) = DC
          RESULT = RESULT + WORK(IP)*DPAR(IP)
          DO 10 K1 = 1,K
              IP     = IP + 1
              WORK(IP) = WORK(IP-1)*VAL(1)
              RESULT   = RESULT + WORK(IP)*DPAR(IP)
   10     CONTINUE
          DC     = DC*VAL(2)
   20 CONTINUE
      POLEV  = RESULT
      RETURN
      END

      SUBROUTINE SPCOP1(N,TID,NC1,IC1,TID1,NC2,IC2,BUF,NACT,
     +                   NORD,INDX1,INDX2,NSTART)
C
C PREPARE WORK SPACE
C
      IMPLICIT NONE
C
      INTEGER NACT, TID, TID1, N, NC1, NC2, NORD
      INTEGER NSTART, NC, LEN, STATUS, NIDENT, NMAX, I, I1
      INTEGER IC1(NC1),IC2(NC2)
      INTEGER INDX1(1),INDX2(1)
      INTEGER DTYPE, NEXT
C
      LOGICAL SEL,NULL1(3),NULL2(3), NULL3
C
      REAL RVAL1,RVAL2,RERR
C
      DOUBLE PRECISION BUF(N,6)
      DOUBLE PRECISION DERR
      DOUBLE PRECISION VAL1(3)
      DOUBLE PRECISION VAL2(3)
      DOUBLE PRECISION REF
C
      CHARACTER FORM*10
C
       INCLUDE 'MID_INCLUDE:TABLES.INC'
       INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      DATA VAL1/0.D0,0.D0,0.D0/
      DATA VAL2/0.D0,0.D0,0.D0/
      DATA REF/1.E30/
      DATA NULL1/.FALSE.,.FALSE.,.FALSE./
      DATA NULL2/.FALSE.,.FALSE.,.FALSE./
C
      I1     = 1
      NACT   = 0
      DERR   = .1D0
      NC     = NC2 - 1
      NORD   = 0
C
C ... CHECK FORMAT IN CATAL
C
      CALL TBFGET(TID1,IC2(1),FORM,LEN,DTYPE,STATUS)
      IF (DTYPE.EQ.D_R4_FORMAT) RERR   = DERR
C
      NSTART = 1
      NIDENT = 0
      NMAX   = 0
      DO 10 I = 1,N
          CALL TBSGET(TID,I,SEL,STATUS)
          IF (SEL) THEN
              CALL TBRRDD(TID,I,NC1,IC1,VAL1,NULL1,STATUS)
              IF ( .NOT. NULL1(2) .AND. .NOT. NULL1(3)) THEN
                  NACT   = NACT + 1  !   ASSIGN X1
                  BUF(NACT,1) = VAL1(2)  !   ASSIGN X2 (OPTIONAL)
                  BUF(NACT,2) = VAL1(3)  !   NEW ORDER
                  IF (VAL1(3).NE.REF) THEN
                      REF    = VAL1(3)
                      NORD   = NORD + 1
                      INDX1(NORD) = NACT  !   CHECK NO. OF IDENTIFICATIONS
                      IF (NIDENT.GT.NMAX) THEN
                          NMAX   = NIDENT
                          NSTART = NORD
                      END IF
                      NIDENT = 0
                      IF (NORD.GT.1) INDX2(NORD-1) = NACT - 1
                  END IF  !   ASSIGN Y1
                  BUF(NACT,3) = 0.D0  !   ASSIGN Y2 (OPTIONAL)
                  BUF(NACT,4) = 0.D0  !   FLAG AS NON-EXISTING
                  BUF(NACT,5) = 0.D0  !   RESIDUAL (X1-FIT(Y1,Y2))
                  BUF(NACT,6) = 0.D0
                  IF ( .NOT. NULL1(1)) THEN
                      NIDENT = NIDENT + 1
                      IF (DTYPE.EQ.D_R4_FORMAT) THEN
                          RVAL1  = VAL1(1)
                          CALL TBESRR(TID1,IC2,RVAL1,RERR,I1,NEXT,
     +                                STATUS)
                          CALL TBERDR(TID1,NEXT,IC2,RVAL2,
     +                                       NULL3,STATUS)
                          VAL2(3) = RVAL2
                      ELSE
                          CALL TBESRD(TID1,IC2,VAL1(1),DERR,I1,NEXT,
     +                                STATUS)
                          CALL TBERDD(TID1,NEXT,IC2,VAL2(3),NULL3,
     +                                STATUS)
                      END IF
                      IF (NEXT.GT.0) THEN
                          CALL TBSGET(TID1,NEXT,SEL,STATUS)
                          IF (SEL) THEN
                              CALL TBRRDD(TID1,NEXT,NC,IC2(2),
     .                        VAL2,NULL2,STATUS)
                              IF (.NOT.NULL2(1).AND..NOT.NULL2(2)) THEN !ASSIGN Y1
                                  BUF(NACT,3) = VAL2(1)  !   ASSIGN Y2 (OPTIONAL)
                                  BUF(NACT,4) = VAL2(2)  !   FLAG AS EXISTING
                                  BUF(NACT,5) = 1.D0
                                  I1 = NEXT              ! ASSUMES SORTED TABLES
                              END IF
                          END IF
                      END IF
                  END IF
              END IF
          END IF
   10 CONTINUE
      IF (NIDENT.GT.NMAX) NSTART = NORD
      NSTART = NSTART-1
      INDX2(NORD) = NACT
      RETURN
      END

      SUBROUTINE SPCOP2(N,BUF,TID,NC,IC)
C
C COPY FINAL IDENTIFICATIONS
C
      IMPLICIT NONE
C
      INTEGER N, TID, NC, I, NACT, STATUS
      INTEGER IC(2)
C
      LOGICAL SEL
C
      DOUBLE PRECISION BUF(N,6),VAL(3)
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
       INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      NACT   = 0
      DO 10 I = 1,N
          CALL TBSGET(TID,I,SEL,STATUS)
          IF (SEL) THEN
              NACT   = NACT + 1
              IF (DABS(BUF(NACT,5)).GT.0.5D0) THEN
                  VAL(1) = BUF(NACT,3)
                  VAL(2) = BUF(NACT,4)
                  CALL TBRWRD(TID,I,NC,IC,VAL,STATUS)
              ELSE
                  CALL TBEDEL(TID,I,IC(1),STATUS)
                  IF (NC.EQ.2) CALL TBEDEL(TID,I,IC(2),STATUS)
              END IF
          END IF
   10 CONTINUE
      RETURN
      END

       SUBROUTINE SEABAD(ARRAY, N, VALUE, ERROR, IFIRST,IPOS)
C
C       BINARY SEARCH OF ELEMENTS IN THE ARRAY
C       DOUBLE  PRECISION VERSION
C       ASCENDING ORDER
C
       IMPLICIT NONE
       INTEGER   N
       DOUBLE PRECISION ARRAY(N)
       DOUBLE PRECISION VALUE
       DOUBLE PRECISION ERROR
       INTEGER   IFIRST
       INTEGER   IPOS
C
       INTEGER I, LOWER, UPPER
C
C ... INITIALIZE COUNTERS
C
       LOWER = IFIRST
       UPPER = N
       IPOS  = 0
C
C ... LOOP
C
5      CONTINUE
         I = (LOWER + UPPER)/2
C
C ... COMPARE KEYS
C
         IF (DABS(VALUE-ARRAY(I)).LE.ERROR) THEN
           IPOS = I
           GO TO 10
         ELSE
           IF (VALUE.LT.ARRAY(I)) THEN
             UPPER = I - 1
           ELSE
             LOWER = I + 1
           ENDIF
         ENDIF
       IF (LOWER.LE.UPPER) GO TO 5
C
C ... FIND FIRST ENTRY IN A SET OF EQUAL KEYS
C
10     IF (IPOS.LE.IFIRST) RETURN
       IF (DABS(VALUE-ARRAY(IPOS-1)).GT.ERROR) RETURN
20     CONTINUE
         IPOS = IPOS - 1
         IF (IPOS.EQ.IFIRST) RETURN
       IF (DABS(VALUE-ARRAY(IPOS-1)).LE.ERROR) GOTO 20
       RETURN
       END

