	FUNCTION curv1(N,X,Y,SLP1,SLPN,YP,TEMP,SIGMA,XF,YF,NF)
	IMPLICIT REAL*8  (A-H,O-Z)
C--	generates spline under tension for monotonic x
C--	NOTE THAT THIS IS A FUNCTION IN ORDER TO RETURN SUCCESS OR FAILURE
	REAL*8 X(N),Y(N),TEMP(N),YP(N),XF(NF),YF(NF)
	INTEGER*4 curv1,n,nf
	IF (SIGMA.EQ.0.) THEN
C	TYPE *,'tension cannot be 0'
	CURV1=0
	RETURN
	ENDIF
	NM1 = N-1
	NP1 = N+1
	DELX1 = X(2) - X(1)
	DX1 = (Y(2) - Y(1))/DELX1
	
	IF (SIGMA.LT.0) GO TO 50
	SLPP1 = SLP1
	SLPPN = SLPN
 
10    SIGMAP = ABS(SIGMA) * FLOAT(N-1)/(X(N) - X(1))
C	FOR INTERACTIVE USE, WE SLOW THIS ROUTINE DOWN BY DOING SOME
C	EXTRA CHECKING FOR THE USER, THIS COULD PREVENT MUCH TIME LOST
C	FROM A CRASH
	DX2=0.0
	DO 12	I=2,N
	DELX2=X(I)-X(I-1)
	IF (DELX2.LE.0) GOTO 990	!ERROR, EXIT
	DELS=SIGMAP*DELX2
	IF (DELS.GT.DX2) DX2=DELS	!FIND MAX
12	CONTINUE
	IF (DX2.GT.50.) THEN
C-- IF IT IS TOO BIG, WE GOT TROUBLE; WE WILL AUTO REDUCE TO FIT
C	TYPE	*,' warning, tension reduced to avoid overflows !'
	SIGMAP=50.*SIGMAP/DX2
	ENDIF
 
	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))
 
	DO 40 I=2,N
	   IBAK = NP1 - I
	   YP(IBAK) = YP(IBAK) - TEMP(IBAK)*YP(IBAK+1)
40    CONTINUE
	GOTO 200
50    IF (N.EQ.2) GO TO 60
 
	DELX2 = X(3) - X(2)
	DELX12 = X(3) - X(1)
	C1 = -(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
 
60    YP(1) = 0.
	YP(2) = 0.
200	CONTINUE
	I1=2
	DO 300	J=1,NF
	T=XF(J)
210    DO 220 I=I1,N
	   IF (X(I)-T) 220,220,230
220    CONTINUE
	I = N
 
230    IF (X(I-1).LE.T .OR. T.LE.X(1)) GO TO 240
 
	I1 = 2
	GO TO 210
 
240    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)
	YF(J) = (YP(I)*SINHD1+YP(I-1)*SINHD2)/SINHS +
     & ((Y(I)-YP(I))*DEL1+(Y(I-1)-YP(I-1))*DEL2)/DELS
	I1 = I
300	CONTINUE
	CURV1=1
	RETURN
990	CONTINUE		!ERROR
C	TYPE *,'independent vector not monotonically increasing'
	CURV1=0
	RETURN
	END
C===========================================================================
	FUNCTION kurv1(N,X,Y,SLP1,SLPN,XP,YP,TEMP,SIGMA,T,XS,YS,NF)
C--	generates spline under tension for any open (x,y) curve
	IMPLICIT REAL*8  (A-H,O-Z)
C--	NOTE THAT THIS IS A FUNCTION IN ORDER TO RETURN SUCCESS OR FAILURE

	REAL*8 X(N),Y(N),XP(N),YP(N),TEMP(N),S,SIGMA,XS(NF),YS(NF),T(NF)
	INTEGER*4 kurv1,n,nf
	DATA DEGRAD/.017453292/
	KURV1=1
	IF (SIGMA.EQ.0.) THEN
C	TYPE *,'tension cannot be 0'
	KURV1=0
	RETURN
	ENDIF
	NM1 = N - 1
	DELX1 = X(2) - X(1)
	DELY1 = Y(2) - Y(1)
	DELS1 = SQRT(DELX1*DELX1+DELY1*DELY1)
	DX1 = DELX1/DELS1
	DY1 = DELY1/DELS1
	S=DELS1
C--	check if slopes were input (yes if sigma>0)
	IF (SIGMA.GT.0) THEN
C--	convert input slopes from degress to radians
	SLPP1 = SLP1*DEGRAD
	SLPPN = SLPN*DEGRAD
	ELSE	!determine slopes
	IF (N.EQ.2) THEN
C--	if only 2 points and no slopes, we have to assume a straight line
	XP(1) = 0.
	XP(2) = 0.
	YP(1) = 0.
	YP(2) = 0.
	SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S
	ELSE
	DELS2 = SQRT((X(3)-X(2))**2+(Y(3)-Y(2))**2)
	DELS12 = DELS1 + DELS2
	C1 = -(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)
	XQ=SLPP1/DEGRAD
	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)
	ENDIF
	ENDIF
C--	slopes are now done
C	TYPE *,'slopes are',SLPP1,SLPPN
	XP(1) = DX1 - COS(SLPP1)
	YP(1) = DY1 - SIN(SLPP1)
	TEMP(1) = DELS1
	IF (N.GT.2) THEN
	DO  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
	   S = S + DELS1
	ENDDO
	ENDIF
	XP(N) = COS(SLPPN) - DX1
	YP(N) = SIN(SLPPN) - DY1
 
	SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S
 
	DELS = SIGMAP*TEMP(1)
	EXPS = EXP(DELS)
C	TYPE *,'DELS,EXPS =',DELS,EXPS
	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.GT.2) THEN
	DO  I=2,NM1
	   DELS = SIGMAP*TEMP(I)
	   EXPS = EXP(DELS)
C	TYPE *,'DELS,EXPS =',DELS,EXPS
	   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
	ENDDO
	ENDIF

	DIAGIN = 1./(DIAG1-SPDIAG*TEMP(NM1))
	XP(N) = DIAGIN*(XP(N)-SPDIAG*XP(NM1))
	YP(N) = DIAGIN*(YP(N)-SPDIAG*YP(NM1))
	
	DO  I=N-1,1,-1
	   XP(I) = XP(I) - TEMP(I)*XP(I+1)
	   YP(I) = YP(I) - TEMP(I)*YP(I+1)
	ENDDO
C--	done with the setup, now do the interpolations
C
	I=2
	SUM=0.0
	   DELX = X(2) - X(1)
	   DELY = Y(2) - Y(1)
	   DELS = SQRT(DELX*DELX+DELY*DELY)
	DO K=1,NF
C--	loop over the nf input values of T
	TN = ABS(T(K)*S)
10    CONTINUE
	IF (TN.LT.SUM) THEN
C--	drop back
	I=I-1
	IF (I.GE.2) THEN
	   DELX = X(I) - X(I-1)
	   DELY = Y(I) - Y(I-1)
	   DELS = SQRT(DELX*DELX+DELY*DELY)
	SUM=SUM-DELS
	GOTO 10
	ELSE
	XS(K)=X(1)
	YS(K)=Y(1)
	GOTO 50
	ENDIF
	ENDIF
	IF (TN.GT.(SUM+DELS)) THEN
C--	go forward
	I1=I
	DO I=I1,N
	   DELX = X(I) - X(I-1)
	   DELY = Y(I) - Y(I-1)
	   DELS = SQRT(DELX*DELX+DELY*DELY)
	   IF ( (SUM+DELS).GT.TN) GOTO 40
	SUM = SUM + DELS
	ENDDO 
	XS(K) = X(N)
	YS(K) = Y(N)
	GOTO 50
40	CONTINUE
	ENDIF
	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(K) = (XP(I)*SINHD1+XP(I-1)*SINHD2)/SINHS +
     & ((X(I)-XP(I))*DEL1+(X(I-1)-XP(I-1))*DEL2)/DELS
	YS(K) = (YP(I)*SINHD1+YP(I-1)*SINHD2)/SINHS +
     & ((Y(I)-YP(I))*DEL1+(Y(I-1)-YP(I-1))*DEL2)/DELS
50	CONTINUE
	enddo
	RETURN
	END
C=============================================================================
	FUNCTION kurvp1(N,X,Y,XP,YP,TEMP,SIGMA,T,XS,YS,NF)
C--	generates spline under tension for any closed loop (x,y) curve
	IMPLICIT REAL*8  (A-H,O-Z)
C--	NOTE THAT THIS IS A FUNCTION IN ORDER TO RETURN SUCCESS OR FAILURE
	INTEGER	N,kurvp1,NF
	REAL*8 X(N),Y(N),XP(N),YP(N),TEMP(1),S,SIGMA,XS(NF),YS(NF),T(NF)
	KURVP1=1
	IF (SIGMA.EQ.0.) THEN
C	TYPE *,'tension cannot be 0'
	KURVP1=0
	RETURN
	ENDIF
	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
	XP(1) = DX1
	YP(1) = DY1
	TEMP(1) = DELS1
	S = DELS1
	DO 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
 
	   S = S + DELS1
	ENDDO
	XP(1) = XP(1) - DX1
	YP(1) = YP(1) - DY1
 
	SIGMAP = ABS(SIGMA)*FLOAT(N)/S
 
	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 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))
C	   TEMP(N+1) = -DIAGIN*TEMP(NM1+I)*SPDIAG !was a bug, n+1 should be n+i
	   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
	ENDDO
30	CONTINUE
	TEMP(NM1) = TEMP(N+NM1) - TEMP(NM1)
	IF (N.EQ.2) GO TO 50
 
	DO 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)
	ENDDO
50    XP(N) =
     & (XP(N)-SPDIG1*XP(1)-SPDIAG*XP(NM1))/(DIAG1+DIAG2+SPDIG1*TEMP(1)
     & +SPDIAG*TEMP(NM1))
	YP(N) =
     & (YP(N)-SPDIG1*YP(1)-SPDIAG*YP(NM1))/(DIAG1+DIAG2+SPDIG1*TEMP(1)
     & +SPDIAG*TEMP(NM1))
 
	DO I=1,NM1
	   XP(I) = XP(I) + TEMP(I)*XP(N)
	   YP(I) = YP(I) + TEMP(I)*YP(N)
	ENDDO
C--	done with setup, now do the interpolations
C
	I=2
	SUM=0.0
	   DELX = X(2) - X(1)
	   DELY = Y(2) - Y(1)
	   DELS = SQRT(DELX*DELX+DELY*DELY)
	DO K=1,NF
C--	loop over the nf input values of T
	TN = ABS(T(K)*S)
	IF (TN.GT.S) TN=S
C	type *,'start, k,tn,sum,dels',k,tn,sum,dels
10    CONTINUE
	IF (TN.LT.SUM) THEN
C--	drop back
C	type *,'drop back, I=',I
	I=I-1
	IF (I.GE.2) THEN
	   DELX = X(I) - X(I-1)
	   DELY = Y(I) - Y(I-1)
	   DELS = SQRT(DELX*DELX+DELY*DELY)
	SUM=SUM-DELS
	GOTO 10
	ELSE
C--	I was 1, we must have dels since we never start with 1
	I=N
	SUM=SUM-DELS
C	type *,'drop back, I was 1'
	GOTO 10
	ENDIF
	ENDIF
	IF (TN.GT.(SUM+DELS)) THEN
C--	go forward
C	type *,'go forward, I=',I
	I1=I
	DO I=I1,N
	   DELX = X(I) - X(I-1)
	   DELY = Y(I) - Y(I-1)
	   DELS = SQRT(DELX*DELX+DELY*DELY)
	   IF ( (SUM+DELS).GT.TN) GOTO 40
	SUM = SUM + DELS
	ENDDO 
C--	must be between N and 1
	I=1
	DELS=S-SUM
	ENDIF
40	CONTINUE
	IM1=I-1
	if (im1.eq.0) im1=N
C	type *,'resultants, i,im1,tn,sum,dels',i,im1,tn,sum,dels
C	type *,'xp,yp',XP(I),YP(I)
	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(K) = (XP(I)*SINHD1+XP(IM1)*SINHD2)/SINHS +
     & ((X(I)-XP(I))*DEL1+(X(IM1)-XP(IM1))*DEL2)/DELS
	YS(K) = (YP(I)*SINHD1+YP(IM1)*SINHD2)/SINHS +
     & ((Y(I)-YP(I))*DEL1+(Y(IM1)-YP(IM1))*DEL2)/DELS
	enddo
	RETURN
	END
