Subroutine Csplot(Xvals,Yvals,Numbr,X,Y,Np) Parameter (Npt_Max=2001) Dimension Xvals(1),Yvals(1),X(1),Y(1),Coefs(Npt_Max) c Call Csplin(X,Y,Np,Coefs) Call Cseval(Xvals,Yvals,Numbr,X,Y,Np,Coefs) Return End Subroutine CSEVAL(Xnew,Ynew,Newpts,Xvals,Yvals,Npoint,Coefs) c c ************************************************************************** c * * c * Subroutine CSEVAL - To interpolate Y values for given X values * c * along a cubic spline. The cubic spline fit was performed by * c * Subroutine CSPLIN. This subroutine assumes that the X values * c * to interpolate by are given in increasing order. * c * * c * Input Values: * c * * c * Xnew - The X values at which to interpolate along the cubic * c * spline * c * Newpts - The number of points to interpolate * c * Xvals - The original X values to which the cubic spline was * c * fitted * c * Yvals - The original Y values to which the cubic spline was * c * fitted * c * Npoint - The number of points to which the cubic spline was * c * fitted * c * Coefs - The second derivative values, computed by CSPLIN, * c * of the cubic spline * c * * c * Output Values: * c * * c * Ynew - The interpolated values of Y, corresponding to the * c * X values in the array Xnew * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c ************************************************************************** c Dimension Xnew(1),Ynew(1),Xvals(1),Yvals(1),Coefs(1) c Inrval = 1 Do 30 I=1,Newpts c c Find the interval along which Xnew(I) resides (interpolation), c or, if Xnew(I) resides outside the values of Xvals (extrapolation), c the nearest interval c Do 10 J = Inrval,Npoint-1 If( Xnew(I) .LT. Xvals(1) .OR. (Xnew(I) .GE. Xvals(J) .AND. $ Xnew(I) .LE. Xvals(J+1))) Go To 20 10 Continue J = Npoint-1 20 Continue c c Increment interval counter so the next search will start in the c same interval c Inrval = J c c Compute the Y value c Dum1 = Xvals(J+1)-Xvals(J) Dum2 = Xvals(J+1)-Xnew(I) Dum3 = Xvals(J)-Xnew(I) Ynew(I) = (Coefs(J)*Dum2*(Dum2*Dum2/Dum1-Dum1)-Coefs(J+1)* $ Dum3*(Dum3*Dum3/Dum1-Dum1))/6.0+(Yvals(J)*Dum2- $ Yvals(J+1)*Dum3)/Dum1 30 Continue Return End Subroutine CSPLIN(Xvals,Yvals,Npoint,Coefs) c c ************************************************************************** c * * c * Subroutine CSPLIN - to fit data values with a cubic spline. * c * This routine assumes that x is nondecreasing, and that the * c * second derivative is zero at the endpoints of the curve. * c * * c * Input variables: * c * * c * Xvals - An array of X values * c * Yvals - An array of Y values corresponding to the X values * c * Npoint - The number of points in the X and Y arrays * c * * c * Output Variables: * c * * c * Coefs - An array containing the second derivative values of * c * the cubic splines * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c ************************************************************************** c Parameter(Npt_Max=2001) Dimension Xvals(1),Yvals(1),Coefs(1),D(Npt_Max) c c First, eliminate the values of X and Y that have deincreasing c X values c Npts = 1 Icount = 1 Do 10 I = 2,Npoint If (Xvals(I) .EQ. Xvals(Npts)) Then Xvals(Npts) = (Xvals(Npts)*Icount + Xvals(I))/(Icount + 1) Yvals(Npts) = (Yvals(Npts)*Icount + Yvals(I))/(Icount + 1) Icount = Icount + 1 Else If (Xvals(I) .GT. Xvals(Npts)) Then Npts = Npts + 1 Xvals(Npts) = Xvals(I) Yvals(Npts) = Yvals(I) Icount = 1 End If 10 Continue Npoint = Npts c c Set up a series of simultaneous equations, with the second c derivatives as the unknowns. This will create a Tri-Diagonal c Matrix. Since the Matrix is Tri-Diagonal, only the non-zero c elements will be saved, in one one dimensional arrays named c C, and the column vector named D. c If (Npoint .EQ. 3) Then c c Solution for Npoint = 3 c Tempvl = (Yvals(3)-Yvals(2))/(Xvals(3)-Xvals(2))/(Xvals(3)- $ Xvals(1)) Coefs(2) = 3.0*(Tempvl-(Yvals(2)-Yvals(1))/(Xvals(2)-Xvals(1))/ $ (Xvals(3)-Xvals(1))) c Else If (Npoint .GT. 3) Then c c For Npoint > 3 c 20 Delta2 = Xvals(2)-Xvals(1) Do 30 I=1,Npoint-2 Delta1 = Xvals(I+2)-Xvals(I+1) B = 2.0*(Delta2+Delta1)/Delta2 Coefs(I) = Delta1 / Delta2 D(I) = (6.0*((Yvals(I+2)-Yvals(I+1))/Delta1-(Yvals(I+1)- $ Yvals(I))/Delta2)/Delta2) c c Solve the Tri-Diagonal Matrix for the second derivative values c If (I .EQ. 1) Then c D(1) = D(1) / B Coefs(1) = Coefs(1) / B Else c B = B - Coefs(I-1) D(I) = (D(I) - D(I-1)) / B Coefs(I) = Coefs(I) / B End If Delta2 = Delta1 30 Continue c Coefs(Npoint-1) = D(Npoint-2) Do 50 I=1,Npoint-3 J = Npoint-2-I D(J) = D(J)-Coefs(J)*D(J+1) Coefs(J+1) = D(J) 50 Continue c End If c Coefs(1) = 0.0 Coefs(Npoint) = 0.0 Return End Subroutine PSPLOT(Xvals,Yvals,Numbr,X,Y,Npt,Itype) c c *********************************************************************** c * * c * Subroutine PSPLOT - To plot parametric spline curves * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Parameter(Npt_Max=2001) Dimension Xvals(1),Yvals(1),X(1),Y(1),Coef(Npt_Max) Dimension Tvals(Npt_Max),Tnew(300) c c Loop twice to get the X and Y values - important because c CSPLIN may remove some points using the pathlength based T values c and duplicate X,Y values c Npt1 = Npt c c Compute the T values c If (Itype .EQ. 1) Then c c Compute Tvals, where Tvals(I) = I c Do 10 I=1,Npt Tvals(I) = I 10 Continue c Else c c Compute Tvals, where Tvals(I) = path length c Tvals(1) = 0.0 Do 20 I=2,Npt Tvals(I) = Tvals(I-1)+SQRT((Xvals(I)-Xvals(I-1))**2+ $ (Yvals(I)-Yvals(I-1))**2) 20 Continue c End If c c Compute the T increments c Tinc = (Tvals(Npt1) - Tvals(1))/(Numbr-1) Do 30 I=1,Numbr Tnew(I) = (I-1)*Tinc + Tvals(1) 30 Continue c c Compute the X and Y values corresponding to the T values c Call Csplin(Tvals,Y,Npt1,Coef) Call Cseval(Tnew,Yvals,Numbr,Tvals,Y,Npt1,Coef) Npt1 = Npt Call Csplin(Tvals,X,Npt1,Coef) Call Cseval(Tnew,Xvals,Numbr,Tvals,X,Npt1,Coef) c Return End