      SUBROUTINE QURV1 (N,X,Y,Z,SLP1X,SLP1Y,SLP1Z,SLPNX,
     *                  SLPNY,SLPNZ,ISLPSW,XP,YP,ZP,TEMP,
     *                  S,SIGMA,IERR)
C
      INTEGER N,ISLPSW,IERR
      REAL X(N),Y(N),Z(N),SLP1X,SLP1Y,SLP1Z,SLPNX,SLPNY,
     *     SLPNZ,XP(N),YP(N),ZP(N),TEMP(N),S(N),SIGMA
C
C                      FROM THE SPLINE UNDER TENSION PACKAGE
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 DETERMINES THE PARAMETERS NECESSARY TO
C COMPUTE A SPLINE UNDER TENSION PASSING THROUGH A SEQUENCE
C OF TRIPLES (X(1),Y(1),Z(1)),...,(X(N),Y(N),Z(N)). THE
C SLOPES AT THE TWO ENDS OF THE CURVE MAY BE SPECIFIED OR
C OMITTED. FOR ACTUAL COMPUTATION OF POINTS ON THE CURVE
C IT IS NECESSARY TO CALL THE SUBROUTINE QURV2.
C
C ON INPUT--
C
C   N IS THE NUMBER OF POINTS TO BE INTERPOLATED (N.GE.2).
C
C   X IS AN ARRAY CONTAINING THE N X-COORDINATES OF THE
C   POINTS.
C
C   Y IS AN ARRAY CONTAINING THE N Y-COORDINATES OF THE
C   POINTS.
C
C   Z IS AN ARRAY CONTAINING THE N Z-COORDINATES OF THE
C   POINTS. (ADJACENT X-Y-Z TRIPLES MUST BE DISTINCT, I. E.
C   EITHER X(I) .NE. X(I+1) OR Y(I) .NE. Y(I+1) OR Z(I) .NE.
C   Z(I+1), FOR I = 1,...,N-1 ).
C
C   SLP1X, SLP1Y, SLP1Z AND SLPNX, SLPNY, SLPNZ CONTAIN THE
C   DESIRED VALUES OF THE COMPONENTS OF TANGENT VECTORS TO
C   THE CURVE AT (X(1),Y(1),Z(1)) AND(X(N),Y(N),Z(N)),
C   RESPECTIVELY. THE POSITIVE SENSE OF THE CURVE IS ASSUMED
C   TO BE THAT MOVING FROM POINT 1 TO POINT N. THE USER MAY
C   OMIT VALUES FOR EITHER OR BOTH OF THESE TRIPLES AND
C   SIGNAL THIS WITH ISLPSW.
C
C   ISLPSW CONTAINS A SWITCH INDICATING WHICH SLOPE DATA
C   SHOULD BE USED AND WHICH SHOULD BE ESTIMATED BY THIS
C   SUBROUTINE,
C          = 0 IF SLP1X, SLP1Y, SLP1Z AND SLPNX, SLPNY,
C              SLPNZ ARE TO BE USED,
C          = 1 IF SLP1X, SLP1Y, SLP1Z ARE TO BE USED BUT
C              NOT SLPNX, SLPNY, SLPNZ,
C          = 2 IF SLPNX, SLPNY, SLPNZ ARE TO BE USED BUT
C              NOT SLP1X, SLP1Y, SLP1Z,
C          = 3 IF BOTH END-TANGENTS ARE TO BE ESTIMATED
C              INTERNALLY.
C
C   XP, YP, AND ZP ARE ARRAYS OF LENGTH AT LEAST N.
C
C   TEMP IS AN ARRAY OF LENGTH AT LEAST N WHICH IS USED
C   FOR SCRATCH STORAGE.
C
C   S IS AN ARRAY OF LENGTH AT LEAST N.
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 CUBIC
C   SPLINE. IF ABS(SIGMA) IS LARGE (E. G. 50.) THE RESULTING
C   CURVE IS NEARLY A POLYGONAL LINE. IF SIGMA EQUALS ZERO A
C   CUBIC SPLINE RESULTS. A STANDARD VALUE FOR SIGMA IS
C   APPROXIMATELY 1. IN ABSOLUTE VALUE.
C
C ON OUTPUT--
C
C   XP, YP, AND ZP CONTAIN INFORMATION ABOUT THE CURVATURE
C   OF THE CURVE AT THE GIVEN NODES.
C
C   S CONTAINS THE POLYGONAL ARCLENGTHS OF THE CURVE.
C
C   IERR CONTAINS AN ERROR FLAG,
C        = 0 FOR NORMAL RETURN,
C        = 1 IF N IS LESS THAN 2,
C        = 2 IF ADJACENT TRIPLES COINCIDE.
C
C AND
C
C   N, X, Y, Z, SLP1X, SLP1Y, SLP1Z, SLPNX, SLPNY, SLPNZ,
C   ISLPSW, AND SIGMA ARE UNALTERED,
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES CEEZ, TERMS,
C AND SNHCSH.
C
C-----------------------------------------------------------
C
      NM1 = N-1
      NP1 = N+1
      IERR = 0
      IF (N .LE. 1) GO TO 9
C
C DETERMINE POLYGONAL ARCLENGTHS
C
      S(1) = 0.
      DO 1 I = 2,N
        IM1 = I-1
    1   S(I) = S(IM1)+SQRT((X(I)-X(IM1))**2+(Y(I)-Y(IM1))**2
     *       +(Z(I)-Z(IM1))**2)
C
C DENORMALIZE TENSION FACTOR
C
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S(N)
C
C APPROXIMATE END SLOPES
C
      IF (ISLPSW .GE. 2) GO TO 2
      SLPP1X = SLP1X
      SLPP1Y = SLP1Y
      SLPP1Z = SLP1Z
      GO TO 3
    2 DELS1 = S(2)-S(1)
      DELS2 = DELS1+DELS1
      IF (N .GT. 2) DELS2 = S(3)-S(1)
      IF (DELS1 .EQ. 0. .OR. DELS2 .EQ. 0.) GO TO 9
      CALL CEEZ (DELS1,DELS2,SIGMAP,C1,C2,C3,N)
      SLPP1X = C1*X(1)+C2*X(2)
      SLPP1Y = C1*Y(1)+C2*Y(2)
      SLPP1Z = C1*Z(1)+C2*Z(2)
      IF (N .EQ. 2) GO TO 3
      SLPP1X = SLPP1X+C3*X(3)
      SLPP1Y = SLPP1Y+C3*Y(3)
      SLPP1Z = SLPP1Z+C3*Z(3)
    3 DELT = SQRT(SLPP1X*SLPP1X+SLPP1Y*SLPP1Y+SLPP1Z*SLPP1Z)
      SLPP1X = SLPP1X/DELT
      SLPP1Y = SLPP1Y/DELT
      SLPP1Z = SLPP1Z/DELT
      IF (ISLPSW .EQ. 1 .OR. ISLPSW .EQ. 3) GO TO 4
      SLPPNX = SLPNX
      SLPPNY = SLPNY
      SLPPNZ = SLPNZ
      GO TO 5
    4 DELSN = S(N)-S(NM1)
      DELSNM = DELSN+DELSN
      IF (N .GT. 2) DELSNM = S(N)-S(N-2)
      IF (DELSN .EQ. 0. .OR. DELSNM .EQ. 0.) GO TO 10
      CALL CEEZ (-DELSN,-DELSNM,SIGMAP,C1,C2,C3,N)
      SLPPNX = C1*X(N)+C2*X(NM1)
      SLPPNY = C1*Y(N)+C2*Y(NM1)
      SLPPNZ = C1*Z(N)+C2*Z(NM1)
      IF (N .EQ. 2) GO TO 5
      SLPPNX = SLPPNX+C3*X(N-2)
      SLPPNY = SLPPNY+C3*Y(N-2)
      SLPPNZ = SLPPNZ+C3*Z(N-2)
    5 DELT = SQRT(SLPPNX*SLPPNX+SLPPNY*SLPPNY+SLPPNZ*SLPPNZ)
      SLPPNX = SLPPNX/DELT
      SLPPNY = SLPPNY/DELT
      SLPPNZ = SLPPNZ/DELT
C
C SET UP RIGHT HAND SIDES AND TRIDIAGONAL SYSTEM FOR XP, YP
C AND ZP AND PERFORM FORWARD ELIMINATION
C
      DX1 = (X(2)-X(1))/S(2)
      DY1 = (Y(2)-Y(1))/S(2)
      DZ1 = (Z(2)-Z(1))/S(2)
      CALL TERMS (DIAG1,SDIAG1,SIGMAP,S(2))
      XP(1) = (DX1-SLPP1X)/DIAG1
      YP(1) = (DY1-SLPP1Y)/DIAG1
      ZP(1) = (DZ1-SLPP1Z)/DIAG1
      TEMP(1) = SDIAG1/DIAG1
      IF (N .EQ. 2) GO TO 7
      DO 6 I = 2,NM1
        DELS2 = S(I+1)-S(I)
        IF (DELS2 .EQ. 0.) GO TO 10
        DX2 = (X(I+1)-X(I))/DELS2
        DY2 = (Y(I+1)-Y(I))/DELS2
        DZ2 = (Z(I+1)-Z(I))/DELS2
        CALL TERMS (DIAG2,SDIAG2,SIGMAP,DELS2)
        DIAG = DIAG1+DIAG2-SDIAG1*TEMP(I-1)
        DIAGIN = 1./DIAG
        XP(I) = (DX2-DX1-SDIAG1*XP(I-1))*DIAGIN
        YP(I) = (DY2-DY1-SDIAG1*YP(I-1))*DIAGIN
        ZP(I) = (DZ2-DZ1-SDIAG1*ZP(I-1))*DIAGIN
        TEMP(I) = SDIAG2*DIAGIN
        DX1 = DX2
        DY1 = DY2
        DZ1 = DZ2
        DIAG1 = DIAG2
    6   SDIAG1 = SDIAG2
    7 DIAG = DIAG1-SDIAG1*TEMP(NM1)
      XP(N) = (SLPPNX-DX1-SDIAG1*XP(NM1))/DIAG
      YP(N) = (SLPPNY-DY1-SDIAG1*YP(NM1))/DIAG
      ZP(N) = (SLPPNZ-DZ1-SDIAG1*ZP(NM1))/DIAG
C
C PERFORM BACK SUBSTITUTION
C
      DO 8 I = 2,N
        IBAK = NP1-I
        T = TEMP(IBAK)
        XP(IBAK) = XP(IBAK)-T*XP(IBAK+1)
        YP(IBAK) = YP(IBAK)-T*YP(IBAK+1)
    8   ZP(IBAK) = ZP(IBAK)-T*ZP(IBAK+1)
      RETURN
C
C TOO FEW POINTS
C
    9 IERR = 1
      RETURN
C
C COINCIDENT ADJACENT POINTS
C
   10 IERR = 2
      RETURN
      END
