      SUBROUTINE KURVP1 (N,X,Y,XP,YP,TEMP,S,SIGMA,IERR)
C
      INTEGER N,IERR
      REAL X(N),Y(N),XP(N),YP(N),TEMP(1),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 FORMING A CLOSED CURVE IN
C THE PLANE AND PASSING THROUGH A SEQUENCE OF PAIRS
C (X(1),Y(1)),...,(X(N),Y(N)). FOR ACTUAL COMPUTATION OF
C POINTS ON THE CURVE IT IS NECESSARY TO CALL THE SUBROUTINE
C KURVP2.
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. (ADJACENT X-Y PAIRS MUST BE DISTINCT, I. E.
C   EITHER X(I) .NE. X(I+1) OR Y(I) .NE. Y(I+1), FOR
C   I = 1,...,N-1.)
C
C   XP AND YP ARE ARRAYS OF LENGTH AT LEAST N.
C
C   TEMP IS AN ARRAY OF LENGTH AT LEAST 2*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 AND YP CONTAIN INFORMATION ABOUT THE CURVATURE OF THE
C   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 COORDINATE PAIRS COINCIDE.
C
C AND
C
C   N, X, Y, AND SIGMA ARE UNALTERED,
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES TERMS AND
C SNHCSH.
C
C-----------------------------------------------------------
C
      NM1 = N-1
      NP1 = N+1
      IERR = 0
      IF (N .LE. 1) GO TO 7
C
C DETERMINE POLYGONAL ARCLENGTHS
C
      S(1) = SQRT((X(N)-X(1))**2+(Y(N)-Y(1))**2)
      DO 1 I = 2,N
        IM1 = I-1
    1   S(I) = S(IM1)+SQRT((X(I)-X(IM1))**2+
     *         (Y(I)-Y(IM1))**2)
C
C DENORMALIZE TENSION FACTOR
C
      SIGMAP = ABS(SIGMA)*FLOAT(N)/S(N)
C
C SET UP RIGHT HAND SIDES OF TRIDIAGONAL (WITH CORNER
C ELEMENTS) LINEAR SYSTEM FOR XP AND YP
C
      DELS1 = S(1)
      IF (DELS1 .EQ. 0.) GO TO 8
      DX1 = (X(1)-X(N))/DELS1
      DY1 = (Y(1)-Y(N))/DELS1
      CALL TERMS(DIAG1,SDIAG1,SIGMAP,DELS1)
      DELS2 = S(2)-S(1)
      IF (DELS2 .EQ. 0.) GO TO 8
      DX2 = (X(2)-X(1))/DELS2
      DY2 = (Y(2)-Y(1))/DELS2
      CALL TERMS(DIAG2,SDIAG2,SIGMAP,DELS2)
      DIAG = DIAG1+DIAG2
      DIAGIN = 1./DIAG
      XP(1) = (DX2-DX1)*DIAGIN
      YP(1) = (DY2-DY1)*DIAGIN
      TEMP(NP1) = -SDIAG1*DIAGIN
      TEMP(1) = SDIAG2*DIAGIN
      DX1 = DX2
      DY1 = DY2
      DIAG1 = DIAG2
      SDIAG1 = SDIAG2
      IF (N .EQ. 2) GO TO 3
      DO 2 I = 2,NM1
        NPI = N+I
        DELS2 = S(I+1)-S(I)
        IF (DELS2 .EQ. 0.) GO TO 8
        DX2 = (X(I+1)-X(I))/DELS2
        DY2 = (Y(I+1)-Y(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
        TEMP(NPI) = -TEMP(NPI-1)*SDIAG1*DIAGIN
        TEMP(I) = SDIAG2*DIAGIN
        DX1 = DX2
        DY1 = DY2
        DIAG1 = DIAG2
    2   SDIAG1 = SDIAG2
    3 DELS2 = S(1)
      DX2 = (X(1)-X(N))/DELS2
      DY2 = (Y(1)-Y(N))/DELS2
      CALL TERMS(DIAG2,SDIAG2,SIGMAP,DELS2)
      XP(N) = DX2-DX1
      YP(N) = DY2-DY1
      TEMP(NM1) = TEMP(2*N-1)-TEMP(NM1)
      IF (N.EQ.2) GO TO 5
C
C PERFORM FIRST STEP OF BACK SUBSTITUTION
C
      DO 4 I = 3,N
        IBAK = NP1-I
        NPIBAK = N+IBAK
        XP(IBAK) = XP(IBAK)-TEMP(IBAK)*XP(IBAK+1)
        YP(IBAK) = YP(IBAK)-TEMP(IBAK)*YP(IBAK+1)
    4   TEMP(IBAK) = TEMP(NPIBAK)-TEMP(IBAK)*TEMP(IBAK+1)
    5 XP(N) = (XP(N)-SDIAG2*XP(1)-SDIAG1*XP(NM1))/
     *        (DIAG1+DIAG2+SDIAG2*TEMP(1)+SDIAG1*TEMP(NM1))
      YP(N) = (YP(N)-SDIAG2*YP(1)-SDIAG1*YP(NM1))/
     *        (DIAG1+DIAG2+SDIAG2*TEMP(1)+SDIAG1*TEMP(NM1))
C
C PERFORM SECOND STEP OF BACK SUBSTITUTION
C
      XPN = XP(N)
      YPN = YP(N)
      DO 6 I = 1,NM1
        XP(I) = XP(I)+TEMP(I)*XPN
    6   YP(I) = YP(I)+TEMP(I)*YPN
      RETURN
C
C TOO FEW POINTS
C
    7 IERR = 1
      RETURN
C
C COINCIDENT ADJACENT POINTS
C
    8 IERR = 2
      RETURN
      END
