C
SUBROUTINE SPSP(NAX,NAY,NAZ,NBX,NBY,NBZ,NX,NY,NZ, * X,Y,Z,VX,VY,VZ,SIGMA,WEIGHT,B,BX,BY,BZ) C INTEGER NAX,NAY,NAZ,NBX,NBY,NBZ,NX,NY,NZ REAL X(NX),Y(NY),Z(NZ),VX(5,NX),VY(5,NY),VZ(5,NZ) REAL SIGMA,WEIGHT,B(*),BX(*),BY(*),BZ(*) C C Complement to FITPACK C by Alan Kaylor Cline C coded -- January 23, 1994 C by Ludek Klimes C Department of Geophysics C Charles University, Prague C C This subroutine evaluates the Sobolev scalar products C of spline under tension basis functions in three variables C (the Sobolev scalar product consists of integrals of the C products of partial derivatives of the two argument functions) C C On input-- C C NXA, NYA, NZA are the orders of partial derivatives of C the first argument function in the scalar product C C NXB, NYB, NZB are the orders of partial derivatives of C the second argument function in the scalar product C C NX, NY, NZ are the numbers of grid points in the C X-, Y-, Z-directions, respectively. (NX, NY, NZ C should be at least 1) 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 VX, VY,VZ are arrays of lengths 5*NX, 5*NY, 5*NZ, C respectively, containing the B-spline basis data for the C X-, Y- and Z-grids. They contain certain coefficients C to be used 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 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 WEIGHT is the weight of the product of NXA,NYA,NZA-partial C derivative of the first argument and NXB,NYB,NZB-partial C derivative of the second argument, in the Sobolev scalar C product. The integral of the product of the partial C derivatives multiplied by WEIGHT is added to matrix B. C C B is the array containing NN*NN matrix B (NN=NX*NY*NZ), C stored as a symmetric matrix ( NN*(NN+1)/2 storage C locations ) if NAX.EQ.NBX and NAY.EQ.NBY and NAZ.EQ.NBZ, C else stored as a general matrix ( NN*NN storage C locations ). The II,JJ-element of the matrix B C will be increased by the integral of the product of C NXA-,NYA-,NZA-partial derivative of the II-th basis C function and NXB-,NYB-,NZB-partial derivative of the C JJ-th basis function, multiplied by WEIGHT. C Here the basis function IX,IY,IZ (1.LE.IX.LE.NX, C 1.LE.IY.LE.NY, 1.LE.IZ.LE.NZ) is indexed by C II=IX+NX*(IY+NY*IZ). C C BX is an auxiliary array of at least NX*(NX+1)/2 C locations for NXA.EQ.NXB, or of at least NX*NX locations C for NXA.NE.NXB. It is used for scratch storage. C C BY is an auxiliary array of at least NY*(NY+1)/2 C locations for NYA.EQ.NYB, or of at least NY*NY locations C for NYA.NE.NYB. It is used for scratch storage. C C BZ is an auxiliary array of at least NZ*(NZ+1)/2 C locations for NZA.EQ.NZB, or of at least NZ*NZ locations C for NZA.NE.NZB. It is used for scratch storage. C C And C C None of the input parameters, except B, BX, BY, BZ, are C altered C C The parameters NX, NY, NZ, X, Y, Z, VX, VY, VZ and SIGMA C should be input unaltered from the output of VAL3B1 C (SURFB1, CURVB1). C C On output-- C C B is the input array increased by the integrals of the C products of NXA-,NYA-,NZA-partial derivatives and C NXB-,NYB-,NZB-partial derivatives of the spline under C tension basis functions, multiplied by WEIGHT. C C This subroutine references package modules QSPL, QINT, C and SNHCSH. C EXTERNAL QSPL C C-------------------------------------------------------------- C C Other variables used inside the subroutine QSPL C INTEGER IX,JX,KX,MX,IY,JY,KY,MY,IZ,JZ,KZ,MZ,II,JJ,KK,MM C C The matrix element B(II,JJ) is located in the array element C B(KK), where C for symmetric matrix B, II.LE.JJ : C KK= (JJ-1)*JJ/2+II C for symmetric matrix B, II.GT.JJ : C KK= (II-1)*II/2+JJ C for nonsymmetric matrix B : C KK= (JJ-1)*NN+II C with NN=NX*NY*NZ being the dimension of the matrix B. C C The matrix element BX(IX,JX) is located in the array element C BX(KX). The meaning of IX,JX,KX is similar as the meaning C of II,JJ,KK in the case of matrix B. C C The matrix element BY(IY,JY) is located in the array element C BZ(KY). The meaning of IY,JY,KY is similar as the meaning C of II,JJ,KK in the case of matrix B. C C The matrix element BZ(IZ,JZ) is located in the array element C BZ(KZ). The meaning of IZ,JZ,KZ is similar as the meaning C of II,JJ,KK in the case of matrix B. C C MM, MX, MY, MZ are auxiliary variables considering the C symmetry of the matrices B, BX, BY, BZ. C C--------------------------------------------------------------- C C Scalar products of B-splines in X-direction KX= 0 MX= NX DO 12 JX=1,NX C Is BX symmetric matrix ? IF(NAX.EQ.NBX) MX=JX DO 11 IX=1,MX KX= KX+1 CALL QSPL(NAX,NBX,IX,JX,NX,X,VX,SIGMA,BX(KX)) C QSPL 11 CONTINUE 12 CONTINUE C C Scalar products of B-splines in Y-direction KY= 0 MY= NY DO 14 JY=1,NY C Is BY symmetric matrix ? IF(NAY.EQ.NBY) MY=JY DO 13 IY=1,MY KY= KY+1 CALL QSPL(NAY,NBY,IY,JY,NY,Y,VY,SIGMA,BY(KY)) C QSPL 13 CONTINUE 14 CONTINUE C C Scalar products of B-splines in Z-direction KZ= 0 MZ= NZ DO 16 JZ=1,NZ C Is BZ symmetric matrix ? IF(NAZ.EQ.NBZ) MZ=JZ DO 15 IZ=1,MZ KZ= KZ+1 CALL QSPL(NAZ,NBZ,IZ,JZ,NZ,Z,VZ,SIGMA,BZ(KZ)) C QSPL 15 CONTINUE 16 CONTINUE C C Scalar products of 3-D B-splines C Is B symmetric matrix ? IF(NAX.EQ.NBX.AND.NAY.EQ.NBY.AND.NAZ.EQ.NBZ) THEN MM= 1 ELSE MM= 0 END IF KK= 0 JJ= 0 DO 27 JZ=1,NZ DO 26 JY=1,NY DO 25 JX=1,NX JJ= JJ+1 II= 0 C Is BZ symmetric matrix ? IF(NAZ.EQ.NBZ) THEN KZ= (JZ-1)*JZ/2 ELSE KZ= (JZ-1)*NZ END IF DO 23 IZ=1,NZ KZ= KZ+1 C Subdiagonal element of matrix BZ IF(NAZ.EQ.NBZ.AND.IZ.GT.JZ) KZ=KZ+IZ-2 C Is BY symmetric matrix ? IF(NAY.EQ.NBY) THEN KY= (JY-1)*JY/2 ELSE KY= (JY-1)*NY END IF DO 22 IY=1,NY KY= KY+1 C Subdiagonal element of matrix BY IF(NAY.EQ.NBY.AND.IY.GT.JY) KY=KY+IY-2 C Is BX symmetric matrix ? IF(NAX.EQ.NBX) THEN KX= (JX-1)*JX/2 ELSE KX= (JX-1)*NX END IF DO 21 IX=1,NX KX= KX+1 C Subdiagonal element of matrix BX IF(NAX.EQ.NBX.AND.IX.GT.JX) KX=KX+IX-2 KK= KK+1 B(KK)= B(KK)+WEIGHT*BX(KX)*BY(KY)*BZ(KZ) II= II+1 IF(MM*II.GE.JJ) GO TO 24 21 CONTINUE 22 CONTINUE 23 CONTINUE 24 CONTINUE 25 CONTINUE 26 CONTINUE 27 CONTINUE C RETURN END C C======================================================================= C C C SUBROUTINE QSPL(NA,NB,IA,IB,N,X,V,SIGMA,Q) C INTEGER NA,NB,IA,IB,N REAL X(N),V(5,N),SIGMA,Q C C Complement to FITPACK C by Alan Kaylor Cline C coded -- January 23, 1994 C by Ludek Klimes C Department of Geophysics C Charles University, Prague C C This subroutine evaluates the Sobolev scalar product C of spline under tension basis functions in one variable C (the Sobolev scalar product consists of integrals of the C products of partial derivatives of the two argument functions) C C On input-- C C NA is the order of the partial derivative of C the first argument function in the scalar product. C C NB is the order of the partial derivative of C the second argument function in the scalar product. C C IA is the index of the first argument function C (1.LE.IA.LE.N). C C IB is the index of the second argument function C (1.LE.IB.LE.N). C C N is the number of grid points. C (N should be at least 1) C C X is the array of the N coordinates of grid points. C These should be strictly increasing. C C V is the array of lengths 5*N, C containing certain coefficients to be used C 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 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 None of the input parameters are altered. C C The parameters N, X, V, and SIGMA C should be input unaltered from the output of VAL3B1 C (SURFB1, CURVB1). C C On output-- C C Q is the integral of the product of NA-th partial C derivative of the IA-th basis function and C NB-th partial derivative of the IB-th spline under C tension basis function. C C This subroutine references package modules QINT, SNHCSH. C EXTERNAL QINT C C--------------------------------------------------------------- C C Other variables used inside the subroutine QSPL: C INTEGER I,J REAL SIGMAP,V1A,V2A,V3A,V4A,V5A,V1B,V2B,V3B,V4B,V5B C C I...Index of the interval. C J...Position of the second B-spline with respect to the C interval I. C SIGMAP...Denormalized tension factor. C V1A,V2A,V3A,V4A,V5A,V1B,V2B,V3B,V4B,V5B...Auxiliary C storage locations for V(1,IA),...,V(5,IB). C C--------------------------------------------------------------- C IF(N.GT.1) GO TO 10 Q = 1. IF(NA.NE.0.OR.NB.NE.0) Q=0. GO TO 90 C 10 SIGMAP= ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1)) V1A= V(1,IA) V2A= V(2,IA) V3A= V(3,IA) V4A= V(4,IA) V5A= V(5,IA) V1B= V(1,IB) V2B= V(2,IB) V3B= V(3,IB) V4B= V(4,IB) V5B= V(5,IB) Q = 0. C I = IA-2 IF(I.LT.1) GO TO 20 J = I-IB+3 IF(J.LT.1) GO TO 20 IF(J.GT.4) GO TO 90 GO TO (11,12,13,14),J 11 CALL QINT(X(I),X(I+1),0. ,0. ,V4A,V1A,NA, * 0. ,0. ,V4B,V1B,NB,SIGMAP,Q) GO TO 20 12 CALL QINT(X(I),X(I+1),0. ,0. ,V4A,V1A,NA, * V4B,V1B,1. ,V2B,NB,SIGMAP,Q) GO TO 20 13 CALL QINT(X(I),X(I+1),0. ,0. ,V4A,V1A,NA, * 1. ,V2B,V5B,V3B,NB,SIGMAP,Q) GO TO 20 14 CALL QINT(X(I),X(I+1),0. ,0. ,V4A,V1A,NA, * V5B,V3B,0. ,0. ,NB,SIGMAP,Q) C QINT C 20 I = IA-1 IF(I.LT.1) GO TO 30 J = I-IB+3 IF(J.LT.1) GO TO 30 IF(J.GT.4) GO TO 90 GO TO (21,22,23,24),J 21 CALL QINT(X(I),X(I+1),V4A,V1A,1. ,V2A,NA, * 0. ,0. ,V4B,V1B,NB,SIGMAP,Q) GO TO 30 22 CALL QINT(X(I),X(I+1),V4A,V1A,1. ,V2A,NA, * V4B,V1B,1. ,V2B,NB,SIGMAP,Q) GO TO 30 23 CALL QINT(X(I),X(I+1),V4A,V1A,1. ,V2A,NA, * 1. ,V2B,V5B,V3B,NB,SIGMAP,Q) GO TO 30 24 CALL QINT(X(I),X(I+1),V4A,V1A,1. ,V2A,NA, * V5B,V3B,0. ,0. ,NB,SIGMAP,Q) C QINT C 30 I = IA IF(I.GE.N) GO TO 90 J = I-IB+3 IF(J.LT.1) GO TO 40 IF(J.GT.4) GO TO 90 GO TO (31,32,33,34),J 31 CALL QINT(X(I),X(I+1),1. ,V2A,V5A,V3A,NA, * 0. ,0. ,V4B,V1B,NB,SIGMAP,Q) GO TO 40 32 CALL QINT(X(I),X(I+1),1. ,V2A,V5A,V3A,NA, * V4B,V1B,1. ,V2B,NB,SIGMAP,Q) GO TO 40 33 CALL QINT(X(I),X(I+1),1. ,V2A,V5A,V3A,NA, * 1. ,V2B,V5B,V3B,NB,SIGMAP,Q) GO TO 40 34 CALL QINT(X(I),X(I+1),1. ,V2A,V5A,V3A,NA, * V5B,V3B,0. ,0. ,NB,SIGMAP,Q) C QINT C 40 I = IA+1 IF(I.GE.N) GO TO 90 J = I-IB+3 IF(J.LT.1) GO TO 90 IF(J.GT.4) GO TO 90 GO TO (41,42,43,44),J 41 CALL QINT(X(I),X(I+1),V5A,V3A,0. ,0. ,NA, * 0. ,0. ,V4B,V1B,NB,SIGMAP,Q) GO TO 90 42 CALL QINT(X(I),X(I+1),V5A,V3A,0. ,0. ,NA, * V4B,V1B,1. ,V2B,NB,SIGMAP,Q) GO TO 90 43 CALL QINT(X(I),X(I+1),V5A,V3A,0. ,0. ,NA, * 1. ,V2B,V5B,V3B,NB,SIGMAP,Q) GO TO 90 44 CALL QINT(X(I),X(I+1),V5A,V3A,0. ,0. ,NA, * V5B,V3B,0. ,0. ,NB,SIGMAP,Q) C QINT C 90 CONTINUE RETURN END C C======================================================================= C C C SUBROUTINE QINT(X1,X2,FA1,DA1,FA2,DA2,NA, * FB1,DB1,FB2,DB2,NB,SIGMAP,Q) C INTEGER NA,NB REAL X1,X2,FA1,DA1,FA2,DA2,FB1,DB1,FB2,DB2,SIGMAP,Q C C Complement to FITPACK C by Alan Kaylor Cline C coded -- January 23, 1994 C by Ludek Klimes C Department of Geophysics C Charles University, Prague C C This subroutine evaluates the integral of the product C of the given derivatives of the two given cubic functions C or spline under tension basis functions in one variable, C over a single specified interval. C C On input-- C C X1, X2 endpoints of the given interval. C C FA1, DA1 function value and second derivative of the C first given function at X1. C C FA2, DA2 function value and second derivative of the C first given function at X2. C C NA is the order of the partial derivative of C the first argument function in the scalar product. C C FB1, DB1, FB2, DB2 the same as FA1, DA1, FA2, DA2, but C for the second given function. C C NB is the order of the partial derivative of C the second argument function in the scalar product. C C SIGMAP is the denormalized tension factor. C C And C C None of the input parameters are altered. C C On output-- C C Q is the integral of the product of NA-th partial C derivative of the first function and C NB-th partial derivative of the second function, C over the interval X1,X2. C C This subroutine references package module SNHCSH. C EXTERNAL SNHCSH C C--------------------------------------------------------------- C C Other variables used inside the subroutine QINT: C INTEGER MA,MB,M REAL QQ,H,SH,CH,SH1,CH1,SIGMA2 REAL A1,A2,A3,A4,B1,B2,B3,B4,AB11,AB21,AB12,AB22 C C--------------------------------------------------------------- C MA= MOD(NA,2) MB= MOD(NB,2) M = MA+MA+MB+1 QQ= 0. C IF(SIGMAP.NE.0.) GO TO 40 C C No tension: H = X2-X1 IF(NA.LE.3.AND.NB.LE.3) GO TO 1 GO TO 91 1 IF(NA.LE.1) GO TO 3 C Coefficients of linear function A3= DA2/H A4=-DA1/H IF(NB.LE.1) GO TO 2 C Coefficients of linear function B3= DB2/H B4=-DB1/H GO TO 80 2 CONTINUE C Coefficients of cubic and linear functions B1= DB2/H B2=-DB1/H B3= FB2/H-DB2*H/6. B4=-FB1/H+DB1*H/6. GO TO 30 3 CONTINUE C Coefficients of cubic and linear functions A1= DA2/H A2=-DA1/H A3= FA2/H-DA2*H/6. A4=-FA1/H+DA1*H/6. IF(NB.LE.1) GO TO 4 C Coefficients of linear function B3= DB2/H B4=-DB1/H GO TO 20 4 CONTINUE C Coefficients of cubic and linear functions B1= DB2/H B2=-DB1/H B3= FB2/H-DB2*H/6. B4=-FB1/H+DB1*H/6. C C Integrals of (cubic function)*(cubic function): GO TO (11,12,13,14),M C (even derivative)*(even derivative) 11 AB11= (H**7)/252. AB21=-(H**7)/5040. AB12= AB21 AB22= AB11 GO TO 15 C (even derivative)*(odd derivative) 12 AB11= (H**6)/72. AB21=-(H**6)/720. AB12=-AB21 AB22=-AB11 GO TO 15 C (odd derivative)*(even derivative) 13 AB11= (H**6)/72. AB21= (H**6)/720. AB12=-AB21 AB22=-AB11 GO TO 15 C (odd derivative)*(odd derivative) 14 AB11= (H**5)/20. AB21= (H**5)/120. AB12= AB21 AB22= AB11 C Accumulation of the computed integral: 15 QQ=QQ+A1*(AB11*B1+AB12*B2)+A2*(AB21*B1+AB22*B2) C C Integrals of (cubic function)*(linear function): 20 GO TO (21,22,23,24),M C (even derivative)*(even derivative) 21 AB11= (H**5)/30. AB21=-(H**5)/120. AB12= AB21 AB22= AB11 GO TO 25 C (even derivative)*(odd derivative) 22 AB11= (H**4)/24. AB21=-(H**4)/24. AB12=-AB21 AB22=-AB11 GO TO 25 C (odd derivative)*(even derivative) 23 AB11= (H**4)/8. AB21= (H**4)/24. AB12=-AB21 AB22=-AB11 GO TO 25 C (odd derivative)*(odd derivative) 24 AB11= (H**3)/6. AB21= (H**3)/6. AB12= AB21 AB22= AB11 C Accumulation of the computed integral: 25 QQ=QQ+A1*(AB11*B3+AB12*B4)+A2*(AB21*B3+AB22*B4) C IF(NB.GT.1) GO TO 80 C C Integrals of (linear function)*(cubic function): 30 GO TO (31,32,33,34),M C (even derivative)*(even derivative) 31 AB11= (H**5)/30. AB21=-(H**5)/120. AB12= AB21 AB22= AB11 GO TO 35 C (even derivative)*(odd derivative) 32 AB11= (H**4)/8. AB21=-(H**4)/24. AB12=-AB21 AB22=-AB11 GO TO 35 C (odd derivative)*(even derivative) 33 AB11= (H**4)/24. AB21= (H**4)/24. AB12=-AB21 AB22=-AB11 GO TO 35 C (odd derivative)*(odd derivative) 34 AB11= (H**3)/6. AB21= (H**3)/6. AB12= AB21 AB22= AB11 C Accumulation of the computed integral: 35 QQ=QQ+A3*(AB11*B1+AB12*B2)+A4*(AB21*B1+AB22*B2) GO TO 80 C C C Nonzero tension: 40 H = SIGMAP*(X2-X1) CALL SNHCSH(SH1,CH1,H,0) C SNHCSH SH= SH1+H CH= CH1+1. SIGMA2= SIGMAP*SIGMAP C C Coefficients of hyperbolic functions (multiplied by SH) A1= DA2/SIGMA2 A2=-DA1/SIGMA2 B1= DB2/SIGMA2 B2=-DB1/SIGMA2 C C Doubled C integrals of (hyperbolic function)*(hyperbolic function): GO TO (51,52,53,54),M C (even derivative)*(even derivative) 51 AB11= CH*SH1+H*CH1 AB21= SH1-H*CH1 AB12= AB21 AB22= AB11 GO TO 55 C (even derivative)*(odd derivative) 52 AB11= SH*SH AB21= -H*SH AB12=-AB21 AB22=-AB11 GO TO 55 C (odd derivative)*(even derivative) 53 AB11= SH*SH AB21= H*SH AB12=-AB21 AB22=-AB11 GO TO 55 C (odd derivative)*(odd derivative) 54 AB11= SH*CH+H AB21= SH+H*CH AB12= AB21 AB22= AB11 C Accumulation of the computed integral: 55 QQ=QQ+(A1*(AB11*B1+AB12*B2)+A2*(AB21*B1+AB22*B2))/(2.*SH*SH) C IF(NB.GT.1) GO TO 70 C C Coefficients of linear function B3= ( FB2-B1)/H B4= (-FB1-B2)/H C C Integrals of (hyperbolic function)*(linear function): GO TO (61,62,63,64),M C (even derivative)*(even derivative) 61 AB11= H*CH1-SH1 AB21= -SH1 AB12= AB21 AB22= AB11 GO TO 65 C (even derivative)*(odd derivative) 62 AB11= CH1 AB21=-CH1 AB12=-AB21 AB22=-AB11 GO TO 65 C (odd derivative)*(even derivative) 63 AB11= H*SH-CH1 AB21= CH1 AB12=-AB21 AB22=-AB11 GO TO 65 C (odd derivative)*(odd derivative) 64 AB11= SH AB21= SH AB12= AB21 AB22= AB11 C Accumulation of the computed integral: 65 QQ=QQ+(A1*(AB11*B3+AB12*B4)+A2*(AB21*B3+AB22*B4))/SH C 70 IF(NA.GT.1) GO TO 90 C C Coefficients of linear function A3= ( FA2-A1)/H A4= (-FA1-A2)/H C C Integrals of (linear function)*(hyperbolic function): GO TO (71,72,73,74),M C (even derivative)*(even derivative) 71 AB11= H*CH1-SH1 AB21= -SH1 AB12= AB21 AB22= AB11 GO TO 75 C (even derivative)*(odd derivative) 72 AB11= H*SH-CH1 AB21= -CH1 AB12=-AB21 AB22=-AB11 GO TO 75 C (odd derivative)*(even derivative) 73 AB11= CH1 AB21= CH1 AB12=-AB21 AB22=-AB11 GO TO 75 C (odd derivative)*(odd derivative) 74 AB11= SH AB21= SH AB12= AB21 AB22= AB11 C Accumulation of the computed integral: 75 QQ=QQ+(A3*(AB11*B1+AB12*B2)+A4*(AB21*B1+AB22*B2))/SH C IF(NB.GT.1) GO TO 90 C C Integrals of (linear function)*(linear function): 80 GO TO (81,82,83,84),M C (even derivative)*(even derivative) 81 AB11= (H**3)/3. AB21=-(H**3)/6. AB12= AB21 AB22= AB11 GO TO 85 C (even derivative)*(odd derivative) 82 AB11= (H**2)/2. AB21=-AB11 AB12=-AB21 AB22=-AB11 GO TO 85 C (odd derivative)*(even derivative) 83 AB11= (H**2)/2. AB21= AB11 AB12=-AB21 AB22=-AB11 GO TO 85 C (odd derivative)*(odd derivative) 84 AB11= H AB21= H AB12= AB21 AB22= AB11 C Accumulation of the computed integral: 85 QQ=QQ+A3*(AB11*B3+AB12*B4)+A4*(AB21*B3+AB22*B4) C C Transformation from independent variable SIGMAP*X to X 90 IF(SIGMAP.NE.0.) QQ=QQ*SIGMAP**(NA+NB-1) 91 Q= Q+QQ RETURN END C C======================================================================= C