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 If (I .NE. Newpts) Then If (Xnew(I+1) .GE. Xnew(I)) Then Inrval = J Else Inrval = 1 End If End If 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 monotonically increasing, and * c * that the second derivative is zero at the endpoints of the * c * 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(Nptmax=2001) Dimension Xvals(1),Yvals(1),Coefs(1),D(Nptmax) 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 three one dimensional arrays named c A, B, and 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 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 PSPLIN(Xvals,Yvals,Npoint,Tvals,Coefx,Coefy,Itype) c c ***************************************************************************** c * * c * Subroutine PSPLIN - To fit parametric cubic splines to a * c * series of X and Y points. X and Y are related as: * c * X = a*T*T + b*T + c ; Y = d*T*T + e*T + f . * c * There are two ways to define T: one way assumes that at Xvals(1), * c * Yvals(1), T=1, and, at Xvals(I), Yvals(I), T=I. This is selected * c * by setting Itype = 1. The other way to define T is to relate it * c * to the path length; at Xvals(1),Yvals(1), T=0, and, at Xvals(N), * c * Yvals(n), T(N)=Sum of path length from Xvals(1),Yvals(1). This * c * is selected by setting Itype = 2. * c * * c * Input Values: * c * * c * Xvals - The X values to fit the Parametric cubic splines to. * c * Yvals - The Y values to fit the Parametric cubic splines to. * c * Npoint - The number of X and Y values to fit. * c * Itype - The type of definition for Tvals - * c * = 1 for Tvals(N) = N * c * = 2 for Tvals(N) = Path length * c * * c * Output values: * c * * c * Tvals - The values of T at Xvals,Yvals * c * Coefx - The values of the T,X second derivatives at * c * Tvals(I),Xvals(I) * c * Coefy - The Values of the T,Y second derivatives at * c * Tvals(I),Yvals(I) * c * * c ***************************************************************************** c Dimension Xvals(1),Yvals(1),Tvals(1),Coefx(1),Coefy(1) c If ( Itype .EQ. 2 ) Go To 20 c c Compute Tvals, where Tvals(I) = I c Do 10 I=1,Npoint Tvals(I) = I 10 Continue Go To 40 c c Compute Tvals, where Tvals(I) = path length c 20 Tvals(1) = 0.0 Do 30 I=2,Npoint Tvals(I) = Tvals(I-1)+SQRT((Xvals(I)-Xvals(I-1))**2+(Yvals(I)- $ Yvals(I-1))**2) 30 Continue c c Now fit the splines to the data: c 40 Continue Call CSPLIN(Tvals,Xvals,Npoint,Coefx) Call CSPLIN(Tvals,Yvals,Npoint,Coefy) Return End Subroutine PSEVAL(Xnew,Ynew,Newpts,Xvals,Yvals,Npoint,Tvals, $ Coefx,Coefy) c c *************************************************************************** c * * c * Subroutine PSEVAL - To compute Y values from known X values * c * along a curve defined using parametric splines. The * c * parametric spline data must come from subroutine PSPLIN. It * c * is assumed that the known X values are arranged in * c * increasing path length. * c * * c * Input values: * c * * c * Xnew - The known X values to interpolate the Y values at * c * Newpts - The number of Y values to find the corresponding X * c * values * c * Xvals - The original X values to which the parametric spline * c * was fitted * c * Yvals - the original Y values to which the parametric spline * c * was fitted * c * Npoint - The number of original points * c * Tvals - The values of T at the original points * c * Coefx - The values of the T,X second derivatives at the * c * original points * c * Coefy - The values of the T,Y second derivatives at the * c * original points * c * * c * NOTE: Tvals, Coefx, and Coefy are computed for the user by PSPLIN * c * * c * Output values: * c * * c * Ynew - The Y values corresponding to the X values in "Xnew" * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *************************************************************************** c Parameter(Nptmax=2001) Dimension Xvals(1),Yvals(1),Tvals(1),Ynew(1),Xnew(1),Coefx(1), $Coefy(1),Tnew(Nptmax) c Npoinx = Npoint c Do 10 I=1,Npoint Ynew(I) = Tvals(I) 10 Continue c c Compute the T values corresponding to the input X values c Call HIGHPT(Tvals,Xvals,Coefx,Npoinx) Call CSYTOX(Tnew,Xnew,Newpts,Tvals,Xvals,Npoinx,Coefx) c c Compute the Y values corresponding to the T values computed above c Do 20 I=1,Npoint Tvals(I) = Ynew(I) 20 Continue Call CSEVAL(Tnew,Ynew,Newpts,Tvals,Yvals,Npoint,Coefy) Return End Subroutine CSYTOX(Xnew,Ynew,Newpts,Xvals,Yvals,Npoint,Coefs) c c *************************************************************************** c * * c * Subroutine CSYTOX - To compute the X value for a given Y value * c * along a cubic spline. It is assumed that the cubic spline was * c * created using subroutine CSPLINE, and that the Y values are * c * increasing. * c * * c * Input values: * c * * c * Ynew - The Y values at which to solve for X at * c * Newpts - The number of points to solve for * c * Xvals - The X coordinates to which the spline was fitted * c * Yvals - The Y coordinates to which the spline was fitted * c * Npoint - The number of points to which the spline was fitted * c * Coefs - The values of the second derivative of the spline * c * at the points to which the spline was fitted * c * * c * Output values: * c * * c * Xnew - The X values corresponding to the Y values in "Ynew" * 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 50 I=1,Newpts c c Find the interval where Ynew(I) exists c Do 30 J=Inrval,Npoint-1 c c Set flags for interpolation and extrapolation: c If (( Ynew(I) .LT. Yvals(1) .AND. Inrval .EQ. 1 )) Then Iflag = 1 Go To 10 Else If (Ynew(I) .GT. Yvals(Npoint) .AND. Inrval .GE. $ (Npoint-1)) Then Iflag = 2 Go To 10 Else Iflag = 0 End If c c For Interpolation: c If ( Ynew(I) .EQ. Yvals(J+1) .AND. J .LT. $ Npoint-1 ) Go To 30 Sign = (Ynew(I)-Yvals(J))*(Ynew(I)-Yvals(J+1)) If ( Sign .GT. 0.0 ) Go To 30 c c When desired interval has been found, come here c 10 Continue Icount = 0 X1 = Xvals(J) X2 = Xvals(J+1) Dx = X2 - X1 c c Set the initial interval c Xa = Xvals(J) Xb = Xvals(J+1) Ya = Yvals(J) Yb = Yvals(J+1) c c Compute the cubic spline coefficients c A = (Coefs(J+1) - Coefs(J))/(6.0 * Dx) B = (Coefs(J+1)/2.0 - 3.0*A*X2) C = (Yvals(J+1) - Yvals(J) - (X2*X2 - X1*X1)*B - (X2*X2* $ X2 - X1*X1*X1)*A)/Dx D = Yvals(J+1) - X2*C - X2*X2*B - X2*X2*X2*A c 20 Continue c c Compute the X value to use c Slope = (Yb - Ya)/(Xb - Xa) If (Slope .EQ. 0.0) Slope = 1.0E-6 X = (Ynew(I)-Ya)/Slope+Xa c c Ensure that X doesn't fall out of our interval c If (Iflag .EQ. 0) Then If (X .LT. Xa .OR. X .GT. Xb) X = (Xa + Xb)/2.0 Else If (Iflag .EQ. 1) Then If (X .GT. Xa) X = Xa Else If (X .LT. Xb) X = Xb End If c c Compute new y value c Y = ((A*X + B)*X + C)*X + D c c Check to see if the result is close enough; if so, then get c the next value of Ynew. Otherwise, repeat with new X guess. c But first, see if we have exceeded the number of iterations allowed. c If (Icount .GT. 15) Go To 40 Icount = Icount + 1 c Error = ABS(Y-Ynew(I))-ABS(Yvals(J+1)-Yvals(J))/1000.0 If ( Error .GT. 1E-10 ) Then c c Set the new interval to use: c If (Y .GT. Ynew(I) .AND. Yb .GE. Ya .OR. Y .LT. Ynew(I) $ .AND. Yb .LE. Ya) Then Yb = Y Xb = X Else Xa = X Ya = Y End if If (Xb .EQ. Xa .OR. Yb .EQ. Ya) Go To 40 c Go To 20 End If c Go To 40 30 Continue J = Npoint - 1 c c Save the last X guess. Then, set the interval counter to the c current interval, so the search for the next value will start c in this interval c 40 Xnew(I) = X c If (I .NE. Newpts) Then If (Xnew(I+1) .GE. Xnew(I)) Then Inrval = J Else Inrval = 1 End If End If c 50 Continue Return End Subroutine HIGHPT(Xvals,Yvals,Coefs,Npoint) c c *********************************************************************** c * * c * Subroutine HIGHPT - To find the maximum and minimum values * c * along a curve fitted with a spline, and then "add" these * c * points to the original points AND the second derivatives. * c * This is useful when trying to compute X values from known * c * Y values. * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Parameter(Nptmax = 2001) Dimension Xvals(1),Yvals(1),Coefs(1) c c Shift the data in Xvals, Yvals, and Coefs to the end of the c arrays. This way we can load intermediate points into c the arrays. (Note: This routine assumes that the arrays are c 200 elements in size.) c Maxpt = Nptmax Do 5 I=1,Npoint Xvals(Maxpt + 1 - I) = Xvals(Npoint+1-I) Yvals(Maxpt + 1 - I) = Yvals(Npoint+1-I) Coefs(Maxpt + 1 - I) = Coefs(Npoint+1-I) 5 Continue c Ii = 0 Do 10 I=1,Npoint-1 c c Save the data for the Ith point c Ii = Ii + 1 Iset = Maxpt - Npoint + I Xvals(Ii) = Xvals(Iset) Yvals(Ii) = Yvals(Iset) Coefs(Ii) = Coefs(Iset) c c Compute some of the cubic spline coefficients c X1 = Xvals(Iset) X2 = Xvals(Iset+1) Dx = X2 - X1 A = (Coefs(Iset+1) - Coefs(Iset))/(6.0 * Dx) B = (Coefs(Iset+1)/2.0 - 3.0*A*X2) C = (Yvals(Iset+1) - Yvals(Iset) - (X2*X2 - X1*X1)*B - (X2*X2* $ X2 - X1*X1*X1)*A)/Dx c c See if a maximum or minimum exists c Aval = 4.0*B*B - 12.0*A*C If (Aval .GE. 0.0) Then c D = Yvals(Iset+1) - X2*C - X2*X2*B - X2*X2*X2*A c c Compute the maximum and minimum X locations c If (A .NE. 0.0) Then Xa = (-2.0*B - Sqrt(Aval))/(6.0*A) Xb = (-2.0*B + Sqrt(Aval))/(6.0*A) If (Xb .EQ. Xa) Xb = 1E20 Else If (B .NE. 0.0) Then Xa = -C/(2.0*B) Xb = 1E20 Else Xa = 0.0 Xb = 1E20 End If c c See if the X values lie on this interval c If ( X1 .LE. Xa .AND. Xa .LE. X2) Then Ii = Ii + 1 Xvals(Ii) = Xa Yvals(Ii) = A*Xa*Xa*Xa + B*Xa*Xa + C*Xa + D Coefs(Ii) = 6.0*A*Xa + 2.0*B End If If ( X1 .LE. Xb .AND. Xb .LE. X2) Then Ii = Ii + 1 Xvals(Ii) = Xb Yvals(Ii) = A*Xb*Xb*Xb + B*Xb*Xb + C*Xb + D Coefs(Ii) = 6.0*A*Xb + 2.0*B End If End If 10 Continue Ii = Ii + 1 Xvals(Ii) = Xvals(Maxpt) Yvals(Ii) = Yvals(Maxpt) Coefs(Ii) = Coefs(Maxpt) c c Save the new values back into the original arrays c Npoint = Ii c Return End