C @(#)necnorm.for	19.1 (ES0-DMD) 02/25/03 14:20:24
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
      PROGRAM ECHNRM
C
C+
C
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  12:14 - 11 JAN 1988
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C.IDENTIFICATION
C
C  program  ECHNORM  version 1.5 830610
C
C.MODIFICATIONS
C
C  910128   P. Ballester  Define and set variables (test mmake -i)
C
C.KEYWORDS
C
C  ECHELLE, CASPEC, RIPPLE NORMALIZATION
C
C.PURPOSE
C
C  FIT A POLYNOMIAL TO THE IMAGE ROWS
C
C.ALGORITHM
C
C
C.INPUT/OUTPUT
C
C  NORMALIZE/ECHELLE INPUT OUTPUT TABLE METHOD PARS
C                                    if METHOD='POLY' then PARS=IDEG
C-
C
      IMPLICIT  NONE
      INTEGER   MADRID, TID
      INTEGER   I,IAV,IDEG
      INTEGER*8 IPNTRA,IPNTRB
      INTEGER*8 IS,ISTAT
      INTEGER   J,NAXISA,NCOL,NITER,IMNO,IMNI
      REAL CONST,ER,RELAX,RIDEN,STEPA1
      REAL WSTR1,XMAX1,XMIN1
C
      INTEGER NPIXA(2),IC(20),KUN,KNUL
      REAL RPAR(5),VALUE(20),P0(3)
      REAL XMIN,XMAX,LHCUTS(4)
      DOUBLE PRECISION STEPA(2),STARTA(2)
      DOUBLE PRECISION A(10),DSTP(2),DSTA(2)
      CHARACTER*8 OUTFRM,INFRM,TABLE
      CHARACTER*16 UNIT(5),LABEL(5),UNIT1(12),LABEL1(12),LABELT,UNITT
      CHARACTER*6 FORM(5),FORM1(12),FORMT
      CHARACTER METH*1
      CHARACTER IDENTA*72,CUNITA*64
C
      INCLUDE   'MID_INCLUDE:ST_DEF.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE   'MID_INCLUDE:ST_DAT.INC'
C
      DATA NCOL/5/
      DATA LABEL(1)/'ORDER'/,UNIT(1)/'UNITLESS'/,FORM(1)/'I5'/
      DATA LABEL(2)/'Y0'/,UNIT(2)/'DN'/,FORM(2)/'E13.3'/
      DATA LABEL(3)/'ALFA'/,UNIT(3)/'UNITLESS'/,FORM(3)/'E13.3'/
      DATA LABEL(4)/'BETA'/,UNIT(4)/'UNITLESS'/,FORM(4)/'E13.4'/
      DATA LABEL(5)/'ERROR'/,UNIT(5)/'UNITLESS'/,FORM(5)/'E13.4'/
      DATA LABEL1(1)/'ORDER'/,UNIT1(1)/'UNITLESS'/,FORM1(1)/'I5'/
      DATA LABEL1(2)/'DEGREE'/,UNIT1(2)/'UNITLESS'/,FORM1(2)/'I4'/
      DATA LABEL1(3)/'A0'/,UNIT1(3)/'UNITLESS'/,FORM1(3)/'E13.3'/
      DATA LABEL1(4)/'A1'/,UNIT1(4)/'UNITLESS'/,FORM1(4)/'E13.3'/
      DATA LABEL1(5)/'A2'/,UNIT1(5)/'UNITLESS'/,FORM1(5)/'E13.3'/
      DATA LABEL1(6)/'A3'/,UNIT1(6)/'UNITLESS'/,FORM1(6)/'E13.3'/
      DATA LABEL1(7)/'A4'/,UNIT1(7)/'UNITLESS'/,FORM1(7)/'E13.3'/
      DATA LABEL1(8)/'A5'/,UNIT1(8)/'UNITLESS'/,FORM1(8)/'E13.3'/
      DATA LABEL1(9)/'A6'/,UNIT1(9)/'UNITLESS'/,FORM1(9)/'E13.3'/
      DATA LABEL1(10)/'A7'/,UNIT1(10)/'UNITLESS'/,FORM1(10)/'E13.3'/
      DATA LABEL1(11)/'A8'/,UNIT1(11)/'UNITLESS'/,FORM1(11)/'E13.3'/
      DATA LABEL1(12)/'A9'/,UNIT1(12)/'UNITLESS'/,FORM1(12)/'E13.3'/
C
C ... INITIALIZE SYSTEM
C
      CALL STSPRO('ECHNRM')
      CALL STKRDC('P1',1,1,8,IAV,INFRM,KUN,KNUL,ISTAT)
      CALL STKRDC('P2',1,1,8,IAV,OUTFRM,KUN,KNUL,ISTAT)
      CALL STKRDC('P3',1,1,8,IAV,TABLE,KUN,KNUL,ISTAT)
      CALL STKRDR('INPUTR',1,3,IAV,RPAR,KUN,KNUL,ISTAT)
      METH   = 'P'
      IF (METH.EQ.'S') THEN
          NITER  = NINT(RPAR(3))
      ELSE
          IDEG   = NINT(RPAR(1))
          IDEG   = MIN(IDEG,9)
          NCOL   = IDEG + 3
          VALUE(2) = IDEG
      END IF
C
C ... MAP INPUT FRM
C
      CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .            2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA,
     .            IPNTRA,IMNI,ISTAT)
C
C ... MAP OUTPUT FRM
C
      IDENTA = 'FITTED VALUES :'//INFRM
      DSTA(1) = STARTA(1)
      DSTA(2) = STARTA(2)
      DSTP(1) = STEPA(1)
      DSTP(2) = STEPA(2)
      CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .            NAXISA,NPIXA,DSTA,DSTP,IDENTA,
     .            CUNITA,IPNTRB,IMNO,ISTAT)
C      CALL SXIPUT(OUTFRM,'O',NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA,
C    +            IPNTRB,ISTAT)
C
C ... INITIALIZE TABLE
C
      CALL TBTINI(TABLE,F_TRANS,F_O_MODE,15,NPIXA(2),TID,ISTAT)
      DO 10 I = 1,NCOL
          IF (METH.EQ.'S') THEN
              UNITT  = UNIT(I)
              LABELT = LABEL(I)
              FORMT  = FORM(I)
          ELSE
              UNITT  = UNIT1(I)
              LABELT = LABEL1(I)
              FORMT  = FORM1(I)
          END IF
          CALL TBCINI(TID,D_R4_FORMAT,1,
     .                FORMT,UNITT,LABELT,IC(I),ISTAT)
   10 CONTINUE
C
C ... ITERATION ON EACH ECHELLE ORDER
C
      IS     = 0  !   DUMPING FACTOR
      RELAX  = RPAR(1)  !   ERROR
      ER     = RPAR(2)
      RIDEN  = 1.0
      STEPA1 = STEPA(1)
      XMIN   = 0.
      XMAX   = XMIN
      DO 30 I = 1,NPIXA(2)
          WSTR1 = WSTR1  !   ORDER NUMBER
          VALUE(1) = STARTA(2) + (I-1)*STEPA(2)
          IF (METH.EQ.'S') THEN  !   INITAL GUESSES
              P0(2)  = 1.
C     P0(3) = 1./(WSTR(I)+NPIXA(1)*STEPA(1)*0.5)
              P0(3)  = 1./ (WSTR1+NPIXA(1)*STEPA1*0.5)
              CONST  = VALUE(1)*3.1415926535897932/10800.0
              CALL FIT1(NPIXA(1),NPIXA(2),MADRID(IPNTRA),MADRID(IPNTRB),
     +                  I,CONST,P0,ER,RELAX,NITER,WSTR1,STEPA1,
     +                  VALUE(2),VALUE(4))
          ELSE
              CALL FIT2(NPIXA(1),NPIXA(2),MADRID(IPNTRA),MADRID(IPNTRB),
     +                  I,STARTA(1),STEPA(1),IDEG,A,XMAX1,XMIN1)
              DO 20 J = 1,IDEG + 1
                  VALUE(J+2) = A(J)
   20         CONTINUE
          END IF
          XMIN   = AMIN1(XMIN,XMIN1)
          XMAX   = AMAX1(XMAX,XMAX1)
          IS     = IS + 1
          CALL TBRWRR(TID,IS,NCOL,IC,VALUE,ISTAT)
   30 CONTINUE
      LHCUTS(1) = XMIN
      LHCUTS(2) = XMAX
      LHCUTS(3) = XMIN
      LHCUTS(4) = XMAX
      CALL STDWRR(IMNO,'LHCUTS',LHCUTS,1,4,KUN,ISTAT)
      CALL TBTCLO(TID,ISTAT)
C
C ... END
C
      CALL STSEPI
      END

      SUBROUTINE FIT1(NX,NY,IN,OUT,I,CONST,P0,ER,RELAX,NITER,X0,DX,P,
     +                ERROR)
C
C NONLINEAR FIT ROUTINE INTERFACE
C      IMPLICIT NONE
C
C INPUT ARGUMENTS
C
C NX INTG*4 NUMBER OF SAMPLES
C  NY INTG*4 NUMBER OF LINES (ORDERS)
C IN REAL*4 INPUT IMAGE
C OUT REAL*4 FITTED VALUES
C I INTG*4 NUMBER OF THE LINE
C M REAL*4 SPECTRAL ORDER
C P0(3) REAL*4 INITIAL GUESSES
C ER REAL*4 ERROR BOUND
C RELAX REAL*4 DUMPING FACTOR
C NITER INTG*4 NUMBER OF ITERATIONS
C X0 REAL*4 STARTING WAVELENGTH
C DX REAL*4 WAVELENGTH STEP
C
C OUTPUT PARAMETERS
C
C P(3) REAL*4 ESTIMATED PARAMETER
C ERROR REAL*4 FINAL ERROR
C
      IMPLICIT NONE

      INTEGER  NX,NY,I,NITER,ISTAT
      REAL     CONST,ER,RELAX,X0,DX,ERROR

      REAL IN(NX,NY),OUT(NX,NY),P0(3),P(3)
      CHARACTER*80 LINE

      P0(1)  = IN(NX/2,I)
      WRITE (LINE,9000) I,P0(1),P0(2),P0(3),ER
      CALL STTPUT(LINE,ISTAT)
      CALL LSQSIN(NX,X0,DX,IN(1,I),OUT(1,I),RELAX,NITER,ER,CONST,P0,P,
     +             ERROR)
      RETURN

 9000 FORMAT (' Line :',I4,' Init.pars :',3E10.4,' Error :',E10.4)
      END

      SUBROUTINE LSQSIN(NP,X0,DX,Y,YFIT,RELAX,NITER,ER,CONST,P0,P,
     +                   ERROR)
C
      IMPLICIT NONE
C     Modified version of the program LSQFIT for the problem :
C
C     find A and B in the expression
C     Y = (sinc(pi*M*A*(X*B-1)))**2
C
      INTEGER    NP,NITER,ITER,INDEP,KZ6901,IPAR,ISTAT,I,J,K,NF
      REAL       X0,DX,RELAX,ER,CONST,ERROR,X,Y2,W,DET,ERR
      REAL       Y(NP),YFIT(NP),P(3),P0(3),D(3)
      DOUBLE PRECISION CHI,PP,Q
      DOUBLE PRECISION AL(3),B(3),A(3,3)
      CHARACTER*80     LINE
C
C ... initialize relevant parameters
C
      ITER   = 0
C     IWEIGH,RELAX  !   one independent variable
      INDEP  = 1  !   three parameters
      IPAR   = 3  !   initial guesses
      P(1)   = P0(1)
      P(2)   = P0(2)
      P(3)   = P0(3)
      NF     = NP - 3
C
C ... iteration loop
C
      DO 80 KZ6901 = 1,2147483647
         IF( .NOT. (ITER.LT.NITER)) GO TO 90
         ITER   = ITER + 1
         WRITE (LINE,9000) ITER,P(1),P(2),P(3)
         CALL STTPUT(LINE,ISTAT)
C
C ... initalize normal equation matrix
C
         DO 20 I = 1,IPAR
             B(I)   = 0.0D0
             DO 10 J = 1,IPAR
                 A(I,J) = 0.0D0
   10        CONTINUE
   20    CONTINUE
         CHI    = 0.0D0
         DO 50 I = 1,NP
C
C ... compute function and derivatives
C
             X      = X0 + (I-1)*DX
             CALL FSINC(X,YFIT(I),CONST,P,D)
C
C ... compute residuals
C
             Y2     = Y(I) - YFIT(I)
             W      = 1.0
             PP     = DBLE(W)
             Q      = DBLE(Y2)
             CHI    = CHI + Q*Q*PP
C
C ... enter derivatives into normal equation matrix
C
             DO 40 J = 1,IPAR
                 AL(1)  = DBLE(D(J))
                 DO 30 K = 1,IPAR
                     A(J,K) = A(J,K) + AL(1)*DBLE(D(K))*PP
   30            CONTINUE
                 B(J)   = B(J) + PP*Q*AL(1)
   40        CONTINUE
   50    CONTINUE
C
C ... invert normal equation
C
         WRITE (*,*) A(1,1),A(1,2),A(1,3),B(1)
         WRITE (*,*) A(2,1),A(2,2),A(2,3),B(2)
         WRITE (*,*) A(3,1),A(3,2),A(3,3),B(3)
C        DO K=1,IPAR
C           KK=IPAR+1-K
C           PP=A(1,1)
C           DO I=2,IPAR
C              Q=A(I,1)
C              AL(I)=Q/PP
C              IF (I.LE.KK) AL(I)=-AL(I)
C              DO J=2,I
C          IF (I.EQ.2.AND.J.EQ.2)WRITE(*,*) A(2,2)
C                 A(I-1,J-1)=A(I,J)+Q*AL(J)
C          IF (I.EQ.2.AND.J.EQ.2)WRITE(*,*) K,Q,AL(J),A(1,1)
C              ENDDO
C           ENDDO
C           A(IPAR,IPAR)=1.0D0/PP
C           DO I=2,IPAR
C              A(IPAR,I-1)=AL(I)
C           ENDDO
C        ENDDO
C
C ... fill matrix
C
C        DO I=2,IPAR
C           DO J=1,I
C              A(J,I)=A(I,J)
C           ENDDO
C        ENDDO
         CALL MATINV(A,3,DET)
         IF (DET.EQ.0.0) THEN
             WRITE (*,*) 'ERROR ------------------------------'
         END IF

         WRITE (*,*) DET
         WRITE (*,*) A(1,1),A(1,2),A(1,3)
         WRITE (*,*) A(2,1),A(2,2),A(2,3)
         WRITE (*,*) A(3,1),A(3,2),A(3,3)
C
C ... compute reduced chi**2 and corrections for unknowns
C
         CHI    = CHI/DBLE(NF)
         DO 70 I = 1,IPAR
             Q      = 0.0D0
             DO 60 J = 1,IPAR
                Q      = Q + A(I,J)*B(J)
   60        CONTINUE
C
C ... compute errors of the corrections
C
            D(I)   = SNGL(DSQRT(DMAX1(A(I,I)*CHI,0.0D0)))
            P(I)   = P(I) + RELAX*SNGL(Q)  !   DEBUG <------------- !
            WRITE (*,*) D(I),P(I)
   70    CONTINUE
         ERROR  = SNGL(DSQRT(DMAX1(CHI,0.0D0)))
         LINE   = ' '
         WRITE (LINE,9010) ERROR
         CALL STTPUT(LINE,ISTAT)
         IF (ABS(ERROR).LE.ERR) RETURN
   80 CONTINUE
   90 CONTINUE
      RETURN

 9000 FORMAT (' Iteration :',I4,' Params :',3E13.4)
 9010 FORMAT (' Error     :',E13.4)
      END

      SUBROUTINE FSINC(X,Y1,CONST,P,D)
C
      IMPLICIT NONE
C     INTERNAL FUNCTION SINC, SINGULARITY AT (X*P(3)->1.)
C
      REAL      CONST,X,Y1,A0,A,A2,SINCA,CA,SA,P,D
      DIMENSION P(3),D(3)

      A0     = X*P(3) - 1.
      IF (ABS(A0).LT.0.1E-20) THEN
          Y1     = 1.
          D(1)   = 1.
          D(2)   = 0.
          D(3)   = 0.

      ELSE
          A      = CONST*P(2)*A0
          A2     = A*A
          SINCA  = SIN(A)/A
          CA     = COS(A) - SA
          D(1)   = SINCA*SINCA
          Y1     = P(1)*D(1)
          D(2)   = (2.*P(1)*SINCA*CA)/P(2)
          D(3)   = (2.*P(1)*SINCA*CA*X)/A0
      END IF

      RETURN

      END

      SUBROUTINE FIT2(NX,NY,IN,OUT,I,X0,DX,IDEG,A,XMIN,XMAX)
C      IMPLICIT NONE
C
C LINEAR FIT ROUTINE INTERFACE
C
C INPUT ARGUMENTS
C
C NX INTG*4 NUMBER OF SAMPLES
C  NY INTG*4 NUMBER OF LINES (ORDERS)
C IN REAL*4 INPUT IMAGE
C OUT REAL*4 FITTED VALUES
C I INTG*4 NUMBER OF THE LINE
C IDEG INTG*4 POLYNOMIAL DEGREE
C
C OUTPUT PARAMETERS
C
C A REAL*8 COEFFS
C
      IMPLICIT NONE
      INTEGER          NX,NY,I,IDEG,NT,J
      REAL             X0,DX,XMIN,XMAX
      REAL             IN(NX,NY),OUT(NX,NY)
      DOUBLE PRECISION A(1)

      NT     = NX - NX/50
      CALL LSQFIT(NT,X0,DX,IN(1,I),OUT(1,I),IDEG,A,XMIN,XMAX)
      DO 10 J = NT + 1,NX
          OUT(J,I) = 0.
   10 CONTINUE
      RETURN

      END

      SUBROUTINE LSQFIT(NX,X0,DX,X,FIT,K,A,XMIN,XMAX)
      IMPLICIT NONE
C
C SOLVE LEAST SQUARE PROBLEM
C
      INTEGER          N,NX,K,L,LL,N1,I,KMIN,K1,J
      REAL             X0,DX,XMIN,XMAX,X1,X2,Y,ERROR
      REAL             X(NX),FIT(NX)
      DOUBLE PRECISION G,XC,A(1)
      COMMON /LS/G(50,50),XC(50),N
C
C ITERATION
C
      L      = 0
      LL     = 0
      N1     = 0
      N      = (K+1)* (L+1)
      DO 20 I = 10,NX - 20
          X1     = X0 + (I-1)*DX
          X2     = X1
          Y      = X(I)
          N1     = N1 + 1
          CALL SETUP(LL+1,X1,X2,Y,K,L)
          IF (LL.NE.0) THEN
              KMIN   = MIN(LL,N+1)
              DO 10 K1 = 1,KMIN
                  CALL HT(K1,LL+1)
   10         CONTINUE
          END IF

          LL     = MIN(LL+1,N+1)
   20 CONTINUE
      CALL SOLVE
      ERROR  = ABS(G(N+1,N+1))
      CALL COMP(NX,X0,DX,FIT,K,L,XMIN,XMAX)
      DO 30 J = 1,K + 1
          A(J)   = XC(J)
   30 CONTINUE
      RETURN

      END

      SUBROUTINE SETUP(LL,X1,X2,Y,K,L)
      IMPLICIT NONE
C
C PREPARE DATA IN WORK SPACE
C
      INTEGER          N,LL,K,L,IP,L1,K1
      REAL             X1,X2,Y
      DOUBLE PRECISION G,X,DC
      COMMON /LS/G(50,50),X(50),N

      IP     = 0
      DC     = 1.D0
C DX1 = DBLE(X1)*0.001D0
C DX2 = DBLE(X2)*0.001D0
C DY  = DBLE(Y)*0.001D0
      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 HT(K,L1)
      IMPLICIT NONE

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.
      DOUBLE PRECISION G,X,S,H,B
      COMMON /LS/G(50,50),X(50),N

C     Construct the transformation
      INTEGER   N,K,L1,KH,J 

      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     Apply the transformation to the rows K and L1 of G.

      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

      SUBROUTINE SOLVE
      IMPLICIT NONE
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)
      DOUBLE PRECISION G,X,S
      COMMON /LS/G(50,50),X(50),N

      INTEGER   N,I,II,J

      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

      SUBROUTINE COMP(NPIX1,START,STEP,OUT,K,L,XMIN,XMAX)
      IMPLICIT NONE
C

      INTEGER          N,NPIX1,K,L,IX,IP,L1,K1
      REAL             START,STEP,XMIN,XMAX
      REAL             OUT(NPIX1)
      DOUBLE PRECISION WORK(50),DC,YVAL,G,A
      DOUBLE PRECISION X1,X2,XSTART,XSTEP
      COMMON /LS/G(50,50),A(50),N
C
C  initialize variables
      XMIN   = 0.
      XMAX   = XMIN
      XSTART = START
      XSTEP  = STEP
C
C  start loop
      DO 30 IX = 1,NPIX1
          X1     = (IX-1)*XSTEP + XSTART
          IP     = 0
          DC     = 1.D0
          YVAL   = 0.D0
          DO 20 L1 = 0,L
              IP     = IP + 1
              WORK(IP) = DC
              YVAL   = YVAL + WORK(IP)*A(IP)
              DO 10 K1 = 1,K
                  IP     = IP + 1
                  WORK(IP) = WORK(IP-1)*X1
                  YVAL   = YVAL + WORK(IP)*A(IP)
   10         CONTINUE
              DC     = DC*X2
   20     CONTINUE
          OUT(IX) = YVAL
          XMIN   = AMIN1(XMIN,OUT(IX))
          XMAX   = AMAX1(XMAX,OUT(IX))
   30 CONTINUE
      RETURN

      END
