C     ALGORITHM 476 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN COMM. ACM, VOL. 17, NO. 04,
C     P. 220.
      SUBROUTINE CURV1(N, X, Y, SLP1, SLPN, YP, TEMP, SIGMA)               A  10
      INTEGER N
      REAL X(N), Y(N), SLP1, SLPN, YP(N), TEMP(N), SIGMA
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO
C COMPUTE AN INTERPOLATORY SPLINE UNDER TENSION THROUGH
C A SEQUENCE OF FUNCTIONAL VALUES. THE SLOPES AT THE TWO
C ENDS OF THE CURVE MAY BE SPECIFIED OR OMITTED. FOR ACTUAL
C COMPUTATION OF POINTS ON THE CURVE IT IS NECESSARY TO CALL
C THE FUNCTION CURV2.
C ON INPUT--
C N IS THE NUMBER OF VALUES TO BE INTERPOLATED (N.GE.2),
C X IS AN ARRAY OF THE N INCREASING ABSCISSAE OF THE
C FUNCTIONAL VALUES,
C Y IS AN ARRAY OF THE N ORDINATES OF THE VALUES,(I.E.Y(K)
C IS THE FUNCTIONAL VALUE CORRESPONDING TO X(K)),
C SLP1 AND SLPN CONTAIN THE DESIRED VALUES FOR THE FIRST
C DERIVATIVE OF THE CURVE AT X(1) AND X(N), RESPECTIVELY.
C IF THE QUANTITY SIGMA IS NEGATIVE THESE VALUES WILL BE
C DETERMINED INTERNALLY AND THE USER NEED ONLY FURNISH
C PLACE-HOLDING PARAMETERS FOR SLP1 AND SLPN. SUCH PLACE-
C HOLDING PARAMETERS WILL BE IGNORED BUT NOT DESTROYED,
C YP IS AN ARRAY OF LENGTH AT LEAST N
C TEMP IS AN ARRAY OF LENGTH AT LEAST N WHICH IS USED FOR
C SCRATCH STORAGE,
C AND
C SIGMA CONTAINS THE TENSION FACTOR. THIS IS NON-ZERO AND
C INDICATES THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY
C ZERO (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. THE SIGN
C OF SIGMA INDICATES WHETHER THE DERIVATIVE INFORMATION
C HAS BEEN INPUT OR NOT. IF SIGMA IS NEGATIVE THE ENDPOINT
C DERIVATIVES WILL BE DETERMINED INTERNALLY. A STANDARD
C VALUE FOR SIGMA IS APPROXIMATELY 1. IN ABSOLUTE VALUE.
C ON OUTPUT--
C YP CONTAINS VALUES PROPORTIONAL TO THE SECOND DERIVATIVE
C OF THE CURVE AT THE GIVEN NODES.
C N,X,Y,SLP1,SLPN AND SIGMA ARE UNALTERED,
      NM1 = N - 1
      NP1 = N + 1
      DELX1 = X(2) - X(1)
      DX1 = (Y(2)-Y(1))/DELX1
C DETERMINE SLOPES IF NECESSARY
      IF (SIGMA.LT.0.) GO TO 50
      SLPP1 = SLP1
      SLPPN = SLPN
C DENORMALIZE TENSION FACTOR
   10 SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
C SET UP RIGHT HAND SIDE AND TRIDIAGONAL SYSTEM FOR YP AND
C PERFORM FORWARD ELIMINATION
      DELS = SIGMAP*DELX1
      EXPS = EXP(DELS)
      SINHS = .5*(EXPS-1./EXPS)
      SINHIN = 1./(DELX1*SINHS)
      DIAG1 = SINHIN*(DELS*.5*(EXPS+1./EXPS)-SINHS)
      DIAGIN = 1./DIAG1
      YP(1) = DIAGIN*(DX1-SLPP1)
      SPDIAG = SINHIN*(SINHS-DELS)
      TEMP(1) = DIAGIN*SPDIAG
      IF (N.EQ.2) GO TO 30
      DO 20 I=2,NM1
        DELX2 = X(I+1) - X(I)
        DX2 = (Y(I+1)-Y(I))/DELX2
        DELS = SIGMAP*DELX2
        EXPS = EXP(DELS)
        SINHS = .5*(EXPS-1./EXPS)
        SINHIN = 1./(DELX2*SINHS)
        DIAG2 = SINHIN*(DELS*(.5*(EXPS+1./EXPS))-SINHS)
        DIAGIN = 1./(DIAG1+DIAG2-SPDIAG*TEMP(I-1))
        YP(I) = DIAGIN*(DX2-DX1-SPDIAG*YP(I-1))
        SPDIAG = SINHIN*(SINHS-DELS)
        TEMP(I) = DIAGIN*SPDIAG
        DX1 = DX2
        DIAG1 = DIAG2
   20 CONTINUE
   30 DIAGIN = 1./(DIAG1-SPDIAG*TEMP(NM1))
      YP(N) = DIAGIN*(SLPPN-DX2-SPDIAG*YP(NM1))
C PERFORM BACK SUBSTITUTION
      DO 40 I=2,N
        IBAK = NP1 - I
        YP(IBAK) = YP(IBAK) - TEMP(IBAK)*YP(IBAK+1)
   40 CONTINUE
      RETURN
   50 IF (N.EQ.2) GO TO 60
C IF NO DERIVATIVES ARE GIVEN USE SECOND ORDER POLYNOMIAL
C INTERPOLATION ON INPUT DATA FOR VALUES AT ENDPOINTS.
      DELX2 = X(3) - X(2)
      DELX12 = X(3) - X(1)
      C1 = -C(DELX12+DELX1)/DELX12/DELX1
      C2 = DELX12/DELX1/DELX2
      C3 = -DELX1/DELX12/DELX2
      SLPP1 = C1*Y(1) + C2*Y(2) + C3*Y(3)
      DELN = X(N) - X(NM1)
      DELNM1 = X(NM1) - X(N-2)
      DELNN = X(N) - X(N-2)
      C1 = (DELNN+DELN)/DELNN/DELN
      C2 = -DELNN/DELN/DELNM1
      C3 = DELN/DELNN/DELNM1
      SLPPN = C3*Y(N-2) + C2*Y(NM1) + C1*Y(N)
      GO TO 10
C IF ONLY TWO POINTS AND NO DERIVATIVES ARE GIVEN, USE
C STRAIGHT LINE FOR CURVE
   60 YP(1) = 0.
      YP(2) = 0.
      RETURN
      END
      FUNCTION CURV2(T, N, X, Y, YP, SIGMA, IT)                            B  10
      INTEGER N, IT
      REAL T, X(N), Y(N), YP(N), SIGMA
C THIS FUNCTION INTERPOLATES A CURVE AT A GIVEN POINT
C USING A SPLINE UNDER TENSION. THE SUBROUTINE CURV1 SHOULD
C BE CALLED EARLIER TO DETERMINE CERTAIN NECESSARY
C PARAMETERS.
C ON INPUT--
C T CONTAINS A REAL VALUE TO BE MAPPED ONTO THE INTERPO-
C LATING CURVE.
C N CONTAINS THE NUMBER OF POINTS WHICH WERE INTERPOLATED
C TO DETERMINE THE CURVE,
C X AND Y ARE ARRAYS CONTAINING THE ORDINATES AND ABCISSAS
C OF THE INTERPOLATED POINTS,
C YP IS AN ARRAY WITH VALUES PROPORTIONAL TO THE SECOND
C DERIVATIVE OF THE CURVE AT THE NODES
C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED)
C IT IS AN INTEGER SWITCH. IF IT IS NOT 1 THIS INDICATES
C THAT THE FUNCTION HAS BEEN CALLED PREVIOUSLY (WITH N,X,
C Y,YP, AND SIGMA UNALTERED) AND THAT THIS VALUE OF T
C EXCEEDS THE PREVIOUS VALUE. WITH SUCH INFORMATION THE
C FUNCTION IS ABLE TO PERFORM THE INTERPOLATION MUCH MORE
C RAPIDLY. IF A USER SEEKS TO INTERPOLATE AT A SEQUENCE
C OF POINTS, EFFICIENCY IS GAINED BY ORDERING THE VALUES
C INCREASING AND SETTING IT TO THE INDEX OF THE CALL.
C IF IT IS 1 THE SEARCH FOR THE INTERVAL (X(K),X(K+1))
C CONTAINING T STARTS WITH K=1.
C THE PARAMETERS N,X,Y,YP AND SIGMA SHOULD BE INPUT
C UNALTERED FROM THE OUTPUT OF CURV1.
C ON OUTPUT--
C CURV2 CONTAINS THE INTERPOLATED VALUE. FOR T LESS THAN
C X(1) CURV2 = Y(1). FOR T GREATER THAN X(N) CURV2 = Y(N).
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
      S = X(N) - X(1)
C DENORMALIZE SIGMA
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S
C IF IT.NE.1 START SEARCH WHERE PREVIOUSLY TERMINATED,
C OTHERWISE START FROM BEGINNING
      IF (IT.EQ.1) I1 = 2
C SEARCH FOR INTERVAL
   10 DO 20 I=I1,N
        IF (X(I)-T) 20, 20, 30
   20 CONTINUE
      I = N
C CHECK TO INSURE CORRECT INTERVAL
   30 IF (X(I-1).LE.T .OR. T.LE.X(1)) GO TO 40
C RESTART SEARCH  AND RESET I1
C ( INPUT ""IT"" WAS INCORRECT )
      I1 = 2
      GO TO 10
C SET UP AND PERFORM INTERPOLATION
   40 DEL1 = T - X(I-1)
      DEL2 = X(I) - T
      DELS = X(I) - X(I-1)
      EXPS1 = EXP(SIGMAP*DEL1)
      SINHD1 = .5*(EXPS1-1./EXPS1)
      EXPS = EXP(SIGMAP*DEL2)
      SINHD2 = .5*(EXPS-1./EXPS)
      EXPS = EXPS1*EXPS
      SINHS = .5*(EXPS-1./EXPS)
      CURV2 = (YP(I)*SINHD1+YP(I-1)*SINHD2)/SINHS +
     * ((Y(I)-YP(I))*DEL1+(Y(I-1)-YP(I-1))*DEL2)/DELS
      I1 = I
      RETURN
      END
      SUBROUTINE KURV1(N, X, Y, SLP1, SLPN, XP, YP, TEMP, S,               C  10
     * SIGMA)
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO
C COMPUTE A SPLINE UNDER TENSION PASSING THROUGH A SEQUENCE
C OF PAIRS  (X(1),Y(1)),...,(X(N),Y(N)) IN THE PLANE. THE
C SLOPES AT THE TWO ENDS OF THE CURVE MAY BE SPECIFIED OR
C OMITTED. FOR ACTUAL COMPUTATION OF POINTS ON THE CURVE IT
C IS NECESSARY TO CALL THE SUBROUTINE KURV2.
C ON INPUT--
C N IS THE NUMBER OF POINTS TO BE INTERPOLATED (N.GE.2),
C X IS AN ARRAY CONTAINING THE N X-COORDINATES OF THE
C POINTS,
C Y IS AN ARRAY CONTAINING THE N Y-COORDINATES OF THE
C POINTS,
C SLP1 AND SLPN CONTAIN THE DESIRED VALUES FOR THE SLOPE
C OF THE CURVE AT  (X(1),Y(1)) AND  (X(N),Y(N)), RESPEC-
C TIVELY. THESE QUANTITIES ARE IN DEGREES AND MEASURED
C COUNTERCLOCKWISE FROM THE POSITIVE X-AXIS. THE POSITIVE
C SENSE OF THE CURVE IS ASSUMED TO BE THAT MOVING FROM THE
C POINT 1 TO POINT N. IF THE QUANTITY SIGMA IS NEGATIVE
C THESE SLOPES WILL BE DETERMINED INTERNALLY AND THE USER
C NEED ONLY FURNISH PLACE-HOLDING PARAMETERS FOR SLP1 AND
C SLPN. SUCH PLACE-HOLDING PARAMETERS WILL BE IGNORED BUT
C NOT DESTROYED,
C XP,YP ARE ARRAYS OF LENGTH AT LEAST N,
C TEMP IS AN ARRAY OF LENGTH AT LEAST N WHICH IS USED FOR
C SCRATCH STORAGE,
C AND
C SIGMA CONTAINS THE TENSION FACTOR. THIS IS NON-ZERO AND
C INDICATES THE CURVINESS DESIRED. IF ABS(SIGMA) IS VERY
C LARGE (E.G. 50.) THE RESULTING CURVE IS VERY NEARLY A
C POLYGONAL LINE. THE SIGN OF SIGMA INDICATES WHETHER
C SLOPE IN FORMATION HAS BEEN INPUT OR NOT. IF SIGMA IS
C NEGATIVE THE END-POINT SLOPES WILL BE DETERMINED
C INTERNALLY. A STANDARD VALUE FOR SIGMA IS APPROXIMATELY
C 1. IN ABSOLUTE VALUE.
C ON OUTPUT--
C N,X,Y,SLP1,SLPN, AND SIGMA ARE UNALTERED,
C XP AND YP CONTAIN INFORMATION ABOUT THE CURVATURE OF THE
C CURVE AT THE GIVEN NODES,
C AND
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE.
      INTEGER N
      REAL X(N), Y(N), XP(N), YP(N), TEMP(N), S, SIGMA
      DEGRAD = 3.1415926535897932/180.
      NM1 = N - 1
      NP1 = N + 1
      DELX1 = X(2) - X(1)
      DELY1 = Y(2) - Y(1)
      DELS1 = SQRT(DELX1*DELX1+DELY1*DELY1)
      DX1 = DELX1/DELS1
      DY1 = DELY1/DELS1
C DETERMINE SLOPES IF NECESSARY
      IF (SIGMA.LT.0.) GO TO 70
      SLPP1 = SLP1*DEGRAD
      SLPPN = SLPN*DEGRAD
C SET UP RIGHT HAND SIDES OF TRIDIAGONAL LINEAR SYSTEM FOR XP
C AND YP
   10 XP(1) = DX1 - COS(SLPP1)
      YP(1) = DY1 - SIN(SLPP1)
      TEMP(1) = DELS1
      S = DELS1
      IF (N.EQ.2) GO TO 30
      DO 20 I=2,NM1
        DELX2 = X(I+1) - X(I)
        DELY2 = Y(I+1) - Y(I)
        DELS2 = SQRT(DELX2*DELX2+DELY2*DELY2)
        DX2 = DELX2/DELS2
        DY2 = DELY2/DELS2
        XP(I) = DX2 - DX1
        YP(I) = DY2 - DY1
        TEMP(I) = DELS2
        DELX1 = DELX2
        DELY1 = DELY2
        DELS1 = DELS2
        DX1 = DX2
        DY1 = DY2
C ACCUMULATE POLYGONAL ARCLENGTH
        S = S + DELS1
   20 CONTINUE
   30 XP(N) = COS(SLPPN) - DX1
      YP(N) = SIN(SLPPN) - DY1
C DENORMALIZE TENSION FACTOR
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S
C PERFORM FORWARD ELIMINATION ON TRIDIAGONAL SYSTEM
      DELS = SIGMAP*TEMP(1)
      EXPS = EXP(DELS)
      SINHS = .5*(EXPS-1./EXPS)
      SINHIN = 1./(TEMP(1)*SINHS)
      DIAG1 = SINHIN*(DELS*.5*(EXPS+1./EXPS)-SINHS)
      DIAGIN = 1./DIAG1
      XP(1) = DIAGIN*XP(1)
      YP(1) = DIAGIN*YP(1)
      SPDIAG = SINHIN*(SINHS-DELS)
      TEMP(1) = DIAGIN*SPDIAG
      IF (N.EQ.2) GO TO 50
      DO 40 I=2,NM1
        DELS = SIGMAP*TEMP(I)
        EXPS = EXP(DELS)
        SINHS = .5*(EXPS-1./EXPS)
        SINHIN = 1./(TEMP(I)*SINHS)
        DIAG2 = SINHIN*(DELS*(.5*(EXPS+1./EXPS))-SINHS)
        DIAGIN = 1./(DIAG1+DIAG2-SPDIAG*TEMP(I-1))
        XP(I) = DIAGIN*(XP(I)-SPDIAG*XP(I-1))
        YP(I) = DIAGIN*(YP(I)-SPDIAG*YP(I-1))
        SPDIAG = SINHIN*(SINHS-DELS)
        TEMP(I) = DIAGIN*SPDIAG
        DIAG1 = DIAG2
   40 CONTINUE
   50 DIAGIN = 1./(DIAG1-SPDIAG*TEMP(NM1))
      XP(N) = DIAGIN*(XP(N)-SPDIAG*XP(NM1))
      YP(N) = DIAGIN*(YP(N)-SPDIAG*YP(NM1))
C PERFORM BACK SUBSTITUTION
      DO 60 I=2,N
        IBAK = NP1 - I
        XP(IBAK) = XP(IBAK) - TEMP(IBAK)*XP(IBAK+1)
        YP(IBAK) = YP(IBAK) - TEMP(IBAK)*YP(IBAK+1)
   60 CONTINUE
      RETURN
   70 IF (N.EQ.2) GO TO 80
C IF NO SLOPES ARE GIVEN, USE SECOND ORDER INTERPOLATION ON
C INPUT DATA FOR SLOPES AT ENDPOINTS
      DELS2 = SQRT((X(3)-X(2))**2+(Y(3)-Y(2))**2)
      DELS12 = DELS1 + DELS2
      C1 = -C(DELS12+DELS1)/DELS12/DELS1
      C2 = DELS12/DELS1/DELS2
      C3 = -DELS1/DELS12/DELS2
      SX = C1*X(1) + C2*X(2) + C3*X(3)
      SY = C1*Y(1) + C2*Y(2) + C3*Y(3)
      SLPP1 = ATAN2(SY,SX)
      DELNM1 = SQRT((X(N-2)-X(NM1))**2+(Y(N-2)-Y(NM1))**2)
      DELN = SQRT((X(NM1)-X(N))**2+(Y(NM1)-Y(N))**2)
      DELNN = DELNM1 + DELN
      C1 = (DELNN+DELN)/DELNN/DELN
      C2 = -DELNN/DELN/DELNM1
      C3 = DELN/DELNN/DELNM1
      SX = C3*X(N-2) + C2*X(NM1) + C1*X(N)
      SY = C3*Y(N-2) + C2*Y(NM1) + C1*Y(N)
      SLPPN = ATAN2(SY,SX)
      GO TO 10
C IF ONLY TWO POINTS AND NO SLOPES ARE GIVEN, USE STRAIGHT
C LINE SEGMENT FOR CURVE
   80 XP(1) = 0.
      XP(2) = 0.
      YP(1) = 0.
      YP(2) = 0.
      RETURN
      END
      SUBROUTINE KURV2(T, XS, YS, N, X, Y, XP, YP, S, SIGMA)               D  10
      INTEGER N
      REAL T, XS, YS, X(N), Y(N), XP(N), YP(N), S, SIGMA
C THIS SUBROUTINE PERFORMS THE MAPPING OF POINTS IN THE
C INTERVAL (0.,1.) ONTO A CURVE IN THE PLANE. THE SUBROUTINE
C KURV1 SHOULD BE CALLED EARLIER TO DETERMINE CERTAIN
C NECESSARY PARAMETERS. THE RESULTING CURVE HAS A PARAMETRIC
C REPRESENTATION BOTH OF WHOSE COMPONENTS ARE SPLINES UNDER
C TENSION AND FUNCTIONS OF THE POLYGONAL ARCLENGTH PARAMETER
C ON INPUT--
C T CONTAINS A REAL VALUE OF ABSOLUTE VALUE LESS THAN OR
C EQUAL TO 1. TO BE MAPPED TO A POINT ON THE CURVE. THE
C SIGN OF T IS IGNORED AND THE INTERVAL (0.,1.) IS MAPPED
C ONTO THE ENTIRE CURVE. IF T IS NEGATIVE THIS INDICATES
C THAT THE SUBROUTINE HAS BEEN CALLED PREVIOUSLY (WITH ALL
C OTHER INPUT VARIABLES UNALTERED) AND THAT THIS VALUE OF
C T EXCEEDS THE PREVIOUS VALUE IN ABSOLUTE VALUE. WITH
C SUCH INFORMATION THE SUBROUTINE IS ABLE TO MAP THE POINT
C MUCH MORE RAPIDLY. THUS IF THE USER SEEKS TO MAP A
C SEQUENCE OF POINTS ONTO THE SAME CURVE, EFFICIENCY IS
C GAINED BY ORDERING THE VALUES INCREASING IN MAGNITUDE
C AND SETTING THE SIGNS OF ALL BUT THE FIRST, NEGATIVE,
C N CONTAINS THE NUMBER OF POINTS WHICH WERE INTERPOLATED
C TO DETERMINE THE CURVE,
C X AND Y ARE ARRAYS CONTAINING THE X- AND Y-COORDINATES
C OF THE INTERPOLATED POINTS,
C XP AND YP ARE THE ARRAYS OUTPUT FROM KURV2 CONTAINING
C CURVATURE INFORMATION,
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE,
C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).
C THE PARAMETERS N,X,Y,XP,YP,S,AND SIGMA SHOULD BE INPUT
C UNALTERED FROM THE OUTPUT OF KURV1.
C ON OUTPUT--
C XS AND YS CONTAIN THE X-ANDY-COORDINATES OF THE IMAGE
C POINT ON THE CURVE.
C T,N,X,Y,XP,YP,S, AND SIGMA ARE UNALTERED.
C DENORMALIZE SIGMA
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S
C STRETCH UNIT INTERVAL INTO ARCLENGTH DISTANCE
      TN = ABS(T*S)
C FOR NEGATIVE T START SEARCH WHERE PREVIOUSLY TERMINATED,
C OTHERWISE START FROM BEGINNING
      IF (T.LT.0.) GO TO 10
      I1 = 2
      XS = X(1)
      YS = Y(1)
      SUM = 0.
      IF (T.LE.0.) RETURN
   10 CONTINUE
C DETERMINE INTO WHICH SEGMENT TN IS MAPPED
      DO 30 I=I1,N
        DELX = X(I) - X(I-1)
        DELY = Y(I) - Y(I-1)
        DELS = SQRT(DELX*DELX+DELY*DELY)
        IF (SUM+DELS-TN) 20, 40, 40
   20   SUM = SUM + DELS
   30 CONTINUE
C IF ABS(T) IS GREATER THAN 1., RETURN TERMINAL POINT ON
C CURVE
      XS = X(N)
      YS = Y(N)
      RETURN
C SET UP AND PERFORM INTERPOLATION
   40 DEL1 = TN - SUM
      DEL2 = DELS - DEL1
      EXPS1 = EXP(SIGMAP*DEL1)
      SINHD1 = .5*(EXPS1-1./EXPS1)
      EXPS = EXP(SIGMAP*DEL2)
      SINHD2 = .5*(EXPS-1./EXPS)
      EXPS = EXPS1*EXPS
      SINHS = .5*(EXPS-1./EXPS)
      XS = (XP(I)*SINHD1+XP(I-1)*SINHD2)/SINHS +
     * ((X(I)-XP(I))*DEL1+(X(I-1)-XP(I-1))*DEL2)/DELS
      YS = (YP(I)*SINHD1+YP(I-1)*SINHD2)/SINHS +
     * ((Y(I)-YP(I))*DEL1+(Y(I-1)-YP(I-1))*DEL2)/DELS
      I1 = I
      RETURN
      END
      SUBROUTINE KURVP1(N, X, Y, XP, YP, TEMP, S, SIGMA)                   E  10
      INTEGER N
      REAL X(N), Y(N), XP(N), YP(N), TEMP(1), S, SIGMA
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 ON INPUT--
C N IS THE NUMBER OF POINTS TO BE INTERPOLATED (N.GE.2),
C X IS AN ARRAY CONTAINING THE N X-COORDINATES OF THE
C POINTS,
C Y IS AN ARRAY CONTAINING THE N Y-COORDINATES OF THE
C POINTS,
C XP,YP ARE ARRAYS OF LENGTH AT LEAST N,
C TEMP IS AN ARRAY OF LENGTH AT LEAST 2*N WHICH IS USED
C FOR SCRATCH STORAGE,
C AND
C SIGMA CONTAINS THE TENSION FACTOR. THIS IS A NON-ZERO
C QUANTITY (WHOSE SIGN IS IGNORED) WHICH INDICATES THE
C CURVINESS DESIRED. IF ABS(SIGMA) IS VERY LARGE (E.G. 50.
C ) THE RESULTING CURVE IS VERY A POLYGON. A STANDARD
C VALUE FOR SIGMA IS APPROXIMATELY 1. IN ABSOLUTE VALUE.
C ON OUTPUT--
C N,X,Y, AND SIGMA ARE UNALTERED,
C XP AND YP CONTAIN INFORMATION ABOUT THE CURVATURE OF THE
C CURVE AT THE GIVEN NODES,
C AND
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE.
      NM1 = N - 1
      NP1 = N + 1
C SET UP RIGHT HAND SIDES OF TRIDIAGONAL (WITH CORNER
C ELEMENTS) LINEAR SYSTEM FOR XP AND YP
      DELX1 = X(2) - X(1)
      DELY1 = Y(2) - Y(1)
      DELS1 = SQRT(DELX1*DELX1+DELY1*DELY1)
      DX1 = DELX1/DELS1
      DY1 = DELY1/DELS1
      XP(1) = DX1
      YP(1) = DY1
      TEMP(1) = DELS1
      S = DELS1
      DO 10 I=2,N
        IP1 = I + 1
        IF (I.EQ.N) IP1 = 1
        DELX2 = X(IP1) - X(I)
        DELY2 = Y(IP1) - Y(I)
        DELS2 = SQRT(DELX2*DELX2+DELY2*DELY2)
        DX2 = DELX2/DELS2
        DY2 = DELY2/DELS2
        XP(I) = DX2 - DX1
        YP(I) = DY2 - DY1
        TEMP(I) = DELS2
        DELX1 = DELX2
        DELY1 = DELY2
        DELS1 = DELS2
        DX1 = DX2
        DY1 = DY2
C ACCUMULATE POLYGONAL ARCLENGTH
        S = S + DELS1
   10 CONTINUE
      XP(1) = XP(1) - DX1
      YP(1) = YP(1) - DY1
C DENORMALIZE TENSION FACTOR
      SIGMAP = ABS(SIGMA)*FLOAT(N)/S
C PERFORM FORWARD ELIMINATION ON TRIDIAGONAL SYSTEM
      DELS = SIGMAP*TEMP(N)
      EXPS = EXP(DELS)
      SINHS = .5*(EXPS-1./EXPS)
      SINHIN = 1./(TEMP(N)*SINHS)
      DIAG1 = SINHIN*(DELS*.5*(EXPS+1./EXPS)-SINHS)
      DIAGIN = 1./DIAG1
      SPDIG1 = SINHIN*(SINHS-DELS)
      SPDIAG = 0.
      DO 20 I=1,N
        DELS = SIGMAP*TEMP(I)
        EXPS = EXP(DELS)
        SINHS = .5*(EXPS-1./EXPS)
        SINHIN = 1./(TEMP(I)*SINHS)
        DIAG2 = SINHIN*(DELS*(.5*(EXPS+1./EXPS))-SINHS)
        IF (I.EQ.N) GO TO 30
        DIAGIN = 1./(DIAG1+DIAG2-SPDIAG*TEMP(I-1))
        XP(I) = DIAGIN*(XP(I)-SPDIAG*XP(I-1))
        YP(I) = DIAGIN*(YP(I)-SPDIAG*YP(I-1))
        TEMP(N+I) = -DIAGIN*TEMP(NM1+I)*SPDIAG
        IF (I.EQ.1) TEMP(NP1) = -DIAGIN*SPDIG1
        SPDIAG = SINHIN*(SINHS-DELS)
        TEMP(I) = DIAGIN*SPDIAG
        DIAG1 = DIAG2
   20 CONTINUE
   30 TEMP(NM1) = TEMP(N+NM1) - TEMP(NM1)
      IF (N.EQ.2) GO TO 50
C PERFORM FIRST STEP OF BACK SUBSTITUTION
      DO 40 I=3,N
        IBAK = NP1 - I
        XP(IBAK) = XP(IBAK) - TEMP(IBAK)*XP(IBAK+1)
        YP(IBAK) = YP(IBAK) - TEMP(IBAK)*YP(IBAK+1)
        TEMP(IBAK) = TEMP(N+IBAK) - TEMP(IBAK)*TEMP(IBAK+1)
   40 CONTINUE
   50 XP(N) =
     * (XP(N)-SPDIG1*XP(1)-SPDIAG*XP(NM1))/(DIAG1+DIAG2+SPDIG1*T
     * EMP(1)+SPDIAG*TEMP(NM1))
      YP(N) =
     * (YP(N)-SPDIG1*YP(1)-SPDIAG*YP(NM1))/(DIAG1+DIAG2+SPDIG1*T
     * EMP(1)+SPDIAG*TEMP(NM1))
C PERFORM SECOND STEP OF BACK SUBSTITUTION
      DO 60 I=1,NM1
        XP(I) = XP(I) + TEMP(I)*XP(N)
        YP(I) = YP(I) + TEMP(I)*YP(N)
   60 CONTINUE
      RETURN
      END
      SUBROUTINE KURVP2(T, XS, YS, N, X, Y, XP, YP, S, SIGMA)              F  10
      INTEGER N
      REAL T, XS, YS, X(N), Y(N), XP(N), YP(N), S, SIGMA
C THIS SUBROUTINE PERFORMS THE MAPPING OF POINTS IN THE
C INTERVAL (0.,1.) ONTO A CLOSED CURVE IN THE PLANE. THE
C SUBROUTINE KURVP1 SHOULD BE CALLED EARLIER TO DETERMINE
C CERTAIN NECESSARY PARAMETERS. THE RESULTING CURVE HAS A
C PARAMETRIC REPRESENTATIONBOTH OF WHOSE COMPONENTS ARE
C PERIODIC SPLINES UNDER TENSION AND FUNCTIONS OF THE POLY-
C GONAL ARCLENGTH PARAMETER.
C ON INPUT--
C T CONTAINS A REAL VALUE OF ABSOLUTE VALUE LESS THAN OR
C EQUAL TO 1. TO BE MAPPED TO A POINT ON THE CURVE. THE
C SIGN OF T IS IGNORED AND THE INTERVAL (0.,1.) IS MAPPED
C ONTO THE ENTIRE CLOSED CURVE. IF T IS NEGATIVE THIS
C INDICATES THAT THE SUBROUTINE HAS BEEN CALLED PREVIOUSLY
C (WITH ALL OTHER INPUT VARIABLES UNALTERED) AND THAT
C THIS VALUE OF T EXCEEDS THE PREVIOUS VALUE IN ABSOLUTE
C VALUE. WITH SUCH INFORMATION THE SUBROUTINE IS ABLE TO
C MAP THE POINT MUCH MORE RAPIDLY. THUS IF THE USER SEEKS
C TO MAP A SEQUENCE OF POINTS ONTO THE SAME CURVE,
C EFFICIENCY IS GAINED BY ORDERING THE VALUES INCREASING
C IN MAGNITUDE AND SETTING THE SIGNS OF ALL BUT THE FIRST,
C NEGATIVE,
C N CONTAINS THE NUMBER OF POINTS WHICH WERE INTERPOLATED
C TO DETERMINE THE CURVE,
C X AND Y ARE ARRAYS CONTAINING THE X- AND Y-COORDINATES
C OF THE INTERPOLATED POINTS,
C XP AND YP ARE THE ARRAYS OUTPUT FROM KURVP1 CONTAINING
C CURVATURE INFORMATION,
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE,
C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).
C THE PARAMETERS N,X,Y,XP,YP,S AND SIGMA SHOULD BE INPUT
C UNALTERED FROM THE OUTPUT OF KURVP1.
C ON OUTPUT--
C XS AND YS CONTAIN THE X- AND Y-COORDINATES OF THE IMAGE
C POINT ON THE CURVE.
C T,N,X,Y,XP,YP,S AND SIGMA ARE UNALTERED.
C DENORMALIZE SIGMA
      SIGMAP = ABS(SIGMA)*FLOAT(N)/S
C STRETCH UNIT INTERVAL INTO ARCLENGTH DISTANCE
      TN = ABS(T*S)
C FOR NEGATIVE T START SEARCH WHERE PREVIOUSLY TERMINATED,
C OTHERWISE START FROM BEGINNING
      IF (T.LT.0.) GO TO 10
      I1 = 2
      SUM = 0.
   10 IF (I1.EQ.1) GO TO 50
C DETERMINE INTO WHICH SEGMENT TN IS MAPPED
      DO 30 I=I1,N
        DELX = X(I) - X(I-1)
        DELY = Y(I) - Y(I-1)
        DELS = SQRT(DELX*DELX+DELY*DELY)
        IF (SUM+DELS-TN) 20, 40, 40
   20   SUM = SUM + DELS
   30 CONTINUE
      I = 1
      IM1 = N
      DELS = S - SUM
      GO TO 50
   40 IM1 = I - 1
C SET UP AND PERFORM INTERPOLATION
   50 DEL1 = TN - SUM
      DEL2 = DELS - DEL1
      EXPS1 = EXP(SIGMAP*DEL1)
      SINHD1 = .5*(EXPS1-1./EXPS1)
      EXPS = EXP(SIGMAP*DEL2)
      SINHD2 = .5*(EXPS-1./EXPS)
      EXPS = EXPS1*EXPS
      SINHS = .5*(EXPS-1./EXPS)
      XS = (XP(I)*SINHD1+XP(IM1)*SINHD2)/SINHS +
     * ((X(I)-XP(I))*DEL1+(X(IM1)-XP(IM1))*DEL2)/DELS
      YS = (YP(I)*SINHD1+YP(IM1)*SINHD2)/SINHS +
     * ((Y(I)-YP(I))*DEL1+(Y(IM1)-YP(IM1))*DEL2)/DELS
      I1 = I
      RETURN
      END