C @(#)caliline.for	17.1.1.1 (ES0-DMD) 01/25/02 17:56:13
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  16:39 - 3 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.IDENTIFICATION
C
C  CALILINE.FOR
C
C.PURPOSE
C
C  WAVELENGTH CALIBRATION. THE PROGRAM COMPUTES DISPERSION COEFFS,
C  ESTIMATED WAVELENGTHS AND RESIDUALS BY MEANS OF A REGRESSION ANALYSIS
C
C  DISPERSION COEFFICIENTS ARE REALLY THE POLYNOMIAL COEFFS OF THE
C  TRANSFORMATION
C  X <- F(WAVE)
C
C.KEYWORDS
C
C  TABLES, 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 OUTPUT HAS TO BE SORTED
C  2. INITIAL GUESS IS DEFINED BY A SET OF COEFFS. WITH THE
C     GEOMETRIC TRANSFORMATION
C
C.INPUT/OUTPUT
C
C  CALIBRATE/LINE error,degree tabline CATALue
C
C  error     - window defining the tolerance in pixels
C  degree    - degree of the polynomial approximation
C  tabline   - table with the positions of the lines. It contains the
C              following information:
C              X        line position in pixels
C              Y        scan line
C              PEAK     maximum value
C              IDENT    initial identification
C              WAVE     actual  identification
C              WAVEC    computed values
C              RESIDUAL residuals WAVE-WAVEC
C              DELTA    relative residuals (a0+a1*X - WAVEC)
C              DELTA1   relative residuals (a0+a1*X - WAVE)
C  CATALue - line CATALue
C
C-----------------------------------------------------------
C
C
C ... PARAMETER DEFINITION
C
      IMPLICIT NONE
C
      CHARACTER*17 COLY
      CHARACTER*16 LABEL(7),UNIT
      CHARACTER*8  FORM
      CHARACTER*80 TABLE,CATAL,LINE
C
      INTEGER   MADRID(1),NAC2,NAR2,NAC1,NAR1
      INTEGER   IDEG,ICOL1(7),ICOL2(2)
      INTEGER   NPIX1(2),NPIX2(2),STATUS,IACT,I
      INTEGER   NCOL1,NROW1,NSC1,NCOL2,NROW2,NSC2
      INTEGER*8 IBUF1,IBUF2
      INTEGER   IST,NACT1,NACT2,NIDENT
      INTEGER   NINCL1,NINCL2,TID,TIC
      INTEGER   NWORK1, NWORK2, KUN, KNUL
C
      INTEGER SIZE,MODE
      PARAMETER (SIZE=30)
      DOUBLE PRECISION ERROR,COEFF(SIZE)
C
      REAL TOL(2)
      REAL START(2)
      REAL STEP(2)
C
      INCLUDE     'MID_INCLUDE:TABLES.INC'
      COMMON/VMR/MADRID
      INCLUDE     'MID_INCLUDE:TABLED.INC'
C
      DATA UNIT/'ANGSTROM'/
      DATA FORM/'F15.7'/
      DATA START/0.,0./
      DATA STEP/1.,1./
      DATA LABEL(1)/'X'/,LABEL(2)/'IDENT'/,LABEL(3)/'WAVE'/,
     +     LABEL(4)/'WAVEC'/,LABEL(5)/'RESIDUAL'/,
     +     LABEL(6)/'DELTA'/,LABEL(7)/'DELTAC'/
      DATA COLY/':WAVE'/
C
C ... GET INTO MIDAS
C
      CALL STSPRO('CALILINE')
C
C ... GET INPUT PARAMETERS + DEFAULT
C
      CALL STKRDC('P2',1,1,64,IACT,TABLE,KUN,KNUL,STATUS)
      CALL STKRDC('P3',1,1,64,IACT,CATAL,KUN,KNUL,STATUS)
      CALL STKRDR('INPUTR',1,2,IACT,TOL,KUN,KNUL,STATUS)
      ERROR  = TOL(1)
      IDEG   = TOL(2)
C
C ... INIT INPUT TABLE
C
      CALL TBTOPN(TABLE,F_U_MODE,TID,STATUS)
      CALL TBIGET(TID,NCOL1,NROW1,NSC1,NAC1,NAR1,STATUS)

      DO 15 I = 1,2
         CALL TBCSER(TID,LABEL(I),ICOL1(I),STATUS)
         IF (ICOL1(I).LT.0) THEN
             WRITE (LINE,17) LABEL(I)
17           FORMAT ('Could not find column label ',A20)
             CALL STETER (9,LINE)
         ENDIF
15    CONTINUE

      DO 10 I = 3,7
         CALL TBCSER(TID,LABEL(I),ICOL1(I),STATUS)
         IF (ICOL1(I).LT.0) THEN
            CALL TBCINI(TID,D_R8_FORMAT,1,FORM,UNIT,LABEL(I),
     +               ICOL1(I),STATUS)
         ENDIF
   10 CONTINUE
C
C ... INIT REFERENCE TABLE
C
      CALL TBTOPN(CATAL,0,TIC,STATUS)
      CALL TBIGET(TIC,NCOL2,NROW2,NSC2,NAC2,NAR2,STATUS)
      CALL TBCSER(TIC,COLY,ICOL2(1),STATUS)
C
C ... PREPARE WORK SPACE
C
      NINCL1 = 2
      NPIX1(1) = NROW1
      NPIX1(2) = 7*2
      NWORK1   = 7*8*NROW1
      CALL TDMGET(NWORK1,IBUF1,STATUS)
      CALL SPCOP1(TID,NROW1,NINCL1,ICOL1,MADRID(IBUF1),NACT1)
      NINCL2 = 1
      NPIX2(1) = NROW2
      NPIX2(2) = 2
      NWORK2   = 8*NROW2
      CALL TDMGET(NWORK2,IBUF2,STATUS)
      CALL SPCOP3(TIC,NROW2,NINCL2,ICOL2,MADRID(IBUF2),NACT2)
C
C ... VERIFY IDENTIFICATIONS
C
      CALL GUESS(MODE,COEFF,SIZE,IDEG)

      IF (MODE.EQ.0)  THEN
      CALL VERIFY(NROW1,NACT1,MADRID(IBUF1),NACT2,
     +            MADRID(IBUF2),NIDENT)
      IF (NIDENT.LE.IDEG) THEN
          CALL STTPUT(' Not enough identifications',STATUS)
          GO TO 20
      END IF
      ENDIF
C
C ... IDENTIFICATION
C
      CALL IDENTF(NROW1,NACT1,MADRID(IBUF1),NROW2,NACT2,
     +            MADRID(IBUF2),ERROR,IDEG,COEFF,IST,MODE)
      IF (IST.NE.0) THEN
          CALL STTPUT(' Wrong identifications in input table',
     +                STATUS)
          GO TO 20
      END IF
C
C ... COMPUTE WAVELENGTHS AND RESIDUALS
C
      CALL COMPRS(NROW1,NACT1,MADRID(IBUF1),IDEG)
C
C ... COPY FINAL RESULTS WITH THE IDENTIFICATIONS IN OUTPUT TABLE
C
      NINCL1 = 5
      CALL SPCOP2(NROW1,MADRID(IBUF1),TID,NINCL1,ICOL1(3))
C
C ... SAVE COEFFICIENTS
C
      CALL STKWRD('OUTPUTD',COEFF,1,20,KUN,STATUS)
C
C ... END
C
20    CONTINUE
      CALL TDMFRE(NWORK1,IBUF1,STATUS)
      CALL TDMFRE(NWORK2,IBUF2,STATUS)
      CALL TBTCLO(TID,STATUS)
      CALL TBTCLO(TIC,STATUS)
      CALL STSEPI
      STOP 
      END

      SUBROUTINE LSIDEN(IDEG,NV,COEFF,NROW1,NACT1,BUF1,NROW2,BUF2,
     +                   ERROR,NTOT,STDEV)
C
C IDENTIFY DETECTED LINES USING DISPERSION COEFFS
C
      IMPLICIT NONE
C
      INTEGER NV,NROW1,NACT1,NROW2
      DOUBLE PRECISION COEFF(NV),BUF1(NROW1,5),ERROR,SAMPLE
      DOUBLE PRECISION BUF2(NROW2)
      DOUBLE PRECISION POLEV,DD(2),STDEV
      DOUBLE PRECISION VALUE
      INTEGER IDEG,NTOT,K,L,I,IPOINT
C
      DATA DD/0.D0,0.D0/
C
C ... LOOP ON THE LINES IN THE CATALUE
C
      NTOT   = 0
      STDEV  = 0.D0
      K      = IDEG + 1
      L      = 1
      DO 10 I = 1,NROW2
C
C ... USE ONLY SUITABLE RANGE !!
C
          DD(1)  = BUF2(I)
          SAMPLE = POLEV(K,L,NV,COEFF,DD)
C
C ... CHECK LIMITS
C
          CALL TZSSD1(BUF1(1,1),NACT1,SAMPLE,ERROR,IPOINT,VALUE)
C
C ... LINE IDENTIFIED
C
          IF (IPOINT.GT.0) THEN
              BUF1(IPOINT,3) = DD(1)
              BUF1(IPOINT,2) = -1.D0
              BUF1(IPOINT,5) = VALUE - SAMPLE
              STDEV  = STDEV + BUF1(IPOINT,5)*BUF1(IPOINT,5)
              NTOT   = NTOT + 1
          END IF
   10 CONTINUE
      IF (NTOT.GT.0) STDEV  = DSQRT(STDEV/NTOT)
      RETURN
      END

      SUBROUTINE LSOLVE(NROW1,NACT1,BUF,IDEG,NV,COEFF)
C
C SOLVE THE L-S FIT TO THE DATA
C
      IMPLICIT NONE
C
      INTEGER NROW1,NACT1,IDEG,NV
      DOUBLE PRECISION BUF(NROW1,6),COEFF(NV)
      DOUBLE PRECISION G(20,20)
C
      INTEGER K,L,NTOTAL,LL,N,N1,I,KMIN,K1,ISTAT
C
C BUILD INPUT ARRAYS
C
      K      = IDEG
      L      = 0
      NTOTAL = 0
      LL     = 0
      N      = NV
      N1     = N + 1
      DO 20 I = 1,NACT1
          IF (DABS(BUF(I,2)).GT.0.5D0) THEN
              NTOTAL = NTOTAL + 1
              CALL TDSET2(LL+1,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,LL+1,G,COEFF,N,20)
   10             CONTINUE
              END IF
              LL     = MIN(LL+1,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
      INTEGER  K,L,NV
      DOUBLE PRECISION DPAR(NV),VAL(2)
      DOUBLE PRECISION DC,RESULT,WORK(30)
C
      INTEGER IP,L1,K1
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(TABLE,N,NC1,IC1,BUF,NACT)
C
C PREPARE WORK SPACE
C
      IMPLICIT NONE
C
      INTEGER TABLE,N,NC1,NACT
      DOUBLE PRECISION BUF(N,NC1),VAL1(2)
C
      INTEGER IC1(NC1),STATUS,I
      LOGICAL SEL
      LOGICAL NULL1(2)
      DATA VAL1/0.D0,0.D0/
      DATA NULL1/.FALSE.,.FALSE./
C
      NACT   = 0
      VAL1(2) = 0.D0
C
C ... READ COLUMNS
C
      DO 10 I = 1,N
C
C ...  READ SELECTION FLAG
C
          CALL TBSGET(TABLE,I,SEL,STATUS)
          IF (SEL) THEN
C
C ...    READ ROW
C
              CALL TBRRDD(TABLE,I,NC1,IC1,VAL1,NULL1,STATUS)
C     IF (.NOT.NULL1(1)) THEN
              NACT   = NACT + 1  !   ASSIGN SAMPLE
              BUF(NACT,1) = VAL1(1)  !   ASSIGN IDENTIFICATION
              BUF(NACT,2) = VAL1(2)  !   CLEAR  WAVE
              BUF(NACT,3) = 0.D0  !   CLEAR  WAVEC
              BUF(NACT,4) = 0.D0  !   CLEAR FLAG
              BUF(NACT,5) = 0.D0  !   CLEAR RESIDUAL
              BUF(NACT,6) = 0.D0
              VAL1(1) = 0.D0
              VAL1(2) = 0.D0
C     ENDIF
          END IF
   10 CONTINUE
      RETURN
      END

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

      SUBROUTINE SPCOP3(TABLE,N,NC1,IC1,BUF,NACT)
C
C PREPARE WORK SPACE
C
      IMPLICIT NONE
C
      INTEGER N,NC1,NACT
      INTEGER TABLE
      DOUBLE PRECISION BUF(N,NC1),VAL1(2)
      INTEGER IC1(NC1),STATUS,I
      LOGICAL SEL
      LOGICAL NULL1(2)
      DATA NULL1/.FALSE.,.FALSE./
      DATA VAL1/0.D0,0.D0/
C
      NACT   = 0
      VAL1(2) = 0.D0
C
C ... READ COLUMNS
C
      DO 10 I = 1,N
C
C ...  READ SELECTION FLAG
C
          CALL TBSGET(TABLE,I,SEL,STATUS)
          IF (SEL) THEN
C
C ...    READ ROW
C
              CALL TBRRDD(TABLE,I,NC1,IC1,VAL1,NULL1,STATUS)
C     IF (.NOT.NULL1(1)) THEN
              NACT   = NACT + 1  !   ASSIGN SAMPLE
              BUF(NACT,1) = VAL1(1)
              VAL1(1) = 0.D0
C     ENDIF
          END IF
   10 CONTINUE
      RETURN
      END

      SUBROUTINE VERIFY(NROW1,NACT1,BUF1,NACT2,BUF2,NIDENT)
C
C VERIFY IDENTIFICATIONS, ASSIGN CLOSEST IDENTIFICATION
C IN THE COMPARISON TABLE
C
C INPUT ARGUMENTS
C NROW1 INTG*4 DIMENSION OF BUF1
C NACT1 INTG*4 NUMBER OF ENTRIES IN BUF1
C BUF1 REAL*8 ARRAY WITH LINES
C NACT2 INTG*4 NUMBER OF ENTRIES IN BUF2
C BUF2 REAL*8 ARRAY WITH REFERENCE LINES
C
C OUTPUT ARGUENTS
C NIDENT INTG*4 NUMBER OF IDENTIFICATIONS
C
      IMPLICIT NONE
C
      INTEGER NROW1,NACT1,NACT2,NIDENT,STAT
      DOUBLE PRECISION BUF1(NROW1,3),BUF2(NACT2),VAL
      REAL   REF,DIF, RATIO, NXTLINE
      CHARACTER*80   LINE
      INTEGER I,J
C
      NIDENT = 0
C
C ... ITERATION ON ROWS
C
      DO 20 I = 1,NACT1
          VAL    = BUF1(I,2)

          IF (VAL.GT.0) THEN
              REF    = VAL
              NIDENT = NIDENT + 1
C
C ... ASSIGN USER PROVIDED VALUE
C     AND FLAG THE ENTRY
C
              BUF1(I,3) = VAL  !   FLAG ENTRY
              BUF1(I,2) = 1.D0
C
C ... CHECK THAT THERE IS AN IDENTIFICATION
C
              DO 10 J = 1,NACT2
                  DIF    = ABS(VAL-BUF2(J))
                  IF (DIF.LT.REF) THEN
                      REF     = DIF
                      NXTLINE = BUF2(J)
                  END IF
   10         CONTINUE

              RATIO = ABS(NXTLINE/VAL - 1.)
              IF (RATIO.GT.1.E-05) THEN
                  WRITE(LINE,610) VAL,NXTLINE
610               FORMAT('Warning: Line ',F12.6,' not found in catalog.'
     +        ,' Closest line:',F12.6)
                  CALL STTPUT(LINE,STAT)
              ENDIF

          END IF
   20 CONTINUE

      RETURN
      END

      SUBROUTINE COMPRS(NROW1,NACT1,BUF,IDEG)
C
C COMPUTE ESTIMATED WAVELENGTHS AND RESIDUALS
C
      IMPLICIT NONE
C
      INTEGER NROW1,NACT1,IDEG
      DOUBLE PRECISION BUF(NROW1,7),COEFF(20)
      DOUBLE PRECISION G(20,20)
C
      INTEGER K,L,LL,NTOTAL,N,N1,NV,I,K1,KMIN,ISTAT
      DOUBLE PRECISION DD(2),LINEAR,POLEV
C
C BUILD INPUT ARRAYS
C
      K      = IDEG
      L      = 0
      NTOTAL = 0
      LL     = 0
      NV     = K + 1
      N      = NV
      N1     = N + 1
      DO 20 I = 1,NACT1
          IF (DABS(BUF(I,2)).GT.0.5D0) THEN
              NTOTAL = NTOTAL + 1
              CALL TDSET2(LL+1,BUF(I,1),BUF(I,1),BUF(I,3),
     +             K,L,G,COEFF,N,20)
              IF (LL.NE.0) THEN
                  KMIN   = MIN(LL,N1)
                  DO 10 K1 = 1,KMIN
                      CALL TDHHTR(K1,LL+1,G,COEFF,N,20)
   10             CONTINUE
              END IF
              LL     = MIN(LL+1,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)
C
C ... COMPUTE ESTIMATED VALUES AND RESIDUALS
C
      DD(2)  = 0.D0
      DO 30 I = 1,NACT1
          DD(1)  = BUF(I,1)
          BUF(I,4) = POLEV(K,L,NV,COEFF,DD)
          BUF(I,5) = BUF(I,3) - BUF(I,4)
   30 CONTINUE
C
C ... REPEAT FOR LINEAR TERMS
C
      K      = 1
      L      = 0
      NTOTAL = 0
      LL     = 0
      NV     = K + 1
      N      = NV
      N1     = N + 1
      DO 50 I = 1,NACT1
          NTOTAL = NTOTAL + 1
          CALL TDSET2(LL+1,BUF(I,1),BUF(I,1),BUF(I,4),
     +                K,L,G,COEFF,N,20)
          IF (LL.NE.0) THEN
              KMIN   = MIN(LL,N1)
              DO 40 K1 = 1,KMIN
                  CALL TDHHTR(K1,LL+1,G,COEFF,N,20)
   40         CONTINUE
          END IF
          LL     = MIN(LL+1,N1)
   50 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)
C
C ... COMPUTE ESTIMATED VALUES AND RESIDUALS
C
      DD(2)  = 0.D0
      DO 60 I = 1,NACT1
          DD(1)  = BUF(I,1)
          LINEAR = POLEV(K,L,NV,COEFF,DD)
          BUF(I,6) = LINEAR - BUF(I,4)
          IF (DABS(BUF(I,2)).GT.0.5D0) BUF(I,7) = LINEAR - BUF(I,3)
   60 CONTINUE
      RETURN
      END

        SUBROUTINE TZSSD1(ARRAY, N, VALUE, ERROR, IPOS, ACTVAL)
C
C        SEQUENTIAL SEARCH FOR ENTRY CLOSE TO VALUE
C        DOUBLE PRECISION VERSION
C
        IMPLICIT NONE
        INTEGER                 N
        DOUBLE PRECISION        ARRAY(N)
        DOUBLE PRECISION        VALUE, ERROR, ACTVAL
        INTEGER                 IPOS
C
        INTEGER                 I, IFIRST, TINULL
        DOUBLE PRECISION        DIST, DIST1, TDNULL
        REAL                    TRNULL
C
C ... INITIALIZE FLAG
C
        IPOS  = 0
        DIST1 = 1.D30
        CALL TBMNUL(TINULL,TRNULL,TDNULL)
C
C ... LOOP ON VALUES
C
        IFIRST = 1
        DO 10 I = IFIRST, N
          IF (ARRAY(I).NE.TDNULL) THEN
            DIST=DABS(ARRAY(I)-VALUE)
            IF (DIST.LT.DIST1.AND.DIST.LE.ERROR) THEN
              IPOS   = I
              DIST1  = DIST
              ACTVAL = ARRAY(I)
            ENDIF
          ENDIF
10      CONTINUE
        RETURN
        END