C SUBROUTINES OF THE SOFTWARE PACKAGE 'FITPACK' BY A.K. CLINE
C USED TO SPECIFY THE MODEL FOR THE COMPLETE RAY TRACING ALGORITHM.
C
C THIS FILE CONSISTS OF THE FOLLOWING PARTS:
C     (0) AUXILIARY SUBROUTINE
C             SNHCSH
C         COMMON TO ALL THE FOLLOWING PARTS.
C     (1) THE SUBROUTINES PREPARING THE PARAMETERS NECESSARY TO COMPUTE
C         AN INTERPOLATORY FUNCTION:
C             CURVN1 (HERMITE REPRESENTATION OF 1-D FUNCTION),
C             CURVB1 (B-SPLINE REPRESENTATION OF 1-D FUNCTION),
C             SURFB1 (B-SPLINE REPRESENTATION OF 2-D FUNCTION),
C             VAL3B1 (B-SPLINE REPRESENTATION OF 3-D FUNCTION),
C             VGEN (AUXILIARY SUBROUTINE),
C             TERMS (AUXILIARY SUBROUTINE),
C             TRIDEC (AUXILIARY SUBROUTINE),
C             TRISOL (AUXILIARY SUBROUTINE).
C         SUBROUTINES CURVN1 AND CURVB1 ARE ALTERNATIVES.
C     (2) THE SUBROUTINES EVALUATING THE VALUE, FIRST AND SECOND PARTIAL
C         DERIVATIVES OF THE INTERPOLATORY FUNCTION AT A GIVEN POINT:
C             CURV2D (HERMITE REPRESENTATION OF 1-D FUNCTION),
C             CURVBD (B-SPLINE REPRESENTATION OF 1-D FUNCTION),
C             SURFBD (B-SPLINE REPRESENTATION OF 2-D FUNCTION),
C             VAL3BD (B-SPLINE REPRESENTATION OF 3-D FUNCTION),
C             DSPLNZ (AUXILIARY SUBROUTINE),
C             INTRVL (AUXILIARY EXTERNAL FUNCTION).
C         SUBROUTINES CURV2D AND CURVBD ARE ALTERNATIVES.
C
C TAKEN FROM:
C     FITPACK - A SOFTWARE PACKAGE FOR CURVE AND SURFACE FITTING
C     EMPLOYING SPLINES UNDER TENSION
C     BY ALAN KAYLOR CLINE, DEPARTMENT OF COMPUTER SCIENCES,
C     THE UNIVERSITY OF TEXAS AT AUSTIN, AUGUST 31, 1981.
C NOTE 1:
C     TO CONFORM WITH THE FORTRAN77 STANDARD, DUMMY ARRAY DIMENSIONS (1)
C     HAVE BEEN CHANGED TO (*), AND SUBROUTINE TRISOL HAS BEEN REVISED.
C NOTE 2:
C     SUBROUTINES CURVB1 AND CURVBD DO NOT BELONG TO THE ORIGINAL
C     VERSION OF FITPACK.
C NOTE 3:
C     TO GET THE ORIGINAL VERSIONS OF THE SUBROUTINES SURFBD AND VAL3BD,
C     THE STATEMENT WITH 'CALL VAR2' MUST BE REMOVED FROM EACH OF THEM.
C     THE STATEMENTS HAVE BEEN ADDED BY L.KLIMES FOR THE SAKE OF INVERSE
C     MODELLING TO THE SUBROUTINES CURVBD, SURFBD, AND VAL3BD.
C     THE THREE APPEARENCES OF THE STATEMENTS 'CALL VAR2' ARE DENOTED BY
C     'V' IN THE FIRST COLUMN.  THE THREE LINES SHOULD BE REMOVED OR
C     MODIFIED BEFORE COMPILATION.
C
C=======================================================================
C
C PART 0:
C
C=======================================================================
C
      SUBROUTINE SNHCSH (SINHM,COSHM,X,ISW)
C
      INTEGER ISW
      REAL SINHM,COSHM,X
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                       CODED BY A. K. CLINE AND R. J. RENKA
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE RETURNS APPROXIMATIONS TO
C       SINHM(X) = SINH(X)-X
C       COSHM(X) = COSH(X)-1
C AND
C       COSHMM(X) = COSH(X)-1-X*X/2
C WITH RELATIVE ERROR LESS THAN 6.16E-6
C
C ON INPUT--
C
C   X CONTAINS THE VALUE OF THE INDEPENDENT VARIABLE.
C
C   ISW INDICATES THE FUNCTION DESIRED
C           = -1 IF ONLY SINHM IS DESIRED,
C           =  0 IF BOTH SINHM AND COSHM ARE DESIRED,
C           =  1 IF ONLY COSHM IS DESIRED,
C           =  2 IF ONLY COSHMM IS DESIRED,
C           =  3 IF BOTH SINHM AND COSHMM ARE DESIRED.
C
C ON OUTPUT--
C
C   SINHM CONTAINS THE VALUE OF SINHM(X) IF ISW .LE. 0 OR
C   ISW .EQ. 3 (SINHM IS UNALTERED IF ISW .EQ.1 OR ISW .EQ.
C   2).
C
C   COSHM CONTAINS THE VALUE OF COSHM(X) IF ISW .EQ. 0 OR
C   ISW .EQ. 1 AND CONTAINS THE VALUE OF COSHMM(X) IF ISW
C   .GE. 2 (COSHM IS UNALTERED IF ISW .EQ. -1).
C
C AND
C
C   X AND ISW ARE UNALTERED.
C
C-----------------------------------------------------------
C
      DATA SP2/5.04850926418006E-04/,
     *     SP1/3.62841692246321E-02/,
     *     SQ1/-1.37157937097122E-02/
      DATA CP2/1.31625490355985E-03/,
     *     CP1/6.57866547762733E-02/,
     *     CQ1/-1.75465241841312E-02/
      DATA ZP2/1.40048186158693E-04/,
     *     ZP1/1.67309141907440E-02/,
     *     ZQ2/9.82154460147143E-05/,
     *     ZQ1/-1.66024148976133E-02/
      XX = X
      AX = ABS(XX)
      XS = XX*XX
      IF ((AX .GE. 2.20) .OR. (AX .GE. 1.17 .AND.
     *     ISW .NE. 2)) EXPX = EXP(AX)
C
C SINHM APPROXIMATION
C
      IF (ISW .EQ. 1 .OR. ISW .EQ. 2) GO TO 2
      IF (AX .GE. 1.17) GO TO 1
      SINHM = (((SP2*XS+SP1)*XS+1.)*XS*XX)/((SQ1*XS+1.)*6.)
      GO TO 2
    1 SINHM = (EXPX-1./EXPX)/2.-AX
      IF (XX .LT. 0.) SINHM = -SINHM
C
C COSHM APPROXIMATION
C
    2 IF (ISW .NE. 0 .AND. ISW .NE. 1) GO TO 4
      IF (AX .GE. 1.17) GO TO 3
      COSHM = (((CP2*XS+CP1)*XS+1.)*XS)/((CQ1*XS+1.)*2.)
      GO TO 4
    3 COSHM = (EXPX+1./EXPX)/2.-1.
C
C COSHMM APPROXIMATION
C
    4 IF (ISW .LE. 1) RETURN
      IF (AX .GE. 2.20) GO TO 5
      COSHM = (((ZP2*XS+ZP1)*XS+1.)*XS*XS)/(((ZQ2*XS+ZQ1)*XS
     *        +1.)*24.)
      RETURN
    5 COSHM = (EXPX+1./EXPX)/2.-1.-XS/2.
      RETURN
      END
C
C=======================================================================
C
C PART 1:
C
C=======================================================================
C
      SUBROUTINE CURVN1 (N,X,Y,YP,TEMP,SIGMA,IERR)
C
      INTEGER N,IERR
      REAL X(N),Y(N),YP(N),TEMP(N),SIGMA
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                    CODED BY A. K. CLINE AND S. E. GALINSKY
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO
C COMPUTE A NATURAL INTERPOLATORY SPLINE UNDER TENSION
C THROUGH A SEQUENCE OF FUNCTIONAL VALUES. FOR ACTUAL
C COMPUTATION OF POINTS ON THE CURVE IT IS NECESSARY TO CALL
C THE FUNCTION CURV2.
C
C ON INPUT--
C
C   N IS THE NUMBER OF VALUES TO BE INTERPOLATED (N.GE.2).
C
C   X IS AN ARRAY OF THE N INCREASING ABSCISSAE OF THE
C   FUNCTIONAL VALUES.
C
C   Y IS AN ARRAY OF THE N ORDINATES OF THE VALUES, (I. E.
C   Y(K) IS THE FUNCTIONAL VALUE CORRESPONDING TO X(K) ).
C
C   YP IS AN ARRAY OF LENGTH AT LEAST N.
C
C   TEMP IS AN ARRAY OF LENGTH AT LEAST N WHICH IS USED FOR
C   SCRATCH STORAGE.
C
C AND
C
C   SIGMA CONTAINS THE TENSION FACTOR. THIS VALUE INDICATES
C   THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY ZERO
C   (E.G. .001) THE RESULTING CURVE IS APPROXIMATELY A
C   CUBIC SPLINE. IF ABS(SIGMA) IS LARGE (E.G. 50.) THE
C   RESULTING CURVE IS NEARLY A POLYGONAL LINE. IF SIGMA
C   EQUALS ZERO A CUBIC SPLINE RESULTS.  A STANDARD VALUE
C   FOR SIGMA IS APPROXIMATELY 1. IN ABSOLUTE VALUE.
C
C ON OUTPUT--
C
C   YP CONTAINS THE VALUES OF THE SECOND DERIVATIVE OF THE
C   CURVE AT THE GIVEN NODES.
C
C   IERR CONTAINS AN ERROR FLAG,
C        = 0 FOR NORMAL RETURN,
C        = 1 IF N IS LESS THAN 2,
C        = 2 IF X-VALUES ARE NOT STRICTLY INCREASING.
C
C AND
C
C   N, X, Y, AND SIGMA ARE UNALTERED.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES SNHCSH.
C
C-----------------------------------------------------------
C
      NM1 = N-1
      NP1 = N+1
      IERR = 0
      IF (N .LE. 1) GO TO 4
      IF (X(N) .LE. X(1)) GO TO 5
C
C DENORMALIZE TENSION FACTOR
C
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
C
C SET UP RIGHT HAND SIDE AND TRIDIAGONAL SYSTEM FOR YP AND
C PERFORM FORWARD ELIMINATION
C
      DELX1 = X(2)-X(1)
      IF (DELX1 .LE. 0.) GO TO 5
      DX1 = (Y(2)-Y(1))/DELX1
      CALL TERMS (DIAG1,SDIAG1,SIGMAP,DELX1)
      SDIAG1 = 0.
      YP(1) = 0.
      TEMP(1) = 0.
      IF (N .EQ. 2) GO TO 2
      DO 1 I = 2,NM1
        DELX2 = X(I+1)-X(I)
        IF (DELX2 .LE. 0.) GO TO 5
        DX2 = (Y(I+1)-Y(I))/DELX2
        CALL TERMS (DIAG2,SDIAG2,SIGMAP,DELX2)
        DIAG = DIAG1+DIAG2-SDIAG1*TEMP(I-1)
        YP(I) = (DX2-DX1-SDIAG1*YP(I-1))/DIAG
        TEMP(I) = SDIAG2/DIAG
        DX1 = DX2
        DIAG1 = DIAG2
    1   SDIAG1 = SDIAG2
   2  YP(N) = 0.
      TEMP(N-1) = 0.
C
C PERFORM BACK SUBSTITUTION
C
      DO 3 I = 2,N
        IBAK = NP1-I
    3   YP(IBAK) = YP(IBAK)-TEMP(IBAK)*YP(IBAK+1)
      RETURN
C
C TOO FEW POINTS
C
    4 IERR = 1
      RETURN
C
C X-VALUES NOT STRICTLY INCREASING
C
    5 IERR = 2
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE CURVB1 (NX,X,W,C,VX,TEMP,SIGMA,IERR)
C
      INTEGER NX,IERR
      REAL X(NX),W(NX),C(NX),VX(5,NX),TEMP(*),SIGMA
C
C                                      COMPLEMENT TO FITPACK
C                                       BY ALAN KAYLOR CLINE
C                                   CODED -- OCTOBER 9, 1986
C                                            BY LUDEK KLIMES
C                                      INST. GEOL. GEOTECHN.
C                               CZECHOSL. ACAD. SCI., PRAGUE
C
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO
C COMPUTE AN INTERPOLATORY FUNCTION ON A ONE DIMENSIONAL
C GRID. THE FUNCTION DETERMINED CAN BE
C REPRESENTED BY SPLINES UNDER TENSION.  FOR ACTUAL
C MAPPING OF POINTS IT IS NECESSARY TO CALL THE SUBROUTINE
C CURVBD, WHICH ALSO RETURNS FIRST AND SECOND DERIVATIVES.
C
C ON INPUT--
C
C   NX IS THE NUMBER OF GRID POINTS.
C   (NX SHOULD BE AT LEAST 2)
C
C   X IS ARRAY OF THE NX COORDINATES OF THE GRID POINTS.
C   THESE SHOULD BE STRICTLY INCREASING.
C
C   W IS AN ARRAY OF THE NX FUNCTIONAL VALUES AT THE
C   THE GRID POINTS, I. E. W(I,J) CONTAINS THE FUNCTIONAL
C   VALUE AT X(I) FOR I = 1,...,NX .
C
C   C IS AN ARRAY OF AT LEAST NX LOCATIONS. THIS
C   PARAMETER MAY COINCIDE WITH W IN WHICH CASE W IS
C   DESTROYED ON OUTPUT.
C
C   VX IS THE ARRAY OF AT LEAST 5 * NX  LOCATIONS.
C
C   TEMP IS AN ARRAY OF AT LEAST 3 * NX LOCATIONS
C   WHICH IS USED FOR SCRATCH STORAGE.
C
C   SIGMA CONTAINS THE TENSION FACTOR. THIS VALUE INDICATE
C   THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY ZERO
C   (E. G. .001) THE RESULTING SURFACE IS APPROXIMATELY THE
C   TENSOR PRODUCT OF CUBIC SPLINES. IF ABS(SIGMA) IS LARGE
C   (E. G. 50.) THE RESULTING SURFACE IS APPROXIMATELY
C   BI-LINEAR. IF SIGMA EQUALS ZERO TENSOR PRODUCTS OF CUBIC
C   SPLINES RESULT. A STANDARD VALUE FOR SIGMA IS
C   APPROXIMATELY 1. IN ABSOLUTE VALUE.
C
C ON OUTPUT--
C
C   C CONTAINS THE COEFFICIENTS OF A REPRESENTATION OF THE
C   INTERPOLATED FUNCTION IN A B-SPLINE FORM.
C
C   VX CONTAINS B-SPLINE UNDER TENSION BASIS DATA.
C
C   IERR CONTAINS AN ERROR FLAG.
C        = 0 FOR NORMAL RETURN,
C        = 1 IF NX IS LESS THAN 2,
C        = 2 IF THE X-ARRAY IS NOT STRICTLY
C            INCREASING.
C
C AND
C
C   NONE OF THE INPUT PARAMETERS ARE ALTERED (EXCEPT W IF
C   THIS PARAMETER AND C ARE IDENTICAL IN THE CALLING
C   SEQUENCE).
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES VGEN, TERMS,
C SNHCSH, TRIDEC, AND TRISOL.
C
C-----------------------------------------------------------
C
C COPY W INTO C
C
      DO 1 I = 1,NX
    1   C(I) = W(I)
C
C GENERATE BASIS FUNCTIONS ALONG X-GRID
C SET UP TRIDIAGONAL SYSTEM AND SOLVE
C
      CALL VGEN (NX,X,SIGMA,VX,IERR)
      IF (IERR .NE. 0) RETURN
      DO 2 I = 2,NX
    2   TEMP(I) = VX(5,I-1)
      NXPI = NX
      DO 3 I = 1,NX
        NXPI = NXPI+1
    3   TEMP(NXPI) = 1.
      DO 4 I = 2,NX
        NXPI = NXPI+1
    4   TEMP(NXPI) = VX(4,I)
      CALL TRIDEC (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),
     *             TEMP(1),TEMP(NX+1),IERR)
      CALL TRISOL (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),C,NX,
     *             1,1)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE SURFB1 (NX,NY,X,Y,W,NW1,C,VX,VY,TEMP,SIGMA,
     *                   IERR)
C
      INTEGER NX,NY,NW1,IERR
      REAL X(NX),Y(NY),W(NW1,NY),C(NX,NY),VX(5,NX),VY(5,NY),
     *     TEMP(*),SIGMA
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                                 CODED BY ALAN KAYLOR CLINE
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO
C COMPUTE AN INTERPOLATORY FUNCTION ON A TWO DIMENSIONAL
C RECTANGULAR GRID. THE FUNCTION DETERMINED CAN BE
C REPRESENTED AS A TENSOR PRODUCT OF SPLINES UNDER TENSION
C FOR ACTUAL MAPPING OF POINTS IT IS NECESSARY TO CALL THE
C SUBROUTINE SURFBD, WHICH ALSO RETURNS FIRST AND SECOND
C PARTIAL DERIVATIVES.
C
C ON INPUT--
C
C   NX AND NY ARE THE NUMBER OF GRID LINES IN THE X- AND Y
C   DIRECTIONS, RESPECTIVELY, OF THE RECTANGULAR GRID. (NX
C   AND NY SHOULD BE AT LEAST 2.)
C
C   X AND Y ARE ARRAYS OF THE NX AND NY COORDINATES OF THE
C   GRID LINES IN X- AND Y-DIRECTIONS, RESPECTIVELY. THESE
C   SHOULD BE STRICTLY INCREASING.
C
C   W IS AN ARRAY OF THE NX * NY FUNCTIONAL VALUES AT THE
C   THE GRID POINTS, I. E. W(I,J) CONTAINS THE FUNCTIONAL
C   VALUE AT (X(I),Y(J)) FOR I = 1,...,NX, AND J = 1,...,NY.
C
C   NW1 IS THE FIRST DIMENSION OF THE ARRAY W USED IN THE
C   CALLING PROGRAM (NW1 .GE. NX).
C
C   C IS AN ARRAY OF AT LEAST NX * NY LOCATIONS. THIS
C   PARAMETER MAY COINCIDE WITH W IN WHICH CASE W IS
C   DESTROYED ON OUTPUT.
C
C   VX AND VY ARE ARRAYS OF AT LEAST 5 * NX AND 5 * NY
C   LOCATIONS, RESPECTIVELY.
C
C   TEMP IS AN ARRAY OF AT LEAST 3 * MAX(NX, NY) LOCATIONS
C   WHICH IS USED FOR SCRATCH STORAGE.
C
C AND
C
C   SIGMA CONTAINS THE TENSION FACTOR. THIS VALUE INDICATE
C   THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY ZERO
C   (E. G. .001) THE RESULTING SURFACE IS APPROXIMATELY THE
C   TENSOR PRODUCT OF CUBIC SPLINES. IF ABS(SIGMA) IS LARGE
C   (E. G. 50.) THE RESULTING SURFACE IS APPROXIMATELY
C   BI-LINEAR. IF SIGMA EQUALS ZERO TENSOR PRODUCTS OF CUBIC
C   SPLINES RESULT. A STANDARD VALUE FOR SIGMA IS
C   APPROXIMATELY 1. IN ABSOLUTE VALUE.
C
C ON OUTPUT--
C
C   C CONTAINS THE COEFFICIENTS OF A REPRESENTATION OF THE
C   INTERPOLATED FUNCTION IN A B-SPLINE TENSOR PRODUCTION
C   FORM.
C
C   VX AND VY CONTAIN B-SPLINE UNDER TENSION BASIS DATA.
C
C   IERR CONTAINS AN ERROR FLAG.
C        = 0 FOR NORMAL RETURN,
C        = 1 IF NX OR NY IS LESS THAN 2,
C        = 2 IF THE X- OR Y-ARRAYS ARE NOT STRICTLY
C            INCREASING.
C
C AND
C
C   NONE OF THE INPUT PARAMETERS ARE ALTERED (EXCEPT W IF
C   THIS PARAMETER AND C ARE IDENTICAL IN THE CALLING
C   SEQUENCE).
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES VGEN, TERMS,
C SNHCSH, TRIDEC, AND TRISOL.
C
C--------------------------------------------------------- -
C
C COPY W INTO C
C
      DO 1 J = 1,NY
        DO 1 I = 1,NX
    1     C(I,J) = W(I,J)
C
C GENERATE BASIS FUNCTIONS ALONG X-GRID
C SET UP TRIDIAGONAL SYSTEM AND SOLVE
C
      CALL VGEN (NX,X,SIGMA,VX,IERR)
      IF (IERR .NE. 0) RETURN
      DO 2 I = 2,NX
    2   TEMP(I) = VX(5,I-1)
      NXPI = NX
      DO 3 I = 1,NX
        NXPI = NXPI+1
    3   TEMP(NXPI) = 1.
      DO 4 I = 2,NX
        NXPI = NXPI+1
    4   TEMP(NXPI) = VX(4,I)
      CALL TRIDEC (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),
     *             TEMP(1),TEMP(NX+1),IERR)
      CALL TRISOL (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),C,NX,
     *            NY,1)
C
C GENERATE BASIS FUNCTIONS ALONG Y-GRID
C SET UP TRIDIAGONAL SYSTEM AND SOLVE
C
      CALL VGEN (NY,Y,SIGMA,VY,IERR)
      IF (IERR .NE. 0) RETURN
      DO 5 J = 2,NY
    5   TEMP(J) = VY(5,J-1)
      NYPJ = NY
      DO 6 J = 1,NY
        NYPJ = NYPJ+1
    6   TEMP(NYPJ) = 1.
      DO 7 J = 2,NY
        NYPJ = NYPJ+1
    7   TEMP(NYPJ) = VY(4,J)
      CALL TRIDEC (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1),
     *             TEMP(1),TEMP(NY+1),IERR)
      CALL TRISOL (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1),C,1,
     *            NX,NX)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE VAL3B1 (NX,NY,NZ,X,Y,Z,W,NW1,NW2,C,VX,VY,
     *                   VZ,TEMP,SIGMA,IERR)
C
      INTEGER NX,NY,NZ,NW1,NW2,IERR
      REAL X(NX),Y(NY),Z(NZ),W(NW1,NW2,NZ),C(NX,NY,NZ),
     *     VX(5,NX),VY(5,NY),VZ(5,NZ),TEMP(*),SIGMA
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                                 CODED BY ALAN KAYLOR CLINE
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO
C COMPUTE AN INTERPOLATORY FUNCTION ON A THREE DIMENSIONAL
C RECTANGULAR GRID. THE FUNCTION DETERMINED CAN BE
C REPRESENTED AS A TENSOR PRODUCT OF SPLINES UNDER TENSION.
C FOR ACTUAL MAPPING OF POINTS IT IS NECESSARY TO CALL THE
C SUBROUTINE VAL3BD, WHICH ALSO RETURNS FIRST AND SECOND
C PARTIAL DERIVATIVES.
C
C ON INPUT--
C
C   NX, NY, AND NZ ARE THE NUMBER OF GRID LINES IN THE X-,
C   Y-, AND Z-DIRECTIONS, RESPECTIVELY, OF THE RECTANGULAR
C   GRID. (NX, NY, AND NZ SHOULD BE AT LEAST 2.)
C
C   X, Y, AND Z ARE ARRAYS OF THE NX, NY, AND NZ COORDINATES
C   OF THE GRID LINES IN THE X-, Y-, AND Z-DIRECTIONS,
C   RESPECTIVELY. THESE SHOULD BE STRICTLY INCREASING.
C
C   W IS AN ARRAY OF THE NX * NY * NZ FUNCTIONAL VALUES AT
C   THE GRID POINTS, I. E. W(I,J,K) CONTAINS THE FUNCTIONAL
C   VALUE AT (X(I),Y(J),Z(K)) FOR I = 1,...,NX,
C   J = 1,...,NY, AND K = 1,...,NZ.
C
C   NW1 AND NW2 ARE THE FIRST TWO DIMENSIONS OF THE ARRAY W
C   USED IN THE CALLING PROGRAM (NW1 .GE. NX AND NW2 .GE.
C   NY).
C
C   C IS AN ARRAY OF AT LEAST NX * NY * NZ LOCATIONS. THIS
C   PARAMETER MAY COINCIDE WITH W IN WHICH CASE W IS
C   DESTROYED ON OUTPUT.
C
C   VX, VY, AND VZ ARE ARRAYS OF AT LEAST 5 * NX, 5 * NY,
C   AND 5 * NZ LOCATIONS, RESPECTIVELY.
C
C   TEMP IS AN ARRAY OF AT LEAST 3 * MAX(NX, NY, NZ)
C   LOCATIONS WHICH IS USED FOR SCRATCH STORAGE.
C
C AND
C
C   SIGMA CONTAINS THE TENSION FACTOR. THIS VALUE INDICATES
C   THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY ZERO
C   (E. G. .001) THE RESULTING SURFACE IS APPROXIMATELY THE
C   TENSOR PRODUCT OF CUBIC SPLINES. IF ABS(SIGMA) IS LARGE
C   (E. G. 50.) THE RESULTING SURFACE IS APPROXIMATELY
C   TRI-LINEAR. IF SIGMA EQUALS ZERO TENSOR PRODUCTS OF
C   CUBIC SPLINES RESULT. A STANDARD VALUE FOR SIGMA IS
C   APPROXIMATELY 1. IN ABSOLUTE VALUE.
C
C ON OUTPUT--
C
C   C CONTAINS THE COEFFICIENTS OF A REPRESENTATION OF THE
C   INTERPOLATED FUNCTION IN A B-SPLINE TENSOR PRODUCTION
C   FORM.
C
C   VX, VY, AND VZ CONTAIN B-SPLINE UNDER TENSION BASIS
C   DATA.
C
C   IERR CONTAINS AN ERROR FLAG.
C        = 0 FOR NORMAL RETURN,
C        = 1 IF NX, NY, OR NZ IS LESS THAN 2,
C        = 2 IF THE X-, Y-, OR Z-ARRAYS ARE NOT STRICTLY
C            INCREASING.
C
C AND
C
C   NONE OF THE INPUT PARAMETERS ARE ALTERED (EXCEPT W IF
C   THIS PARAMETER AND C ARE IDENTICAL IN THE CALLING
C   SEQUENCE).
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES VGEN, TERMS,
C SNHCSH, TRIDEC, AND TRISOL.
C
C-----------------------------------------------------------
C
C COPY W INTO C
C
      DO 1 K = 1,NZ
        DO 1 J = 1,NY
          DO 1 I = 1,NX
    1       C(I,J,K) = W(I,J,K)
C
C GENERATE BASIS FUNCTIONS ALONG X-GRID
C SET UP TRIDIAGONAL SYSTEM AND SOLVE
C
      CALL VGEN (NX,X,SIGMA,VX,IERR)
      IF (IERR .NE. 0) RETURN
      DO 2 I = 2,NX
    2   TEMP(I) = VX(5,I-1)
      NXPI = NX
      DO 3 I = 1,NX
        NXPI = NXPI+1
    3   TEMP(NXPI) = 1.
      DO 4 I = 2,NX
        NXPI = NXPI+1
    4   TEMP(NXPI) = VX(4,I)
      CALL TRIDEC (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),
     *             TEMP(1),TEMP(NX+1),IERR)
      CALL TRISOL (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),C,NX,
     *            NY*NZ,1)
C
C GENERATE BASIS FUNCTIONS ALONG Y-GRID
C SET UP TRIDIAGONAL SYSTEM AND SOLVE
C
      CALL VGEN (NY,Y,SIGMA,VY,IERR)
      IF (IERR .NE. 0) RETURN
      DO 5 J = 2,NY
    5   TEMP(J) = VY(5,J-1)
      NYPJ = NY
      DO 6 J = 1,NY
        NYPJ = NYPJ+1
    6   TEMP(NYPJ) = 1.
      DO 7 J = 2,NY
        NYPJ = NYPJ+1
    7   TEMP(NYPJ) = VY(4,J)
      CALL TRIDEC (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1),
     *             TEMP(1),TEMP(NY+1),IERR)
      DO 8 K = 1,NZ
    8   CALL TRISOL (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1),C(1,1,K),
     *               1,NX,NX)
C
C GENERATE BASIS FUNCTIONS ALONG Z-GRID
C SET UP TRIDIAGONAL SYSTEM AND SOLVE
C
      CALL VGEN (NZ,Z,SIGMA,VZ,IERR)
      IF (IERR .NE. 0) RETURN
      DO 9 K = 2,NZ
    9   TEMP(K) = VZ(5,K-1)
      NZPK = NZ
      DO 10 K = 1,NZ
        NZPK = NZPK+1
   10   TEMP(NZPK) = 1.
      DO 11 K = 2,NZ
        NZPK = NZPK+1
   11   TEMP(NZPK) = VZ(4,K)
      CALL TRIDEC (NZ,TEMP(1),TEMP(NZ+1),TEMP(2*NZ+1),
     *             TEMP(1),TEMP(NZ+1),IERR)
      CALL TRISOL (NZ,TEMP(1),TEMP(NZ+1),TEMP(2*NZ+1),C,1,
     *            NX*NY,NX*NY)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE VGEN (N,X,SIGMA,V,IERR)
C
      INTEGER N,IERR
      REAL X(N),SIGMA,V(5,N)
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                                 CODED BY ALAN KAYLOR CLINE
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE GENERATES AN ARRAY OF COEFFICIENTS USED BY
C OTHER SUBROUTINES FOR THE DETERMINATION OF A B-SPLINE
C UNDER TENSION BASIS.
C
C ON INPUT--
C
C   N IS THE NUMBER OF KNOTS DEFINING THE BASIS (N .GE. 2).
C
C   X IS THE ARRAY OF THE N INCREASING KNOTS. ANY LINEAR
C   COMBINATION OF THE RESULTING BASIS WILL HAVE THIRD
C   DERIVATIVE DISCONTINUITIES ONLY AT THE INTERIOR KNOTS,
C   (I. E. X(2),...,X(N-1) ).
C
C   SIGMA CONTAINS THE TENSION FACTOR. THIS VALUE INDICATES
C   THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY ZERO
C   (E. G. .001) THE BASIS FUNCTIONS ARE APPROXIMATELY CUBIC
C   SPLINES. IF ABS(SIGMA) IS LARGE (E. G. 50.) THE BASIS
C   FUNCTIONS ARE NEARLY PIECEWISE LINEAR. IF SIGMA EQUALS
C   ZERO A CUBIC SPLINE BASIS RESULTS. A STANDARD VALUE FOR
C   SIGMA IS APPROXIMATELY 1. IN ABSOLUTE VALUE.
C
C AND
C
C   V IS AN ARRAY OF AT LEAST 5*N LOCATIONS.
C
C ON OUTPUT--
C
C   V CONTAINS CERTAIN COEFFICIENTS TO BE USED BY OTHER
C   SUBPROGRAMS FOR THE DETERMINATION OF THE B-SPLINE UNDER
C   TENSION BASIS. CONSIDERED AS A 5 BY N ARRAY, FOR I = 1,
C   ... , N, B-SPLINE BASIS FUNCTION I IS SPECIFIED BY--
C     V(1,I) = SECOND DERIVATIVE AT X(I-1), FOR I .NE. 1,
C     V(2,I) = SECOND DERIVATIVE AT X(I),   FOR ALL I,
C     V(3,I) = SECOND DERIVATIVE AT X(I+1), FOR I .NE. N,
C     V(4,I) = FUNCTION VALUE AT X(I-1),    FOR I .NE. 1,
C     V(5,I) = FUNCTION VALUE AT X(I+1),    FOR I .NE. N,
C   AND THE PROPERTIES THAT IT HAS--
C     1. FUNCTION VALUE 1 AT X(I),
C     2. FUNCTION VALUE AND SECOND DERIVATIVE = 0 AT
C        X(1), ... , X(I-2), AND X(I+2), ... , X(N).
C   IN V(5,N) AND V(3,N) ARE CONTAINED FUNCTION VALUE AND
C   SECOND DERIVATIVE OF BASIS FUNCTION ZERO AT X(1),
C   RESPECTIVELY. IN V(4,1) AND V(1,1) ARE CONTAINED
C   FUNCTION VALUE AND SECOND DERIVATIVE OF BASIS FUNCTION
C   N+1 AT X(N), RESPECTIVELY. FUCTION VALUE AND SECOND
C   DERIVATIVE OF THESE TWO BASIS FUNCTIONS ARE ZERO AT ALL
C   OTHER KNOTS. ONLY BASIS FUNCTION ZERO HAS NON-ZERO
C   SECOND DERIVATIVE VALUE AT X(1) AND ONLY BASIS
C   FUNCTION N+1 HAS NON-ZERO SECOND DERIVATIVE AT X(N).
C
C   IERR CONTAINS AN ERROR FLAG,
C        = 0 FOR NORMAL RETURN,
C        = 1 IF N IS LESS THAN 2,
C        = 2 IF X-VALUES ARE NOT STRICTLY INCREASING.
C
C AND
C
C   N, X, AND SIGMA ARE UNALTERED.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES TERMS AND
C SNHCSH.
C
C-----------------------------------------------------------
C
      NM1 = N-1
      IERR = 0
      IF (N .LE. 1) GO TO 3
      IF (X(N) .LE. X(1)) GO TO 4
C
C DENORMALIZE TENSION FACTOR
C
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
C
C GENERATE COEFFICIENTS FOR LEFT END BASIS FUNCTIONS
C
      D3 = X(2)-X(1)
      IF (D3 .LE. 0.) GO TO 4
      CALL TERMS (DIAG3,SDIAG3,SIGMAP,D3)
      D4 = D3
      IF (N .GE. 3) D4 = X(3)-X(2)
      IF (D4 .LE. 0.) GO TO 4
      CALL TERMS (DIAG4,SDIAG4,SIGMAP,D4)
      A22 = DIAG3+SDIAG3
      A23 = DIAG3+DIAG4+SDIAG3+SDIAG4
      V(2,1) = 0.
      V(3,1) = 1./(D3*(DIAG3+DIAG4)+(D3+D4)*SDIAG4)
      V(5,1) = SDIAG4*D4*V(3,1)
      IF (N .EQ. 2) GO TO 2
      A22 = 2.*A22
      D1 = D3
      D2 = D3
      D3 = D4
      DIAG1 = DIAG3
      DIAG2 = DIAG3
      DIAG3 = DIAG4
      SDIAG1 = SDIAG3
      SDIAG2 = SDIAG3
      SDIAG3 = SDIAG4
C
C GENERATE COEFFICIENTS FOR INTERIOR BASIS FUNCTIONS
C
      DO 1 I = 2,NM1
        IF (I .NE. NM1) D4 = X(I+2)-X(I+1)
        IF (D4 .LE. 0.) GO TO 4
        IF (D4 .NE. D3) CALL TERMS (DIAG4,SDIAG4,SIGMAP,D4)
        A11 = DIAG1+DIAG2+SDIAG1*(1.+D1/D2)
        A12 = SDIAG2/A11
        B1 = 1./(D2*A11)
        A33 = DIAG3+DIAG4+SDIAG4*(1.+D4/D3)
        A32 = SDIAG3/A33
        B3 = 1./(D3*A33)
        A21 = A22
        A22 = A23
        A23 = DIAG3+DIAG4+SDIAG3+SDIAG4
        V(2,I) = -(A21*B1+A23*B3)/(A22-A21*A12-A23*A32)
        V(1,I) = B1-A12*V(2,I)
        V(3,I) = B3-A32*V(2,I)
        V(4,I) = SDIAG1*D1*V(1,I)
        V(5,I) = SDIAG4*D4*V(3,I)
C
C SAVE CONSTANTS FOR NEXT ITERATION
C
        D1 = D2
        D2 = D3
        D3 = D4
        DIAG1 = DIAG2
        DIAG2 = DIAG3
        DIAG3 = DIAG4
        SDIAG1 = SDIAG2
        SDIAG2 = SDIAG3
    1   SDIAG3 = SDIAG4
C
C GENERATE COEFFICIENTS FOR RIGHT END BASIS FUNCTIONS
C
      V(2,N) = 0.
      V(1,N) = 1./(D2*(DIAG1+DIAG2)+(D2+D1)*SDIAG1)
      V(4,N) = SDIAG1*D1*V(1,N)
      V(3,N) = V(1,3)
      V(5,N) = V(4,3)
      V(1,1) = V(3,N-2)
      V(4,1) = V(5,N-2)
C
C ADJUST BASIS FOR NATURAL END CONDITIONS
C
      V(4,2) = V(4,2)-V(1,2)*V(5,N)/V(3,N)
      V(1,2) = 0.
      V(5,NM1) = V(5,NM1)-V(3,NM1)*V(4,1)/V(1,1)
      V(3,NM1) = 0.
      RETURN
C
C N EQUAL TO 2
C
    2 V(4,1) = V(5,1)
      V(1,1) = V(3,1)
      V(3,1) = 0.
      V(5,1) = 0.
      V(1,2) = 0.
      V(2,2) = 0.
      V(3,2) = V(1,1)
      V(4,2) = 0.
      V(5,2) = V(4,1)
      RETURN
C
C TOO FEW KNOTS
C
    3 IERR = 1
      RETURN
C
C X-VALUES NOT STRICTLY INCREASING
C
    4 IERR = 2
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TERMS (DIAG,SDIAG,SIGMA,DEL)
C
      REAL DIAG,SDIAG,SIGMA,DEL
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                       CODED BY A. K. CLINE AND R. J. RENKA
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE COMPUTES THE DIAGONAL AND SUPERDIAGONAL
C TERMS OF THE TRIDIAGONAL LINEAR SYSTEM ASSOCIATED WITH
C SPLINE UNDER TENSION INTERPOLATION.
C
C ON INPUT--
C
C   SIGMA CONTAINS THE TENSION FACTOR.
C
C AND
C
C   DEL CONTAINS THE STEP SIZE.
C
C ON OUTPUT--
C
C               (SIGMA*DEL*COSH(SIGMA*DEL) - SINH(SIGMA*DEL)
C   DIAG = DEL*--------------------------------------------.
C                     (SIGMA*DEL)**2 * SINH(SIGMA*DEL)
C
C                   SINH(SIGMA*DEL) - SIGMA*DEL
C   SDIAG = DEL*----------------------------------.
C                (SIGMA*DEL)**2 * SINH(SIGMA*DEL)
C
C AND
C
C   SIGMA AND DEL ARE UNALTERED.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULE SNHCSH.
C
C-----------------------------------------------------------
C
      IF (SIGMA .NE. 0.) GO TO 1
      DIAG = DEL/3.
      SDIAG = DEL/6.
      RETURN
    1 SIGDEL = SIGMA*DEL
      CALL SNHCSH (SINHM,COSHM,SIGDEL,0)
      DENOM = DEL/((SINHM+SIGDEL)*SIGDEL*SIGDEL)
      DIAG = DENOM*(SIGDEL*COSHM-SINHM)
      SDIAG = DENOM*SINHM
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TRIDEC (N,SUBDI,DIAGI,SUPD,SUBDO,DIAGO,
     *                   IERR)
C
      INTEGER N,IERR
      REAL SUBDI(N),DIAGI(N),SUPD(N),SUBDO(N),DIAGO(N)
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                                 CODED BY ALAN KAYLOR CLINE
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE FACTORIZES A TRIDIAGONAL MATRIX IN ORDER
C TO SOLVE SYSTEMS OF LINEAR EQUATIONS. THE FACTORIZATION
C EMPLOYS GAUSSIAN ELIMINATION WITHOUT ANY INTERCHANGING OF
C COLUMNS OR ROWS. THE SUBROUTINE TRISOL MAY BE CALLED TO
C ACTUALLY SOLVE THE SYSTEM ONCE THE FACTORIZATION HAS BEEN
C PERFORMED.
C
C ON INPUT--
C
C   N CONTAINS THE ORDER OF THE MATRIX (N .GE. 1).
C
C   SUBDI IS AN ARRAY CONTAINING THE SUBDIAGONAL ELEMENTS OF
C   THE MATRIX IN POSITIONS 2, ... , N.
C
C   DIAGI IS AN ARRAY CONTAINING THE DIAGONAL ELEMENTS OF
C   THE MATRIX.
C
C   SUPD IS AN ARRAY CONTAINING THE SUPERDIAGONAL ELEMENTS
C   OF THE MATRIX IN POSITIONS 1, ... , N-1.
C
C AND
C
C   SUBDO AND DIAGO ARE ARRAYS OF LENGTH N. (THE STORAGE
C   FOR THESE MAY COINCIDE WITH THAT FOR SUBDI AND DIAGI,
C   RESPECTIVELY, IN WHICH CASE THE ORIGINAL CONTENTS OF
C   SUBDI AND DIAGI WILL BE DESTROYED.)
C
C ON OUTPUT--
C
C   SUBDO AND DIAGO CONTAIN THE SUBDIAGONAL AND DIAGONAL OF
C   THE FACTORIZATION MATRIX.
C
C   IERR CONTAINS AN ERROR FLAG,
C        = 0 FOR NORMAL RETURN,
C        = 1 IF N IS LESS THAN 1,
C        = 2 IF THE SYSTEM IS SINGULAR.
C
C AND
C
C   N, SUBDI, DIAGI, AND SUPD ARE UNALTERED (UNLESS STORAGE
C   FOR SUBDI OR DIAGI COINCIDED WITH THAT FOR SUBDO
C   OR DIAGO, RESPECTIVELY).
C
C-----------------------------------------------------------
C
      IF (N .LE. 0) GO TO 3
      IERR = 2
      DIAGO(1) = DIAGI(1)
      IF (N .EQ. 1) GO TO 2
C
C FORWARD ELIMINATION
C
      DO 1 I = 2,N
        IM1 = I-1
        IF (DIAGO(IM1) .EQ. 0.) RETURN
        DIAGO(IM1) = 1./DIAGO(IM1)
        SUBDO(I) = SUBDI(I)*DIAGO(IM1)
    1   DIAGO(I) = DIAGI(I)-SUBDO(I)*SUPD(IM1)
    2 IF (DIAGO(N) .EQ. 0.) RETURN
      DIAGO(N) = 1./DIAGO(N)
      IERR = 0
      RETURN
C
C N LESS THAN 1
C
    3 IERR = 1
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TRISOL (N,SUBD,DIAG,SUPD,RHS,MRHS,NUMRHS,
     *                   INCRHS)
C
      INTEGER N,MRHS,NUMRHS,INCRHS
      REAL SUBD(N),DIAG(N),SUPD(N)
      REAL RHS(1+INCRHS*(N-1)+MRHS*(NUMRHS-1))
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                                 CODED BY ALAN KAYLOR CLINE
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C                               REVISED -- DECEMBER 31, 1992
C                                            BY LUDEK KLIMES
C                                   INSTITUTE OF GEOTECHNICS
C                               CZECHOSL. ACAD. SCI., PRAGUE
C
C THIS SUBROUTINE SOLVES TRIDIAGONAL SYSTEMS OF LINEAR
C EQUATIONS WITH MULTIPLE RIGHT HAND SIDES. THE RIGHT HAND
C SIDES MAY BE STORED ROW-WISE OR COLUMN-WISE. THE
C SUBROUTINE TRIDEC SHOULD BE CALLED EARLIER TO DETERMINE A
C FACTORIZATION OF THE TRIDIAGONAL MATRIX. THE SOLUTION
C VECTORS OVER-WRITE THE RIGHT HAND SIDES.
C
C ON INPUT--
C
C   N CONTAINS THE ORDER OF THE MATRIX (N .GE. 1).
C
C   SUBD, DIAG, AND SUPD ARE ARRAYS OF LENGTH N CONTAINING
C   THE SUBDIAGONAL, DIAGONAL, AND SUPERDIAGONAL OF THE
C   FACTORIZATION, RESPECTIVELY.
C
C   RHS IS AN ARRAY CONTAINING THE RIGHT HAND SIDES OF THE
C   TRIDIAGONAL SYSTEM.
C
C   MRHS IS THE INCREMENT BETWEEN THE FIRST COMPONENTS OF
C   EACH OF THE RIGHT HAND SIDE VECTORS IN STORAGE (MRHS
C   .GE. 1).
C
C   NUMRHS IS THE NUMBER OF RIGHT HAND SIDES TO BE SOLVED.
C
C AND
C
C   INCRHS IS THE INCREMENT BETWEEN COMPONENTS WITHIN EACH
C   OF THE RIGHT HAND SIDE VECTORS IN STORAGE (INCRHS .GE.
C   1).
C
C THE PARAMETERS N, SUBD, DIAG, AND SUPD MAY BE INPUT AS THE
C PARAMETERS N, SUBDO, DIAGO, AND SUPD OUTPUT BY SUBROUTINE
C TRIDEC, RESPECTIVELY.
C
C ON OUTPUT--
C
C   RHS CONTAINS THE SOLUTION VECTORS IN THE SAME STORAGE
C   STRUCTURE AS FOR THE RIGHT HAND SIDES.
C
C AND
C
C   N, SUBD, DIAG, SUPD, MRHS, NUMRHS, AND INCRHS ARE
C   UNALTERED.
C
C-----------------------------------------------------------
C
      NP1 = N+1
C
C LOOP ON RIGHT HAND SIDES
C
      DO 4 K = 1,NUMRHS
C
C FORWARD ELIMINATION
C
        IRHS = 1+MRHS*(K-1)
        IF (N .EQ. 1) GO TO 2
        DO 1 I = 2,N
          IM1RHS = IRHS
          IRHS = IRHS+INCRHS
    1     RHS(IRHS) = RHS(IRHS)-SUBD(I)*RHS(IM1RHS)
C
C BACK SUBSTITUTION
C
    2   RHS(IRHS) = DIAG(N)*RHS(IRHS)
        IF (N .EQ. 1) GO TO 4
        DO 3 IBAK = 2,N
          I = NP1-IBAK
          RHS(IM1RHS) = DIAG(I)*(RHS(IM1RHS)-SUPD(I)
     *                                           *RHS(IRHS))
          IRHS = IM1RHS
    3     IM1RHS = IM1RHS-INCRHS
    4   CONTINUE
      RETURN
      END
C
C=======================================================================
C
C PART 2:
C
C=======================================================================
C
      SUBROUTINE CURV2D (T,YY,YX,YXX,N,X,Y,YP,SIGMA)
C
      INTEGER N
      REAL T,YY,YX,YXX,X(N),Y(N),YP(N),SIGMA
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                                 CODED BY ALAN KAYLOR CLINE
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE DETERMINES FUNCTION VALUE, FIRST, AND
C SECOND DERIVATIVES OF A CURVE AT A GIVEN POINT USING A
C SPLINE UNDER TENSION. THE SUBROUTINE CURV1 SHOULD BE
C CALLED EARLIER TO DETERMINE CERTAIN NECESSARY PARAMETERS.
C
C ON INPUT--
C
C   T CONTAINS A REAL VALUE AT WHICH THE FUNCTION AND
C   DERIVATIVES ARE TO BE EVALUATED.
C
C   N CONTAINS THE NUMBER OF POINTS WHICH WERE SPECIFIED TO
C   DETERMINE THE CURVE.
C
C   X AND Y ARE ARRAYS CONTAINING THE ABSCISSAE AND
C   ORDINATES, RESPECTIVELY, OF THE SPECIFIED POINTS.
C
C   YP IS AN ARRAY OF SECOND DERIVATIVE VALUES OF THE CURVE
C   AT THE NODES.
C
C AND
C
C   SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).
C
C THE PARAMETERS N, X, Y, YP, AND SIGMA SHOULD BE INPUT
C UNALTERED FROM THE OUTPUT OF CURV1.
C
C ON OUTPUT--
C
C   YY, YX, AND YXX CONTAIN THE FUNCTION VALUE, FIRST AND
C   SECOND DERIVATIVES, RESPECTIVELY.
C
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES INTRVL AND
C SNHCSH.
C
C-----------------------------------------------------------
C
C DETERMINE INTERVAL
C
      IM1 = INTRVL(T,X,N)
      I = IM1+1
C
C DENORMALIZE TENSION FACTOR
C
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
C
C SET UP AND PERFORM INTERPOLATION
C
      DEL1 = T-X(IM1)
      DEL2 = X(I)-T
      DELS = X(I)-X(IM1)
      YY = (Y(I)*DEL1+Y(IM1)*DEL2)/DELS
      YX = (Y(I)-Y(IM1))/DELS
      IF (SIGMAP .NE. 0.) GO TO 1
      YY = YY-DEL1*DEL2*(YP(I)*(DEL1+DELS)+YP(IM1)*
     *        (DEL2+DELS))/(6.*DELS)
      YX = YX+(YP(I)*(2.*DEL1*DEL1-DEL2*(DEL1+DELS))-
     *         YP(IM1)*(2.*DEL2*DEL2-DEL1*(DEL2+DELS)))
     *         /(6.*DELS)
      YXX = (YP(I)*DEL1+YP(IM1)*DEL2)/DELS
      RETURN
    1 DELP1 = SIGMAP*(DEL1+DELS)/2.
      DELP2 = SIGMAP*(DEL2+DELS)/2.
      CALL SNHCSH (SINHM1,COSHM1,SIGMAP*DEL1,0)
      CALL SNHCSH (SINHM2,COSHM2,SIGMAP*DEL2,0)
      CALL SNHCSH (SINHMS,DUMMY,SIGMAP*DELS,-1)
      CALL SNHCSH (SINHP1,DUMMY,SIGMAP*DEL1/2.,-1)
      CALL SNHCSH (SINHP2,DUMMY,SIGMAP*DEL2/2.,-1)
      CALL SNHCSH (DUMMY,COSHP1,DELP1,1)
      CALL SNHCSH (DUMMY,COSHP2,DELP2,1)
      YY = YY+(YP(I)*(SINHM1*DEL2-DEL1*(2.*(COSHP1+1.)*
     *        SINHP2+SIGMAP*COSHP1*DEL2))+YP(IM1)*(SINHM2*
     *        DEL1-DEL2*(2.*(COSHP2+1.)*SINHP1+SIGMAP*
     *        COSHP2*DEL1)))/(SIGMAP*SIGMAP*DELS*(SINHMS+
     *        SIGMAP*DELS))
      YX = YX+(YP(I)*(DELS*SIGMAP*COSHM1-SINHMS)-
     *     YP(IM1)*(DELS*SIGMAP*COSHM2-SINHMS))/
     *     (SIGMAP*SIGMAP*DELS*(SINHMS+SIGMAP*DELS))
      YXX = (YP(I)*(SINHM1+SIGMAP*DEL1)+YP(IM1)*(SINHM2+
     *       SIGMAP*DEL2))/(SINHMS+SIGMAP*DELS)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE CURVBD (XX,W,WX,WXX,NX,X,C,VX,SIGMA)
C
      INTEGER NX
      REAL XX,W,WX,WXX,X(NX),VX(5,NX),C(NX),SIGMA
C
C                                      COMPLEMENT TO FITPACK
C                                       BY ALAN KAYLOR CLINE
C                                   CODED -- OCTOBER 9, 1986
C                                            BY LUDEK KLIMES
C                                      INST. GEOL. GEOTECHN.
C                               CZECHOSL. ACAD. SCI., PRAGUE
C
C THIS SUBROUTINE EVALUATES THE FUNCTION VALUE, THE
C FIRST PARTIAL DERIVATIVE, AND THE SECOND PARTIAL
C DERIVATIVE OF A SPLINE UNDER TENSION IN ONE VARIABLE.
C
C ON INPUT--
C
C   XX CONTAINS THE X-COORDINATE OF THE POINT
C   AT WHICH THE INTERPOLATION IS TO BE PERFORMED
C
C   NX IS THE NUMBER OF GRID POINTS
C
C   X IS ARRAY CONTAINING THE X-GRID VALUES.
C
C   C IS AN ARRAY OF COEFFICIENTS DESCRIBING THE FUNCTION IN
C   TERMS OF A B-SPLINE UNDER TENSION BASIS. IN THE
C   EXPANSION OF THE FUNCTION, FOR I = 1,...,NX ,
C   THE COEFFICIENT MULTIPLYING THE BASIS
C   FUNCTION I IS STORED IN C(I).
C
C   VX IS THE ARRAY OF LENGTH 5*NX
C   CONTAINING THE B-SPLINE BASIS DATA
C
C   SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).
C
C THE PARAMETERS NX, X, C, VX, AND SIGMA
C SHOULD BE INPUT UNALTERED FROM THE OUTPUT OF CURVB1.
C
C ON OUTPUT--
C
C   W CONTAINS THE INTERPOLATED FUNCTION VALUE.
C
C   WX CONTAINS THE FIRST DERIVATIVE .
C
C   WXX CONTAINS THE SECOND DERIVATIVE .
C
C AND
C
C   NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES DSPLNZ, INTRVL,
C AND SNHCSH.
C
C--------------------------------------------------------------
C
      REAL BX(3,4)
C
C EVALUATE BASIS FUNCTIONS AT XX
C
      CALL DSPLNZ (XX,NX,X,VX,SIGMA,ISTART,BX)
C
C ACCUMULATE BASIS FUNCTIONS
C
      SUM = 0.
      SUMX = 0.
      SUMXX = 0.
      DO 1 I = 1,4
        II = ISTART+I-1
        IF (II .EQ. 0 .OR. II .GT. NX) GO TO 1
        BX1I = BX(1,I)
        CI  = C(II)
        SUM = SUM+CI*BX1I
        SUMX = SUMX+CI*BX(2,I)
        SUMXX = SUMXX+CI*BX(3,I)
V       CALL VAR2(II,BX1I,BX(2,I),0.,0.)
    1   CONTINUE
      W = SUM
      WX = SUMX
      WXX = SUMXX
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE SURFBD (XX,YY,W,WX,WY,WXX,WXY,WYY,NX,NY,X,
     *                   Y,C,VX,VY,SIGMA)
C
      INTEGER NX,NY
      REAL XX,YY,W,WX,WY,WXX,WXY,WYY,X(NX),Y(NY),VX(5,NX),
     *     VY(5,NY),C(NX,NY),SIGMA
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                                 CODED BY ALAN KAYLOR CLINE
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE EVALUATES THE FUNCTION VALUE, THE TWO
C FIRST PARTIAL DERIVATIVES, AND THE SIX SECOND PARTIAL
C DERIVATIVES OF A TENSOR PRODUCT SPLINE UNDER TENSION IN
C TWO VARIABLES.
C
C ON INPUT--
C
C   XX AND YY CONTAIN THE X- AND Y-COORDINES OF THE POINT
C   AT WHICH THE INTERPOLATION IS TO BE PERFORME .D.
C
C   NX AND NY ARE THE NUMBER OF GRID LINES IN THE X- AND Y-
C   DIRECTIONS, RESPECTIVELY, OF THE RECTANGULAR GRID WHICH
C   SPECIFIED THE FUNCTION.
C
C   X AND Y ARE ARRAYS CONTAINING THE X- AND Y-GRID VALUES,
C   RESPECTIVELY.
C
C   C IS AN ARRAY OF COEFFICIENTS DESCRIBING THE FUNCTION IN
C   TERMS OF A B-SPLINE UNDER TENSION BASIS. IN THE
C   EXPANSION OF THE FUNCTION, FOR I = 1,...,NX AND J = 1,
C   ...,NY, THE COEFFICIENT MULTIPLYING THE PRODUCT OF BASIS
C   FUNCTION I IN X AND BASIS FUNCTION J IN Y IS STORED IN
C   C(I,J).
C
C   VX AND VY VZ ARE ARRAYS OF LENGTH 5*NX AND 5*NY,
C   RESPECTIVELY, CONTAINING THE B-SPLINE BASIS DATA FOR THE
C   X- AND Y-GRIDS.
C
C AND
C
C   SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).
C
C THE PARAMETERS NX, NY, X, Y, Z, C, VX, VY, AND SIGMA
C SHOULD BE INPUT UNALTERED FROM THE OUTPUT OF SURFB1.
C
C ON OUTPUT--
C
C   W CONTAINS THE INTERPOLATED FUNCTION VALUE.
C
C   WX AND WY CONTAIN THE X- AND Y-PARTIAL DERIVATIVES,
C   RESPECTIVELY.
C
C   WXX, WXY, AND WYY CONTAIN THE XX-, XY-, AND YY-PARTIAL
C   DERIVATIVES, RESPECTIVELY.
C
C AND
C
C   NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES DSPLNZ, INTRVL,
C AND SNHCSH.
C
C--------------------------------------------------------- ----
C
      REAL BX(3,4),BY(3,4)
C
C EVALUATE BASIS FUNCTIONS AT XX AND YY
C
      CALL DSPLNZ (XX,NX,X,VX,SIGMA,ISTART,BX)
      CALL DSPLNZ (YY,NY,Y,VY,SIGMA,JSTART,BY)
C
C ACCUMULATE TENSOR PRODUCTS
C
      SUM = 0.
      SUMX = 0.
      SUMY = 0.
      SUMXX = 0.
      SUMXY = 0.
      SUMYY = 0.
      DO 2 J = 1,4
        JJ = JSTART+J-1
        IF (JJ .EQ. 0 .OR. JJ .GT. NY) GO TO 2
        BY1J = BY(1,J)
        BY2J = BY(2,J)
        BY3J = BY(3,J)
        DO 1 I = 1,4
          II = ISTART+I-1
          IF (II .EQ. 0 .OR. II .GT. NX) GO TO 1
          BX1I = BX(1,I)
          BX2I =  BX(2,I)
          CIJ = C(II,JJ)
          SUM = SUM+CIJ*BX1I*BY1J
          SUMX = SUMX+CIJ*BX2I*BY1J
          SUMY = SUMY+CIJ*BX1I*BY2J
          SUMXX = SUMXX+CIJ*BX(3,I)*BY1J
          SUMXY = SUMXY+CIJ*BX2I*BY2J
          SUMYY = SUMYY+CIJ*BX1I*BY3J
V         CALL VAR2(II+NX*(JJ-1),BX1I*BY1J,BX2I*BY1J,BX1I*BY2J,0.)
    1     CONTINUE
    2   CONTINUE
      W = SUM
      WX = SUMX
      WY = SUMY
      WXX = SUMXX
      WXY = SUMXY
      WYY = SUMYY
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE VAL3BD (XX,YY,ZZ,W,WX,WY,WZ,WXX,WXY,WYY,
     *                   WYZ,WZZ,WXZ,NX,NY,NZ,X,Y,Z,C,VX,VY,
     *                   VZ,SIGMA)
C
      INTEGER NX,NY,NZ
      REAL XX,YY,ZZ,W,WX,WY,WZ,WXX,WXY,WYY,WYZ,WZZ,WXZ,
     *     X(NX),Y(NY),Z(NZ),VX(5,NX),VY(5,NY),VZ(5,NZ),
     *     C(NX,NY,NZ),SIGMA
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                                 CODED BY ALAN KAYLOR CLINE
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE EVALUATES THE FUNCTION VALUE, THE THREE
C FIRST PARTIAL DERIVATIVES, AND THE SIX SECOND PARTIAL
C DERIVATIVES OF A TENSOR PRODUCT SPLINE UNDER TENSION IN
C THREE VARIABLES.
C
C ON INPUT--
C
C   XX, YY, AND ZZ CONTAIN THE X-, Y-, AND Z-COORDINATES OF
C   THE POINT AT WHICH THE INTERPOLATION IS TO BE PERFORMED.
C
C   NX, NY, AND NZ ARE THE NUMBER OF GRID LINES IN THE X-,
C   Y-, AND Z-DIRECTIONS, RESPECTIVELY, OF THE RECTANGULAR
C   GRID WHICH SPECIFIED THE FUNCTION.
C
C   X, Y, AND Z ARE ARRAYS CONTAINING THE X-, Y-, AND Z-GRID
C   VALUES, RESPECTIVELY.
C
C   C IS AN ARRAY OF COEFFICIENTS DESCRIBING THE FUNCTION IN
C   TERMS OF A B-SPLINE UNDER TENSION BASIS. IN THE
C   EXPANSION OF THE FUNCTION, FOR I = 1,...,NX, J = 1,...,
C   NY, AND K = 1,...,NZ, THE COEFFICIENT MULTIPLYING THE
C   PRODUCT OF BASIS FUNCTION I IN X, BASIS FUNCTION J IN Y,
C   AND BASIS FUNCTION K IN Z IS STORED IN C(I,J,K).
C
C   VX, VY, AND VZ ARE ARRAYS OF LENGTH 5*NX, 5*NY, AND
C   5*NZ, RESPECTIVELY, CONTAINING THE B-SPLINE BASIS DATA
C   FOR THE X-, Y-, AND Z-GRIDS.
C
C AND
C
C   SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).
C
C THE PARAMETERS NX, NY, NZ, X, Y, Z, C, VX, VY, VZ, AND
C SIGMA SHOULD BE INPUT UNALTERED FROM THE OUTPUT OF
C VAL3B1.
C
C ON OUTPUT--
C
C   W CONTAINS THE INTERPOLATED FUNCTION VALUE.
C
C   WX, WY, AND WZ CONTAIN THE X-, Y-, AND Z-PARTIAL
C   DERIVATIVES, RESPECTIVELY.
C
C   WXX, WXY, WYY, WYZ, WZZ, AND WXZ CONTAIN THE XX-, XY-
C   YY-, YZ-, ZZ-, AND XZ-PARTIAL DERIVATIVES, RESPECTIVELY.
C
C AND
C
C   NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES DSPLNZ, INTRVL,
C AND SNHCSH.
C
C--------------------------------------------------------------
C
      REAL BX(3,4),BY(3,4),BZ(3,4)
C
C EVALUATE BASIS FUNCTIONS AT XX, YY, AND ZZ
C
      CALL DSPLNZ (XX,NX,X,VX,SIGMA,ISTART,BX)
      CALL DSPLNZ (YY,NY,Y,VY,SIGMA,JSTART,BY)
      CALL DSPLNZ (ZZ,NZ,Z,VZ,SIGMA,KSTART,BZ)
C
C ACCUMULATE TENSOR PRODUCTS
C
      SUM = 0.
      SUMX = 0.
      SUMY = 0.
      SUMZ = 0.
      SUMXX = 0.
      SUMXY = 0.
      SUMYY = 0.
      SUMYZ = 0.
      SUMZZ = 0.
      SUMXZ = 0.
      DO 3 K = 1,4
        KK = KSTART+K-1
        IF (KK .EQ. 0 .OR. KK .GT. NZ) GO TO 3
        BZ1K = BZ(1,K)
        BZ2K = BZ(2,K)
        BZ3K = BZ(3,K)
        DO 2 J = 1,4
          JJ = JSTART+J-1
          IF (JJ .EQ. 0 .OR. JJ .GT. NY) GO TO 2
          BY1J = BY(1,J)
          BY2J = BY(2,J)
          BY3J = BY(3,J)
          DO 1 I = 1,4
            II = ISTART+I-1
            IF (II .EQ. 0 .OR. II .GT. NX) GO TO 1
            BX1I = BX(1,I)
            BX2I =  BX(2,I)
            CIJK = C(II,JJ,KK)
            SUM = SUM+CIJK*BX1I*BY1J*BZ1K
            SUMX = SUMX+CIJK*BX2I*BY1J*BZ1K
            SUMY = SUMY+CIJK*BX1I*BY2J*BZ1K
            SUMZ = SUMZ+CIJK*BX1I*BY1J*BZ2K
            SUMXX = SUMXX+CIJK*BX(3,I)*BY1J*BZ1K
            SUMXY = SUMXY+CIJK*BX2I*BY2J*BZ1K
            SUMYY = SUMYY+CIJK*BX1I*BY3J*BZ1K
            SUMYZ = SUMYZ+CIJK*BX1I*BY2J*BZ2K
            SUMZZ = SUMZZ+CIJK*BX1I*BY1J*BZ3K
            SUMXZ = SUMXZ+CIJK*BX2I*BY1J*BZ2K
V           CALL VAR2(II+NX*(JJ-1+NY*(KK-1)),BX1I*BY1J*BZ1K,
V    *                     BX2I*BY1J*BZ1K,BX1I*BY2J*BZ1K,BX1I*BY1J*BZ2K)
    1       CONTINUE
    2     CONTINUE
    3   CONTINUE
      W = SUM
      WX = SUMX
      WY = SUMY
      WZ = SUMZ
      WXX = SUMXX
      WXY = SUMXY
      WYY = SUMYY
      WYZ = SUMYZ
      WZZ = SUMZZ
      WXZ = SUMXZ
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE DSPLNZ (T,N,X,V,SIGMA,ISTART,B)
C
      INTEGER N,ISTART
      REAL T,X(N),V(5,N),SIGMA,B(3,4)
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                                 CODED BY ALAN KAYLOR CLINE
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE EVALUATES AT A GIVEN POINT THE FOUR NON-
C ZERO BASIS FUNCTIONS OF A B-SPLINE UNDER TENSION BASIS AND
C THEIR FIRST AND SECOND DERIVATIVES. THE INDEX OF THE FIRST
C NON-ZERO BASIS FUNCTION IS ALSO DETERMINED. (THE SENSE OF
C THE WORD NON-ZERO IS EXTENDED TO INCLUDE THE SPECIAL CASE
C WHERE THE GIVEN POINT COINCIDES WITH A KNOT IN WHICH CASE
C THE LAST OF THE FOUR RETURNED FUNCTION VALUES MAY BE ZERO.
C ) THE SUBROUTINE VGEN SHOULD BE CALLED EARLIER TO
C DETERMINE CERTAIN NECESSARY COEFFICIENTS.
C
C ON INPUT--
C
C   T CONTAINS A REAL VALUE AT WHICH THE BASIS FUNCTIONS ARE
C   TO BE EVALUATED.
C
C   N CONTAINS THE NUMBER OF KNOTS DEFINING THE BASIS.
C
C   X CONTAINS THE ARRAY OF KNOTS.
C
C   V CONTAINS THE ARRAY OF COEFFICIENTS DETERMINED BY VGEN
C   FOR CALCULATION OF BASIS FUNCTIONS.
C
C   SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).
C
C   ISTART IS AN INTEGER VARIABLE.
C
C AND
C
C   B IS A REAL ARRAY WITH 3 ROWS AND 4 COLUMNS.
C
C THE PARAMETERS N, X, V, AND SIGMA SHOULD BE INPUT
C UNALTERED FROM THE OUTPUT OF VGEN.
C
C ON OUTPUT--
C
C   ISTART CONTAINS THE INDEX OF THE FIRST NON-ZERO BASIS
C   FUNCTION. THUS 0 .LE. ISTART .LE. N-2 AND THE N0N-ZERO
C   BASIS FUNCTIONS HAVE INDICES ISTART, ... , ISTART+3.
C
C   B CONTAINS THE VALUES AT T OF BASIS FUNCTIONS ISTART,
C   ... , ISTART+3 IN B(1,1), ... , B(1,4), RESPECTIVELY.
C   FIRST AND SECOND DERIVATIVES OF THE CORRESPONDING
C   FUNCTIONS ARE CONTAINED IN B(2,1), ... , B(2,4), AND
C   B(3,1), ... , B(3,4), RESPECTIVELY.
C
C   T, N, X, V, AND SIGMA ARE UNALTERED.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES INTRVL AND
C SNHCSH.
C
C-----------------------------------------------------------
C
C DENORMALIZE TENSION FACTOR
C
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
C
C DETERMINE INDEX OF FIRST NON-ZERO BASIS FUNCTION
C
      I = INTRVL (T,X,N)-1
C
C COMPUTE DISTANCES TO ADJACENT KNOTS AND LAGRANGIAN
C WEIGHTS
C
      DEL1 = T-X(I+1)
      DEL2 = X(I+2)-T
      DELS = X(I+2)-X(I+1)
      C10 = DEL2/DELS
      C20 = DEL1/DELS
      C11 = -1./DELS
      C21 = 1./DELS
      IF (SIGMAP .NE. 0.) GO TO 1
      FAC = -DEL1*DEL2/(6.*DELS)
      CP10 = FAC*(DEL2+DELS)
      CP20 = FAC*(DEL1+DELS)
      CP11 = -(2.*DEL2*DEL2-DEL1*(DEL2+DELS))/(6.*DELS)
      CP21 = (2.*DEL1*DEL1-DEL2*(DEL1+DELS))/(6.*DELS)
      CP12 = C10
      CP22 = C20
      GO TO 2
    1 DELP1 = SIGMAP*(DEL1+DELS)/2.
      DELP2 = SIGMAP*(DEL2+DELS)/2.
      CALL SNHCSH (SINHM1,COSHM1,SIGMAP*DEL1,0)
      CALL SNHCSH (SINHM2,COSHM2,SIGMAP*DEL2,0)
      CALL SNHCSH (SINHMS,DUMMY,SIGMAP*DELS,-1)
      CALL SNHCSH (SINHP1,DUMMY,SIGMAP*DEL1/2.,-1)
      CALL SNHCSH (SINHP2,DUMMY,SIGMAP*DEL2/2.,-1)
      CALL SNHCSH (DUMMY,COSHP1,DELP1,1)
      CALL SNHCSH (DUMMY,COSHP2,DELP2,1)
      SINHS = SINHMS+SIGMAP*DELS
      DENOM = SIGMAP*SIGMAP*DELS*SINHS
      CP10 = (SINHM2*DEL1-DEL2*(2.*(COSHP2+1.)*SINHP1
     *       +SIGMAP*COSHP2*DEL1))/DENOM
      CP20 = (SINHM1*DEL2-DEL1*(2.*(COSHP1+1.)*SINHP2
     *       +SIGMAP*COSHP1*DEL2))/DENOM
      CP11 = -(DELS*SIGMAP*COSHM2-SINHMS)/DENOM
      CP21 = (DELS*SIGMAP*COSHM1-SINHMS)/DENOM
      CP12 = (SINHM2+SIGMAP*DEL2)/SINHS
      CP22 = (SINHM1+SIGMAP*DEL1)/SINHS
C
C COMPUTE BASIS FUNCTION VALUES
C
    2 II = I
      IF (II .EQ. 0) II = N
      IIP1 = I+1
      IIP2 = I+2
      IIP3 = I+3
      IF (IIP2 .EQ. N) IIP3 = 1
      B(1,1) = C10*V(5,II)+CP10*V(3,II)
      B(1,2) = C10+C20*V(5,IIP1)+CP10*V(2,IIP1)+
     *         CP20*V(3,IIP1)
      B(1,3) = C10*V(4,IIP2)+C20+CP10*V(1,IIP2)+
     *         CP20*V(2,IIP2)
      B(1,4) = C20*V(4,IIP3)+CP20*V(1,IIP3)
      B(2,1) = C11*V(5,II)+CP11*V(3,II)
      B(2,2) = C11+C21*V(5,IIP1)+CP11*V(2,IIP1)+
     *         CP21*V(3,IIP1)
      B(2,3) = C11*V(4,IIP2)+C21+CP11*V(1,IIP2)+
     *         CP21*V(2,IIP2)
      B(2,4) = C21*V(4,IIP3)+CP21*V(1,IIP3)
      B(3,1) = CP12*V(3,II)
      B(3,2) = CP12*V(2,IIP1)+CP22*V(3,IIP1)
      B(3,3) = CP12*V(1,IIP2)+CP22*V(2,IIP2)
      B(3,4) = CP22*V(1,IIP3)
      ISTART = I
      RETURN
      END
C
C=======================================================================
C
      FUNCTION INTRVL (T,X,N)
C
      INTEGER N
      REAL T,X(N)
C
C                            FROM FITPACK -- AUGUST 31, 1981
C                       CODED BY A. K. CLINE AND R. J. RENKA
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS FUNCTION DETERMINES THE INDEX OF THE INTERVAL
C (DETERMINED BY A GIVEN INCREASING SEQUENCE) IN WHICH
C A GIVEN VALUE LIES.
C
C ON INPUT--
C
C   T IS THE GIVEN VALUE.
C
C   X IS A VECTOR OF STRICTLY INCREASING VALUES.
C
C AND
C
C   N IS THE LENGTH OF X (N .GE. 2).
C
C ON OUTPUT--
C
C   INTRVL RETURNS AN INTEGER I SUCH THAT
C
C          I = 1       IF             T .LT. X(2)  ,
C          I = N-1     IF X(N-1) .LE. T            ,
C          OTHERWISE       X(I)  .LE. T .LT. X(I+1),
C
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C-----------------------------------------------------------
C
      TT = T
      IF (TT .LT. X(2)) GO TO 4
      IF (TT .GE. X(N-1)) GO TO 5
      IL = 2
      IH = N-1
C
C LINEAR INTERPOLATION
C
    1 I = MIN0(IL+IFIX(FLOAT(IH-IL)*(TT-X(IL))/(X(IH)-X(IL))),
     *         IH-1)
      IF (TT .LT. X(I)) GO TO 2
      IF (TT .LT. X(I+1)) GO TO 3
C
C TOO HIGH
C
      IL = I+1
      GO TO 1
C
C TOO LOW
C
    2 IH = I
      GO TO 1
    3 INTRVL = I
      RETURN
C
C LEFT END
C
    4 INTRVL = 1
      RETURN
C
C RIGHT END
C
    5 INTRVL = N-1
      RETURN
      END
C
C=======================================================================
C