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 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 CURVN1 C CURVB1 C SURFB1 C VAL3B1 C VGEN C TERMS C TRIDEC C TRISOL 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 CURV2D C CURVBD C SURFBD C VAL3BD C DSPLNZ C INTRVL 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 The lines denoted by '*V' in the first two columns of file C 'fit.for' calculate the model variations with respect to the model C parameters. C File 'fitv.for', intended for the model inversion, is created C from 'fit.for' by replacing each '*V' in the first two columns C by spaces using program 'clean.for'. Subroutines VAR4 and VAR5 C of file 'var.for' may then be called to handle the variations. C Note 4: 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 appearances 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 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 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 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 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 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 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. Function 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 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 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 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 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 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 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-coordinates of the point C at which the interpolation is to be performed. 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 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 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 non-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 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