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 2 columns. 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)
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
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
CALL VAR2(II+NX*(JJ-1+NY*(KK-1)),BX1I*BY1J*BZ1K,
* 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