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