C===========================================================================
C Copyright (C) 1995-2005 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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  13:19 - 15 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C. Modif 30-Jul-90 P.Ballester
C 
C.VERSION
C 051017		last modif
C 
C-------------------------------------------------------------------

      SUBROUTINE SMOOT(X,Y,W,Q,M,XB,XE,K,S,N,T,C,FP,IOPT,IER,NRDATA,
     +                 FP0,FPOLD,NPLUS)
C
       IMPLICIT NONE
C
C  GIVEN THE SET OF DATA POINTS (X(I),Y(I)) AND THE SET OF POSITIVE
C  NUMBERS W(I),I=1,2,...M, SMOOT DETERMINES A SMOOTH SPLINE
C  APPROXIMATION OF DEGREE K FOR THE INTERVAL XB <= X <= XE.
C  THE NUMBER OF KNOTS OF THIS SPLINE S(X) AND THE POSITION
C  T(J),J=1,...N IS CHOSEN AUTOMATICALLY BY THE ROUTINE.
C  THE SMOOTHNESS OF S(X) IS ACHIEVED BY MINIMALIZING THE SUM(D(R)**2)
C  WHERE D(R) STANDS FOR THE DISCONTINUITY JUMP OF THE KTH DERIVATIVE
C  OF S(X) AT THE KNOT T(R),R=K+2,...,N-K-1. THE AMOUNT OF SMOOTH-
C  NESS IS DETERMINED BY THE CONDITION THAT F(P)=SUM(W(I)*(Y(I)-
C  S(X(I)))**2) BE <= S, WITH S A GIVEN NON-NEGATIVE CONSTANT.
C  THE FIT S(X) IS GIVEN IN ITS B-SPLINE REPRESENTATION ( B-SPLINE
C  COEFFICIENTS C(J),J=1,2,...N-K-1) AND CAN BE EVALUATED BY MEANS
C  OF FUNCTION DERIV.
C
C  CALLING SEQUENCE:
C     CALL SMOOT(X,Y,W,Q,M,XB,XE,K,S,N,T,C,FP,IOPT,IER)
C
C  INPUT PARAMETERS:
C    X   : ARRAY,MINIMUM LENGTH M,CONTAINING THE DATA POSITIONS X(I).
C    Y   : ARRAY,MINIMUM LENGTH M,CONTAINING THE DATA VALUES Y(I).
C    W   : ARRAY,MINIMUM LENGTH M,CONTAINING THE WEIGHTS W(I).
C    M   : INTEGER VALUE,CONTAINING THE NUMBER OF DATA POINTS.
C    XB,XE:REAL VALUES, CONTAINING THE BOUNDARIES OF THE APPR. INTERVAL
C    K   : INTEGER VALUE,CONTAINING THE DEGREE OF S(X).
C    S   : REAL VALUE,CONTAINING THE SMOOTHING FACTOR.
C    IOPT: INTEGER VALUE WHICH TAKES THE VALUE 0 OR 1.
C      IOPT=0: THE ROUTINE WILL RESTART ALL COMPUTATIONS.
C      IOPT=1: THE ROUTINE WILL START WITH THE KNOTS FOUND AT THE
C               LAST CALL OF THE ROUTINE. IF IOPT=1 THE OUTPUT
C               PARAMETERS T AND N ARE INPUT PARAMETERS AS WELL.
C
C  OUTPUT PARAMETERS:
C    T   : ARRAY,LENGTH NEST (SEE DATA INITIALIZATION STATEMENT),
C          WHICH CONTAINS THE POSITION OF THE KNOTS,I.E. THE POSITION
C          OF THE INTERIOR KNOTS T(K+2),...T(N-K-1), AS WELL AS THE
C          POSITION OF THE KNOTS T(1)=T(2)=...=T(K+1)=XB AND XE =
C          T(N-K)=...=T(N) WHICH ARE NEEDED FOR THE B-SPLINE REPRESENT.
C    C   : ARRAY,LENGTH NEST,WHICH CONTAINS THE B-SPLINE COEFFICIENTS.
C    N   : INTEGER VALUE,CONTAINING THE TOTAL NUMBER OF KNOTS.
C    FP  : REAL VALUE,WHICH CONTAINS THE SUM(WI*(YI-S(XI))**2),I=1,...M
C    IER : ERROR CODE
C      IER=0: NORMAL RETURN.
C      IER=-1:NORMAL RETURN, S(X) IS AN INTERPOLATING SPLINE.
C      IER=-2:NORMAL RETURN, S(X) IS A POLYNOMIAL OF DEGREE K.
C      IER>0 :ABNORMAL TERMINATION.
C        IER=1:THE REQUIRED STORAGE SPACE EXCEEDS THE AVAILABLE
C              STORAGE SPACE,SPECIFIED BY THE PARAMETER NEST.
C              PROBABLY CAUSES:NEST OR S TOO SMALL.
C        IER=2:A THEORETICALLY IMPOSSIBLE RESULT WAS FOUND DURING
C              THE ITERATION PROCESS.
C              PROBABLY CAUSES:TOL TOO SMALL
C        IER=3:THE MAXIMAL NUMBER OF ITERATIONS MAXIT HAS BEEN
C              REACHED.
C              PROBABLY CAUSES:MAXIT OR TOL TOO SMALL.
C        IER=10:SOME OF THE INPUTDATA ARE INVALID(SEE RESTRICTIONS).
C
C
C   INPUT/OUTPUT PARAMETERS:
C     NRDATA : ARRAY, LENGTH NEST KEEP INFO BETWEEN ITERATIONS
C     FPO    : REAL, KEEP INFO
C     FPOLD  : REAL, KEEP INFO
C     NPLUS  : REAL, KEEP INFO
C
C  RESTRICTIONS:
C    1) M > K > 0
C    2) XB <= X(R) < X(R+1) <= XE, R=1,2,...M-1.
C    3) W(R) > 0, R=1,2,...M.
C    4) S >= 0.
C    5) NEST >= 2*K+2.
C
C  OTHER SUBROUTINES REQUIRED:
C    BSPLIN,COSSIN,ROTATE,BACK,NKNOT,DISCO AND RATION.
      INTEGER MXP
      PARAMETER (MXP=1000)
      INTEGER M,N,IOPT,IER,K,I,J,L,NMAX,NEST,MAXIT,NEW,ICHECK,N8
      INTEGER NPLUS,MK1,K2,K3,K1,NMIN,ITER,NRINT,NK1,IT,I1,I2,I3
      INTEGER NPL1,L0,NRDATA(MXP)
      REAL X(M),Y(M),W(M),T(MXP),C(MXP),FPINT(MXP),
     +          Z(MXP),DIAG(MXP),DPRIME(MXP),A(MXP,5),G(MXP,6),B(MXP,7),
     +          Q(M,6),H(7),FP,XE,XB,S,FP0,FPOLD,TOL
      REAL ACC,XI,YI,WI,PIV,FPART,TERM,STORE,P1
      REAL PINV,P,COS,SIN,FPMS,F1,F2,F3,P2,P3,RATION
C
C    NRDATA: INTEGER ARRAY,LENGTH NEST,WHICH GIVES THE NUMBER OF
C            DATA POINTS INSIDE EACH KNOT INTERVAL.
C    FP0   : REAL VALUE,WHICH CONTAINS THE WEIGHTED SUM OF SQUARES
C            OF RESIDUAL RIGHT HAND SIDES, FOR THE LEAST SQUARES
C            POLYMOMIAL OF DEGREE K.
C    FPOLD : REAL VALUE,WHICH CONTAINS THE WEIGHTED SUM OF RESIDUAL
C            RIGHT HAND SIDES FOR THE LEAST SQUARES SPLINE WHICH
C            CORRESPONDS TO THE LAST  FOUND SET OF KNOTS BUT ONE.
C    NPLUS : INTEGER VALUE,GIVING THE NUMBER OF KNOTS OF THE LAST
C            SET MINUS THE NUMBER OF THE LAST SET BUT ONE.
C
C  DATA INITIALIZATION STATEMENT TO SPECIFY
C    TOL  : THE REQUESTED RELATIVE ACCURACY FOR THE ROOT OF F(P) = S.
C    MAXIT: THE MAXIMAL NUMBER OF ITERATIONS ALLOWED.
C    NEST : AN OVER-ESTIMATE OF THE NUMBER OF KNOTS N. THIS PARAMETER
C          MUST BE SET BY THE USER TO INDICATE THE STORAGE SPACE
C          AVAILABLE TO THE SUBROUTINE. THE DIMENSION SPECIFICATIONS
C          OF THE ARRAYS   T,C,NRDATA,FPINT,Z,DIAG,DPRIME(N),A(N,K),
C          Q(M,K+1),G(N,K+1),B(N,K+2) AND H(K+2)   DEPEND
C          ON N,M AND K. SINCE N IS UNKNOWN AT THE TIME THE
C          USER SETS UP THE DIMENSION INFORMATION AN OVER-ESTIMATE
C          OF THESE ARRAYS WILL GENERALLY BE MADE. THE FOLLOWING
C          REMARKS ARE INTENDED TO HELP THE USER
C            1) 2*K+2 <= N <= M+K+1
C            2) THE SMALLER THE VALUE OF S, THE GREATER N WILL BE.
C            3) NORMALLY N = M/2 IS AN OVER-ESTIMATE.
      DATA TOL/0.001/,MAXIT/40/,NEST/MXP/ ! Old values TOL=0.001 and MAXIT=40
C  BEFORE STARTING COMPUTATIONS A DATA CHECK IS MADE. IF THE INPUT
C  DATA ARE INVALID CONTROLE IS IMMEDIATELY REPASSED TO THE DRIVER
C  PROGRAM (IER=10).
 
      FPMS = 0.0
   10 IER    = 0
      K1     = K + 1
      K2     = K1 + 1
      NMIN   = 2*K1
      IF (K.LE.0) IER    = 10
      IF (M.LT.K1 .OR. NEST.LT.NMIN) IER    = 10
      IF (S.LT.0.) IER    = 10
      IF (XB.GT.X(1) .OR. XE.LT.X(M) .OR. W(1).LE.0.) IER    = 10
      IF (IER.NE.0) GO TO 400
C  CALCULATION OF ACC, THE ABSOLUTE TOLERANCE FOR THE ROOT OF F(P)=S.
      ACC    = TOL*S
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  PART 1: DETERMINATION OF THE NUMBER OF KNOTS AND THEIR POSITION     C
C  **************************************************************      C
C  GIVEN A SET OF KNOTS WE COMPUTE THE LEAST-SQUARES SPLINE SINF(X).   C
C  IF THE SUM F(P=INF) <= S WE ACCEPT THE CHOICE OF KNOTS.             C
C  OTHERWISE WE HAVE TO INCREASE THEIR NUMBER.                         C
C  THE INITIAL CHOICE OF KNOTS DEPENDS ON THE VALUE OF S AND IOPT.     C
C    IF S=0 WE HAVE SPLINE INTERPOLATION; IN THAT CASE THE NUMBER OF   C
C    KNOTS EQUALS NMAX = M+K+1.                                        C
C    IF S > 0 AND                                                      C
C      IOPT=0 WE FIRST COMPUTE THE LEAST-SQUARES POLYNOMIAL OF         C
C      DEGREE K; N = NMIN = 2*K+2                                      C
C      IOPT=1 WE START WITH THE SET OF KNOTS FOUND AT THE LAST         C
C      CALL OF THE ROUTINE, EXCEPT FOR THE CASE THAT S > FP0; THEN     C
C      WE COMPUTE DIRECTLY THE LEAST-SQUARES POLYNOMIAL OF DEGREE K.   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  DETERMINE NMAX, THE NUMBER OF KNOTS FOR SPLINE INTERPOLATION.
      NMAX   = M + K1
      IF (S.GT.0.) GO TO 50
C  IF S=0, S(X) IS AN INTERPOLATING SPLINE.
      N      = NMAX
C  TEST WHETHER THE REQUIRED STORAGE SPACE EXCEEDS THE AVAILABLE ONE.
      IF (N.GT.NEST) GO TO 380
C  FIND THE POSITION OF THE INTERIOR KNOTS IN CASE OF INTERPOLATION.
      MK1    = M - K1
      IF (MK1.EQ.0) GO TO 70
      K3     = K/2
      I      = K2
      J      = K3 + 2
      IF (K3*2.EQ.K) GO TO 30
      DO 20 L = 1,MK1
          T(I)   = X(J)
          I      = I + 1
          J      = J + 1
   20 CONTINUE
      GO TO 70

   30 DO 40 L = 1,MK1
          T(I)   = (X(J)+X(J-1))*0.5    ! 30-JUl-90
          I      = I + 1
          J      = J + 1
   40 CONTINUE
      GO TO 70
C  IF S>0 OUR INITIAL CHOICE OF KNOTS DEPENDS ON THE VALUE OF IOPT.
C  IF IOPT=0 OR IOPT=1 AND S>=FP0, WE START COMPUTING THE LEAST-SQUARES
C  POLYNOMIAL OF DEGREE K WHICH IS A SPLINE WITHOUT INTERIOR KNOTS.
C  IF IOPT=1 AND FP0>S WE START COMPUTING THE LEAST SQUARES SPLINE
C  ACCORDING TO THE SET OF KNOTS FOUND AT THE LAST CALL OF THE ROUTINE.
   50 IF (IOPT.LE.0) GO TO 60
      IF (FP0.GT.S) GO TO 70
   60 N      = NMIN
      NRDATA(1) = M - 2
C  MAIN LOOP FOR THE DIFFERENT SETS OF KNOTS. M IS A SAVE UPPER BOUND
C  FOR THE NUMBER OF TRIALS.
   70 DO 220 ITER = 1,M
          IF (N.EQ.NMIN) IER    = -2
C  FIND NRINT, TNE NUMBER OF KNOT INTERVALS.
          NRINT  = N - NMIN + 1
C  FIND THE POSITION OF THE ADDITIONAL KNOTS WHICH ARE NEEDED FOR
C  THE B-SPLINE REPRESENTATION OF S(X).
          NK1    = N - K1
          I      = N
          DO 80 J = 1,K1
              T(J)   = XB
              T(I)   = XE
              I      = I - 1
   80     CONTINUE
C  COMPUTE THE B-SPLINE COEFFICIENTS OF THE LEAST-SQUARES SPLINE
C  SINF(X). THE OBSERVATION MATRIX A IS BUILT UP ROW BY ROW AND
C  REDUCED TO UPPER TRIANGULAR FORM BY GIVENS TRANSFORMATIONS
C  WITHOUT SQUARE ROOTS. AT THE SAME TIME FP=F(P=INF) IS COMPUTED.
          FP     = 0.
C  INITIALIZE THE OBSERVATION MATRIX A.
          DO 100 I = 1,NK1
              DIAG(I) = 0.
              Z(I)   = 0.
              DO 90 J = 1,K
                  A(I,J) = 0.
   90         CONTINUE
  100     CONTINUE
          L      = K1
          DO 150 IT = 1,M
C  FETCH THE CURRENT DATA POINT X(IT),Y(IT).
              XI     = X(IT)
              YI     = Y(IT)
              WI     = W(IT)
C  SEARCH FOR KNOT INTERVAL T(L) <= XI < T(L+1).
              IF (XI.GE.T(L+1) .AND. L.NE.NK1) L      = L + 1
C  EVALUATE THE (K+1) NON-ZERO B-SPLINES AT XI AND STORE THEM IN Q.
              CALL BSPLIN(T,N,K,XI,L,H)
              DO 110 I = 1,K1
                  IF (H(I).LT.0.1E-07) H(I)   = 0.
                  Q(IT,I) = H(I)
  110         CONTINUE
C  ROTATE THE NEW ROW OF THE OBSERVATION MATRIX INTO TRIANGLE BY
C  GIVENS TRANSFORMATIONS WITHOUT SQUARE ROOTS.
              J      = L - K1
              DO 130 I = 1,K1
                  IF (WI.LE.1.E-30) GO TO 150
                  J      = J + 1
                  PIV    = H(I)
                  IF (PIV.EQ.0.) GO TO 130
C  CALCULATE THE PARAMETERS OF THE GIVENS TRANSFORMATION.
                  CALL COSSIN(PIV,WI,DIAG(J),COS,SIN)
C  TRANSFORMATIONS TO RIGHT HAND SIDE.
                  CALL ROTATE(PIV,COS,SIN,YI,Z(J))
                  IF (I.EQ.K1) GO TO 140
                  I2     = 0
                  I3     = I + 1
                  DO 120 I1 = I3,K1
                      I2     = I2 + 1
C  TRANSFORMATIONS TO LEFT HAND SIDE.
                      CALL ROTATE(PIV,COS,SIN,H(I1),A(J,I2))
  120             CONTINUE
  130         CONTINUE
C  ADD CONTRIBUTION OF THIS ROW TO THE SUM OF SQUARES OF RESIDUAL
C  RIGHT HAND SIDES.
  140         FP     = FP + WI*YI**2
  150     CONTINUE
          IF (IER.EQ.-2) FP0    = FP
C  BACKWARD SUBSTITUTION TO OBTAIN THE B-SPLINE COEFFICIENTS.
          CALL BACK(A,Z,NK1,K,C)
C  TEST WHETHER THE APPROXIMATION SINF(X) IS AN ACCEPTABLE SOLUTION.
          FPMS   = FP - S
          IF (ABS(FPMS).LT.ACC) GO TO 400
C  IF F(P=INF) < S ACCEPT THE CHOICE OF KNOTS.
          IF (FPMS.LT.0.) GO TO 230
C  IF N = NMAX, SINF(X) IS AN INTERPOLATING SPLINE.
          IF (N.EQ.NMAX) GO TO 390
C  INCREASE THE NUMBER OF KNOTS.
C  IF N=NEST WE CANNOT INCREASE THE NUMBER OF KNOTS BECAUSE OF
C  THE STORAGE CAPACITY LIMITATION.
          IF (N.EQ.NEST) GO TO 380
C  DETERMINE THE NUMBER OF KNOTS NPLUS WE ARE GOING TO ADD.
          IF (IER.EQ.0) GO TO 160
          NPLUS  = 1
          IER    = 0
          GO TO 170

  160     NPL1   = NPLUS*2
          IF (FPOLD-FP.GT.ACC) NPL1   = FLOAT(NPLUS)*FPMS/ (FPOLD-FP)
          NPLUS  = MIN0(NPLUS*2,MAX0(NPL1,NPLUS/2,1))
  170     FPOLD  = FP
C  COMPUTE THE SUM(W(I)*(Y(I)-S(X(I)))**2) FOR EACH KNOT INTERVAL
C  T(J+K) <= X(I) <= T(J+K+1) AND STORE IT IN FPINT(J),J=1,2,...NRINT.
          FPART  = 0.
          I      = 1
          L      = K2
          NEW    = 0
          DO 200 IT = 1,M
              IF (X(IT).LT.T(L) .OR. L.GT.NK1) GO TO 180
              NEW    = 1
              L      = L + 1
  180         TERM   = 0.
              L0     = L - K2
              DO 190 J = 1,K1
                  L0     = L0 + 1
                  TERM   = TERM + C(L0)*Q(IT,J)
  190         CONTINUE
              TERM   = W(IT)* (TERM-Y(IT))**2
              FPART  = FPART + TERM
              IF (NEW.EQ.0) GO TO 200
              STORE  = TERM*0.5
              FPINT(I) = FPART - STORE
              I      = I + 1
              FPART  = STORE
              NEW    = 0
  200     CONTINUE
          FPINT(NRINT) = FPART
          DO 210 L = 1,NPLUS
C  ADD A NEW KNOT.
              CALL NKNOT(X,M,T,N,FPINT,NRDATA,NRINT)
C  TEST WHETHER WE CANNOT FURTHER INCREASE THE NUMBER OF KNOTS.
              IF (N.EQ.NMAX .OR. N.EQ.NEST) GO TO 220
  210     CONTINUE
C  RESTART THE COMPUTATIONS WITH THE NEW SET OF KNOTS.
  220 CONTINUE
C  TEST WHETHER THE LEAST-SQUARES KTH DEGREE POLYNOMIAL IS A SOLUTION
C  OF OUR APPROXIMATION PROBLEM.
  230 IF (IER.EQ.-2) GO TO 400
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  PART 2: DETERMINATION OF THE SMOOTHING SPLINE SP(X).                C
C  ***************************************************                 C
C  WE HAVE DETERMINED THE NUMBER OF KNOTS AND THEIR POSITION.          C
C  WE NOW COMPUTE THE B-SPLINE COEFFICIENTS OF THE SMOOTHING SPLINE    C
C  SP(X). THE OBSERVATION MATRIX A IS EXTENDED BY THE ROWS OF MATRIX   C
C  B EXPRESSING THAT THE KTH DERIVATIVE DISCONTINUITIES OF SP(X) AT    C
C  THE INTERIOR KNOTS T(K+2),...T(N-K-1) MUST BE ZERO. THE CORRES-     C
C  PONDING WEIGHTS OF THESE ADDITIONAL ROWS ARE SET TO 1/SQRT(P).      C
C  ITERATIVELY WE THEN HAVE TO DETERMINE THE VALUE OF P SUCH THAT      C
C  F(P)=SUM(W(I)*(Y(I)-SP(X(I)))**2) BE = S. WE ALREADY KNOW THAT      C
C  THE LEAST-SQUARES KTH DEGREE POLYNOMIAL CORRESPONDS TO P=0, AND     C
C  THAT THE LEAST-SQUARES SPLINE CORRESPONDS TO P=INFINITY. THE        C
C  ITERATION PROCESS WHICH IS PROPOSED HERE, MAKES USE OF RATIONAL     C
C  INTERPOLATION. SINCE F(P) IS A CONVEX AND STRICTLY DECREASING       C
C  FUNCTION OF P, IT CAN BE APPROXIMATED BY A RATIONAL FUNCTION        C
C  R(P) = (U*P+V)/(P+W). THREE VALUES OF P(P1,P2,P3) WITH CORRESPOND-  C
C  ING VALUES OF F(P) (F1=F(P1)-S,F2=F(P2)-S,F3=F(P3)-S) ARE USED      C
C  TO CALCULATE THE NEW VALUE OF P SUCH THAT R(P)=S. CONVERGENCE IS    C
C  GUARANTEED BY TAKING F1>0 AND F3<0.                                 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  EVALUATE THE DISCONTINUITY JUMP OF THE KTH DERIVATIVE OF THE
C  B-SPLINES AT THE KNOTS T(L),L=K+2,...N-K-1 AND STORE IN B.
      CALL DISCO(T,N,K2,B)
C  INITIAL VALUE FOR P.
      P1     = 0.
      P2     = 0.
      F1     = FP0 - S
      P3     = -1.
      F3     = FPMS
      P      = -F1/F3
      ICHECK = 0
      N8     = N - NMIN
C  ITERATION PROCESS TO FIND THE ROOT OF F(P) = S.
      DO 350 ITER = 1,MAXIT
C  THE ROWS OF MATRIX B WITH WEIGHT 1/SQRT(P) ARE ROTATED INTO
C  THE TRIANGULARISED OBSERVATION MATRIX A WHICH IS STORED IN G.
          IF (ABS(P).LT.1.0E-30) THEN     !P = 0   => PINV = 1/P is in trouble
             GOTO 360                     !indicate that iteration failed
          ENDIF   
C 
          PINV   = 1.0/P
          DO 250 I = 1,NK1
              DPRIME(I) = DIAG(I)
              C(I)   = Z(I)
              G(I,K1) = 0.
              DO 240 J = 1,K
                  G(I,J) = A(I,J)
  240         CONTINUE
  250     CONTINUE
          DO 290 IT = 1,N8
C*PL*ERROR* Comment line too long
C   THE ROW OF MATRIX B IS ROTATED INTO TRIANGLE BY GIVENS TRANSFORMATIONS
              DO 260 I = 1,K2
                  H(I)   = B(IT,I)
  260         CONTINUE
              YI     = 0.
              WI     = PINV
              DO 280 J = IT,NK1
                  IF (WI.LE.1E-30) GO TO 290
                  PIV    = H(1)
C  CALCULATE THE PARAMETERS OF THE GIVENS TRANSFORMATION.
                  CALL COSSIN(PIV,WI,DPRIME(J),COS,SIN)
C  TRANSFORMATIONS TO RIGHT HAND SIDE.
                  CALL ROTATE(PIV,COS,SIN,YI,C(J))
                  IF (J.EQ.NK1) GO TO 290
                  I2     = K1
                  IF (J.GT.N8) I2     = NK1 - J
                  DO 270 I = 1,I2
C  TRANSFORMATIONS TO LEFT HAND SIDE.
                      CALL ROTATE(PIV,COS,SIN,H(I+1),G(J,I))
                      H(I)   = H(I+1)
  270             CONTINUE
                  H(I2+1) = 0.
  280         CONTINUE
  290     CONTINUE
C  BACKWARD SUBSTITUTION TO OBTAIN THE B-SPLINE COEFFICIENTS.
          CALL BACK(G,C,NK1,K1,C)
C  COMPUTATION OF F(P).
          FP     = 0.
          L      = K2
          DO 320 IT = 1,M
              IF (X(IT).LT.T(L) .OR. L.GT.NK1) GO TO 300
              L      = L + 1
  300         L0     = L - K2
              TERM   = 0.
              DO 310 J = 1,K1
                  L0     = L0 + 1
                  TERM   = TERM + C(L0)*Q(IT,J)
  310         CONTINUE
              FP     = FP + W(IT)* (TERM-Y(IT))**2
  320     CONTINUE
C  TEST WHETHER THE APPROXIMATION SP(X) IS AN ACCEPTABLE SOLUTION.
          FPMS   = FP - S
C          TYPE*,ITER,FPMS,ACC
          IF (ABS(FPMS).LT.ACC) GO TO 400
C  TEST WHETHER THE MAXIMAL NUMBER OF ITERATIONS IS REACHED.
          IF (ITER.EQ.MAXIT) GO TO 360
C  CARRY OUT ONE MORE STEP OF THE ITERATION PROCESS.
          P2     = P
          F2     = FPMS
          IF (ICHECK.NE.0) GO TO 340
          IF ((F2-F3).GT.ACC) GO TO 330
C  OUR INITIAL CHOICE OF P IS TOO LARGE.
          P      = P*0.1E-02
          P3     = P2
          F3     = F2
          GO TO 350

  330     IF ((F1-F2).GT.ACC) GO TO 340
C  OUR INITIAL CHOICE OF P IS TOO SMALL
          P      = P*0.1E+04
          P1     = P2
          F1     = F2
          GO TO 350
C  TEST WETHER THE ITERATION PROCESS PROCEEDS AS THEORETICALLY
C  EXPECTED.
  340     IF (F2.GE.F1 .OR. F2.LE.F3) THEN
              TOL    = TOL*4.
              IF (TOL.GT.1.5) GO TO 370
              GO TO 10
          END IF
C 
          ICHECK = 1
C  FIND THE NEW VALUE FOR P.
          P = RATION(P1,F1,P2,F2,P3,F3)
C 
  350 CONTINUE
C  ERROR CODES AND MESSAGES.
  360 IER    = 3
      GO TO 400

  370 IER    = 2
      GO TO 400

  380 IER    = 1
      GO TO 400

  390 IER    = -1
  400 RETURN
C 
C 
      END

      SUBROUTINE BSPLIN(T,N,K,X,L,H)
C  SUBROUTINE BSPLIN EVALUATES THE (K+1) NON-ZERO B-SPLINES OF
C  DEGREE K AT T(L) <= X < T(L+1) USING THE STABLE RECURRENCE
C  RELATION OF DE BOOR AND COX.
C  THE DIMENSION SPECIFICATIONS OF THE FOLLOWING ARRAYS MUST BE
C  AT LEAST H(K+1),HH(K).

      IMPLICIT NONE
      INTEGER N,J,K,I,L,LJ,LI
      REAL T(N),H(6),HH(5),X,F

      H(1)   = 1.
      DO 30 J = 1,K
          DO 10 I = 1,J
              HH(I)  = H(I)
   10     CONTINUE
          H(1)   = 0.
          DO 20 I = 1,J
              LI     = L + I
              LJ     = LI - J
              F      = HH(I)/ (T(LI)-T(LJ))
              H(I)   = H(I) + F* (T(LI)-X)
              H(I+1) = F* (X-T(LJ))
   20     CONTINUE
   30 CONTINUE
      RETURN

      END

      SUBROUTINE COSSIN(PIV,WI,WW,COS,SIN)
C  SUBROUTINE COSSIN CALCULATES THE PARAMETERS OF A GIVENS
C  TRANSFORMATION WITHOUT SQUARE ROOTS.
      IMPLICIT NONE
      REAL COS,SIN,PIV,WI,WW,STORE,DD
      STORE  = PIV*WI
      DD     = WW + STORE*PIV
      COS    = WW/DD
      SIN    = STORE/DD
      WW     = DD
      WI     = COS*WI
      RETURN

      END

      SUBROUTINE ROTATE(PIV,COS,SIN,A,B)
C  SUBROUTINE ROTATE APPLIES A GIVENS ROTATION TO A AND B.
      IMPLICIT NONE
      REAL PIV,COS,SIN,A,B,STORE

      STORE  = B
      B      = COS*STORE + SIN*A
      A      = A - PIV*STORE
      RETURN

      END

      SUBROUTINE BACK(A,Z,N,K,C)
      IMPLICIT NONE
      INTEGER MXP
      PARAMETER (MXP=1000)
      INTEGER I,J,I1,M,L,N,K
      REAL A,Z,C,STORE
C  SUBROUTINE BACK CALCULATES THE SOLUTION OF THE SYSTEM OF
C  EQUATIONS A*C = Z WITH A A N X N UNIT UPPER TRIANGULAR MATRIX
C  OF BANDWIDTH K+1.
C  ATTENTION: THE FIRST DIMENSION SPECIFICATION OF MATRIX A MUST
C  BE THE SAME AS IN THE CALLING PROGRAM.
      DIMENSION A(MXP,K),Z(N),C(N)

      C(N)   = Z(N)
      I      = N - 1
      IF (I.EQ.0) GO TO 30
      DO 20 J = 2,N
          STORE  = Z(I)
          I1     = K
          IF (J.LE.K) I1     = J - 1
          M      = I
          DO 10 L = 1,I1
              M      = M + 1
              STORE  = STORE - C(M)*A(I,L)
   10     CONTINUE
          C(I)   = STORE
          I      = I - 1
   20 CONTINUE
   30 RETURN

      END

      SUBROUTINE NKNOT(X,M,T,N,FPINT,NRDATA,NRINT)
      IMPLICIT NONE
      INTEGER MXP
      PARAMETER (MXP=1000)
      INTEGER K,J,M,N,NRINT,JBEGIN,JPOINT,NUMBER,NRDATA(MXP)
      INTEGER MAXPT,MAXBEG,NEXT,JJ,JK,IHALF,NRX
      REAL X(M),T(MXP),FPINT(MXP),FPMAX
C  SUBROUTINE NKNOT LOCATES AN ADDITIONAL KNOT FOR A SPLINE OF DEGREE
C  K AND ADJUSTS THE CORRESPONDING PARAMETERS,I.E.
C    T     : THE POSITION OF THE KNOTS.
C    N     : THE NUMBER OF KNOTS.
C    NRINT : THE NUMBER OF KNOTINTERVALS.
C    FPINT : THE SUM OF SQUARES OF RESIDUAL RIGHT HAND SIDES
C            FOR EACH KNOT INTERVAL.
C    NRDATA: THE NUMBER OF DATA POINTS INSIDE EACH KNOT INTERVAL.
C  THE ARRAYS T,FPINT AND NRDATA MUST HAVE THE SAME DIMENSION
C  SPECIFICATIONS AS IN THE CALLING PROGRAM.

      K      = (N-NRINT-1)/2
C  SEARCH FOR KNOT INTERVAL T(NUMBER+K) <= X <= T(NUMBER+K+1) WHERE
C  FPINT(NUMBER) IS MAXIMAL ON THE CONDITION THAT NRDATA(NUMBER)
C  NOT EQUALS ZERO.

      FPMAX  = 0.
      MAXPT = 0
      MAXPBEG = 0
      NUMBER = 0
      JBEGIN = 1
      DO 20 J = 1,NRINT
          JPOINT = NRDATA(J)
          IF (FPMAX.GE.FPINT(J) .OR. JPOINT.EQ.0) GO TO 10
          FPMAX  = FPINT(J)
          NUMBER = J
          MAXPT  = JPOINT
          MAXBEG = JBEGIN
   10     JBEGIN = JBEGIN + JPOINT + 1
   20 CONTINUE
C  LET COINCIDE THE NEW KNOT T(NUMBER+K+1) WITH A DATA POINT X(NRX)
C  INSIDE THE OLD KNOT INTERVAL T(NUMBER+K) <= X <= T(NUMBER+K+1).
      IHALF  = MAXPT/2 + 1
      NRX    = MAXBEG + IHALF
      NEXT   = NUMBER + 1
      IF (NEXT.GT.NRINT) GO TO 40
C  ADJUSTS THE DIFFERENT PARAMETERS.
      DO 30 J = NEXT,NRINT
          JJ     = NEXT + NRINT - J
          FPINT(JJ+1) = FPINT(JJ)
          NRDATA(JJ+1) = NRDATA(JJ)
          JK     = JJ + K
          T(JK+1) = T(JK)
   30 CONTINUE
   40 NRDATA(NUMBER) = IHALF - 1
      NRDATA(NEXT) = MAXPT - IHALF
      FPINT(NUMBER) = FPMAX*FLOAT(NRDATA(NUMBER))/FLOAT(MAXPT)
      FPINT(NEXT) = FPMAX*FLOAT(NRDATA(NEXT))/FLOAT(MAXPT)
      JK     = NEXT + K
      T(JK)  = X(NRX)
      N      = N + 1
      NRINT  = NRINT + 1
      RETURN

      END

      SUBROUTINE DISCO(T,N,K2,B)
      IMPLICIT NONE
      INTEGER MXP
      PARAMETER (MXP=1000)
      INTEGER N,K1,K2,K,L,NK1,J,LP,LMK,IK,LJ,LK,JK,I
      REAL T(N),B(MXP,K2),H(12),PROD
C  SUBROUTINE DISCO CALCULATES THE DISCONTINUITY JUMPS OF THE KTH
C  DERIVATIVE OF THE B-SPLINES OF DEGREE K AT THE KNOTS T(K+2)..T(N-K-1)
C  THE FIRST DIMENSION SPECIFICATION OF THE MATRIX B MUST BE THE SAME AS
C  IN THE CALLING PROGRAM; H MUST HAVE A DIMENSION SPECIFICATION AT
C  LEAST 2*K+2.

      K1     = K2 - 1
      K      = K1 - 1
      NK1    = N - K1
      DO 40 L = K2,NK1
          LMK    = L - K1
          DO 10 J = 1,K1
              IK     = J + K1
              LJ     = L + J
              LK     = LJ - K2
              H(J)   = T(L) - T(LK)
              H(IK)  = T(L) - T(LJ)
   10     CONTINUE
          LP     = LMK
          DO 30 J = 1,K2
              JK     = J + K
              PROD   = 1.
              DO 20 I = J,JK
                  PROD   = PROD*H(I)
   20         CONTINUE
              LK     = LP + K1
              B(LMK,J) = (T(LK)-T(LP))/PROD
              LP     = LP + 1
   30     CONTINUE
   40 CONTINUE
      RETURN

      END

      REAL FUNCTION RATION(P1,F1,P2,F2,P3,F3)
      IMPLICIT NONE
      REAL P,P1,F1,P2,F2,P3,F3,H1,H2,H3
C  GIVEN THREE POINTS (P1,F1),(P2,F2) AND (P3,F3), FUNCTION RATION
C  GIVES THE VALUE OF P SUCH THAT THE RATIONAL INTERPOLATING FUNCTION
C  OF THE FORM R(P) = (U*P+V)/(P+W) EQUALS ZERO AT P.
      IF (P3.GT.0.) GO TO 10
C  VALUE OF P IN CASE P3 = INFINITY.
      P      = (P1* (F1-F3)*F2-P2* (F2-F3)*F1)/ ((F1-F2)*F3)
      GO TO 20
C  VALUE OF P IN CASE P3 ^= INFINITY.
   10 H1     = F1* (F2-F3)
      H2     = F2* (F3-F1)
      H3     = F3* (F1-F2)
      P      = - (P1*P2*H3+P2*P3*H1+P3*P1*H2)/ (P1*H1+P2*H2+P3*H3)
C  ADJUST THE VALUE OF P1,F1,P3 AND F3 SUCH THAT F1 > 0 AND F3 < 0.
   20 IF (F2.LT.0.) GO TO 30
      P1     = P2
      F1     = F2
      GO TO 40

   30 P3     = P2
      F3     = F2
   40 RATION = P
      RETURN

      END


      REAL FUNCTION DERIV(T,N,C,NK1,NU,ARG,L)
C  FUNCTION DERIV EVALUATES A SPLINE S(X) OF DEGREE K WHICH IS
C  GIVEN IN ITS NORMALIZED B-SPLINE REPRESENTATION OR CALCULATES
      IMPLICIT NONE
C  DERIVATIVES OF ANY SPECIFIED ORDER NU.
C
C  CALLING SEQUENCE
C     VALUE = DERIV(T,N,C,NK1,NU,ARG,L)
C
C  INPUT PARAMETERS:
C    T    : ARRAY,MINIMUM LENGTH N, WHICH CONTAINS THE POSITION
C           OF THE KNOTS OF S(X),I.E. THE POSITION OF THE INTERIOR
C           KNOTS T(K+2),...T(N-K-1) AS WELL AS THE POSITION OF THE
C           KNOTS T(1),...T(K+1) AND T(N-K),...T(N) WHICH ARE NEEDED
C           FOR THE B-SPLINE REPRESENTATION.
C    N    : INTEGER VALUE GIVING THE TOTAL NUMBER OF KNOTS OF S(X).
C    C    : ARRAY,LENGTH NK1, CONTAINING THE B-SPLINE COEFFICIENTS.
C    NK1  : INTEGER VALUE,GIVING THE DIMENSION OF S(X),I.E.NK1 = N-K-1.
C    NU   : INTEGER VALUE WHICH GIVES THE ORDER OF THE DERIVATIVE.
C    ARG  : REAL VALUE,GIVING THE VALUE OF THE ARGUMENT.
C    L    : INTEGER VALUE WHICH SPECIFIES THE POSITION OF THE ARGUMENT
C           I.E.   T(L) <= ARG <  T(L+1)   OR
C                  L = NK1 IF ARG = T(NK1+1).
C
C  OUTPUT PARAMETER:
C    VALUE: REAL VALUE,GIVING THE VALUE OF THE NUTH DERIVATIVE OF
C           S(X) AT X = ARG.
C
C  OTHER SUBROUTINES REQUIRED: NONE.
C
C  RESTRICTIONS:
C    1) NU >= 0
C    2) T(K+1) <= ARG <= T(NK1+1)
C  THE DIMENSION SPECIFICATION OF THE ARRAY H MUST BE AT LEAST K+1.
      INTEGER N,NK1,NU,L,I,J,K1,IK,JJ,NU1,NU2,LI,LJ
      REAL T(N),C(NK1),H(6),ARG

      DERIV  = 0.
      K1     = N - NK1
      IF (NU.LT.0 .OR. NU.GE.K1) RETURN
      DO 10 I = 1,K1
          IK     = L + I - K1
          H(I)   = C(IK)
   10 CONTINUE
      IF (NU.EQ.0) GO TO 40
      NU1    = NU + 1
      DO 30 J = 2,NU1
          DO 20 JJ = J,K1
              I      = J + K1 - JJ
              LI     = L + I - K1
              LJ     = L + I - J + 1
              H(I)   = (H(I)-H(I-1))/ (T(LJ)-T(LI))
   20     CONTINUE
   30 CONTINUE
      IF (NU.EQ.K1-1) GO TO 70
   40 NU2    = NU + 2
      DO 60 J = NU2,K1
          DO 50 JJ = J,K1
              I      = J + K1 - JJ
              LI     = L + I - K1
              LJ     = L + I - J + 1
              H(I)   = ((ARG-T(LI))*H(I)+ (T(LJ)-ARG)*H(I-1))/
     +                 (T(LJ)-T(LI))
   50     CONTINUE
   60 CONTINUE
   70 DERIV  = H(K1)
      IF (NU.EQ.0) RETURN
      DO 80 I = 1,NU
          DERIV  = DERIV*FLOAT(K1-I)
   80 CONTINUE
      RETURN

      END
 
      SUBROUTINE IERSPL (IER,LEVEL)
C
C. IDENTIFICATION
C
C     Author : P. BALLESTER    ESO-Garching
C
C     Versions :      1.00     27.06.90
C
C. PURPOSE
C
C     This subroutine displays error messages using MIDAS interface
C        STTPUT, depending on the error code returned by SMOSPL (IER)
C        and a level (LEVEL)
C
C. INPUTS :
C
C     IER   : Error code returned by the subroutine SMOSPL
C
C     LEVEL : = 0     No messages
C             = 1     Only severe error messages
C             = 2     Severe error messages and diagnosis
C             = 3     Maximum of information
C
C --- 
      IMPLICIT NONE
C
      INTEGER IER,LEVEL,STAT
      CHARACTER*80 STRING
C
      IF (IER.GT.0) THEN

         IF (LEVEL.GT.0) THEN
           STRING = 'Error while interpolating:'
           CALL STTPUT (STRING,STAT)
         ENDIF

         IF (LEVEL.GT.1) THEN

         IF (IER.EQ.1) THEN
            STRING = '--> Insufficient storage space'
            CALL STTPUT (STRING,STAT)
         ELSE IF (IER.EQ.2) THEN
            STRING = '--> Theoretically impossible result'
            CALL STTPUT (STRING,STAT)
            STRING = 
     + '    The requested relative accuracy is probably too small'
            CALL STTPUT (STRING,STAT)
         ELSE IF (IER.EQ.3) THEN
            STRING = '--> Maximal number of iterations has been reached'
            CALL STTPUT (STRING,STAT)
         ELSE IF (IER.EQ.10) THEN
            STRING = '--> Invalid input data'
            CALL STTPUT (STRING,STAT)
         ENDIF

         ENDIF

       ENDIF

       IF (IER.LE.0.AND.LEVEL.GT.2) THEN

          STRING = 'Warning : Normal return from spline interpolation'
          CALL STTPUT (STRING,STAT)

          IF (IER.EQ.-1) THEN
             STRING = 'S(x) is an interpolating spline'
             CALL STTPUT (STRING,STAT)
          ELSE IF (IER.EQ.-2) THEN
             STRING = 'S(x) is a polynomial of degree k'
             CALL STTPUT (STRING,STAT)
          ENDIF

       ENDIF
       RETURN
       END




       SUBROUTINE SORSPL (TAB,AUX1,AUX2,NP,SENSE)

C
C. IDENTIFICATION
C
C      Author : P. BALLESTER      ESO-Garching
C
C      Versions :        1.00       27.06.90
C                        1.10       17.07.90
C
C. PURPOSE
C
C      Sort values by increasing values in a 1D real array (TAB).
C      Elements of the auxiliary arrays AUX1 and AUX2 are also moved.
C      This subroutine is intended to be used before the SMOSPL routine
C      to sort the data, if necessary. The syntax is :
C                  CALL SORSPL (X,Y,W,NP,SENS)
C
C. INPUT/OUTPUT :
C
C      TAB     Real array   Address of the reference array
C      AUX1    Real array   Address of the first auxiliary array
C      AUX2    Real array   Address of the second auxiliary array
C      NP      Integer      Number of significant values in TAB,AUX1,AUX2
C      SENSE   Integer      >=0 = sort by increasing values
C                           <0 = sort by decreasing values
C
C. VERSIONS
C
C    1.00    Sort reference and auxiliary arrays by increasing values
C    1.10    Sort by increasing or decreasing values.
C
C. VERSION
C 051017	last modif
C -------------------------------------------------------------------------
C
       IMPLICIT NONE
C
       INTEGER FLAG,NP,PTR,SENSE
       REAL    TAB(NP),AUX1(NP),AUX2(NP),BUFFER,SS

       IF (SENSE.GE.0) SS =  1.0	! Sort by increasing values
       IF (SENSE.LT.0) SS = -1.0	! Sort by decreasing values

100    FLAG = 0
       DO 200 PTR = 1,NP-1

          IF ((SS*TAB(PTR)).GT.(SS*TAB(PTR+1))) THEN

              FLAG = 1

              BUFFER = TAB (PTR)
              TAB (PTR) = TAB (PTR+1)
              TAB (PTR+1) = BUFFER

              BUFFER = AUX1 (PTR)
              AUX1 (PTR) = AUX1 (PTR+1)
              AUX1 (PTR+1) = BUFFER

              BUFFER = AUX2 (PTR)
              AUX2 (PTR) = AUX2 (PTR+1)
              AUX2 (PTR+1) = BUFFER

          ENDIF
200    CONTINUE

       IF (FLAG.NE.0) GOTO 100

       RETURN
       END         




         SUBROUTINE ARGSPL (ARG,KNOTS,NBKNOT,DEGREE,PREPOS,POSARG,NK1)

         IMPLICIT NONE
         INTEGER  NBKNOT,DEGREE,POSARG
         REAL     ARG,KNOTS(NBKNOT)
         INTEGER   NK1,PREPOS,STAT
         CHARACTER STRING*80

         NK1 = NBKNOT - DEGREE - 1

         IF (ARG.EQ.KNOTS(NK1+1)) THEN
             POSARG = NK1
             GOTO 550
         ENDIF

         DO 500 POSARG = PREPOS , NK1
            IF (KNOTS(POSARG).LE.ARG.AND.ARG.LT.KNOTS(POSARG+1)) THEN
cd               WRITE (STRING,5500) POSARG,KNOTS(POSARG),ARG,KNOTS(POSARG+1)
cd               CALL STTPUT (STRING,STAT)
                 PREPOS = POSARG  
                GOTO 550
            ENDIF
500      CONTINUE

C
C --- 5.1.2. Cannot find argument position
C
         STRING = '*** Argument position not defined'
         CALL STTPUT (STRING,STAT)
         WRITE (STRING,5500) POSARG,KNOTS(POSARG),ARG,KNOTS(POSARG+1)
5500     FORMAT ('Pos : ',I5,' Relation : ',F6.1,' <',F6.1,' <',F6.1)
         CALL STTPUT (STRING,STAT)

550      RETURN
         END
