Subroutine Cocomp(Xvals,Yvals,Coefs,I,A,B,C,D) c c *********************************************************************** c * * c * Subroutine COCOMP - To compute the spline coefficients * c * A * (X^3) + B * (X^2) + C * (X) + D * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Dimension Xvals(1),Yvals(1),Coefs(1) c X1 = Xvals(I) X2 = Xvals(I+1) Dx = X2 - X1 A = (Coefs(I+1) - Coefs(I))/(6.0 * Dx) B = (Coefs(I+1)/2.0 - 3.0*A*X2) C = (Yvals(I+1) - Yvals(I) - (X2*X2 - X1*X1)*B - (X2*X2* $ X2 - X1*X1*X1)*A)/Dx D = Yvals(I+1) - X2*C - X2*X2*B - X2*X2*X2*A c Return End Subroutine CSINT(Xnew,Ynew,Newpts,Xvals,Yvals,Npoint,Coefs) c c *********************************************************************** c * * c * Subroutine CSINT - To Integrate Cubic Splines * 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 Arsum = 0.0 Iflag = 0 c c Start Loop over the number of points to evaluate: c Do 30 I=1,Newpts 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 c If (Xnew(1) .LE. Xvals(J+1)) Then If (Iflag .EQ. 1) Go To 10 c c Compute the area over an interval with no "new" points c Call Cocomp(Xvals,Yvals,Coefs,J,A,B,C,D) Y2 = (((A/4.0 * Xvals(J+1) + B/3.0) * Xvals(J+1) + $ C/2.0) * Xvals(J+1) + D) * Xvals(J+1) Y1 = (((A/4.0 * Xvals(J) + B/3.0) * Xvals(J) + $ C/2.0) * Xvals(J) + D) * Xvals(J) Arsum = Arsum + (Y2 - Y1) If (J .GE. Npoint-1) Iflag = 1 End If 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 Coefficients to the spline segment c Call Cocomp(Xvals,Yvals,Coefs,J,A,B,C,D) c c Set the integration limits c If (Xnew(I) .LT. Xvals(1) .OR. Xnew(1) .GT. Xvals(J)) Then X1 = Xnew(1) Arsum = 0.0 Else If (Xnew(I) .GT. Xvals(J+1)) Then X1 = Xvals(J+1) Else X1 = Xvals(J) End If X2 = Xnew(I) c c Integrate the spline from its initial point to the new X location c Ynew(I) = (((A/4.0 * X2 + B/3.0) * X2 + $ C/2.0) * X2 + D) * X2 - $ (((A/4.0 * X1 + B/3.0) * X1 + $ C/2.0) * X1 + D) * X1 + Arsum c If (Xnew(I) .LT. Xvals(1) .OR. Xnew(1) .GT. Xvals(J)) Then Arsum = (((A/4.0 * Xvals(J) + B/3.0) * Xvals(J) + $ C/2.0) * Xvals(J) + D) * Xvals(J) - $ (((A/4.0 * Xnew(1) + B/3.0) * Xnew(J) + $ C/2.0) * Xnew(1) + D) * Xnew(1) End If c 30 Continue Return End Subroutine PSINT(Xnew,Ynew,Tnew,Newpts,Xvals,Yvals,Tvals,Npoint, $ Coefx,Coefy) c c *********************************************************************** c * * c * Subroutine PSINT - To Integrate Parametric Splines * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Dimension Xvals(1),Yvals(1),Tvals(1),Ynew(1),Xnew(1),Coefx(1), $Coefy(1),Tnew(1) c Inrval = 1 Arsum = 0.0 Iflag = 0 Tneww = Tnew(1) c c Start Loop over the number of points to evaluate: c Do 30 I=1,Newpts Do 10 J = Inrval,Npoint-1 If( Tnew(I) .LT. Tvals(1) .OR. (Tnew(I) .GE. Tvals(J) .AND. $ Tnew(I) .LE. Tvals(J+1))) Go To 20 c If (Tneww .LE. Tvals(J+1)) Then If (Iflag .EQ. 1) Go To 10 c c Compute the area over an interval with no "new" points c Call Cocomp(Tvals,Xvals,Coefx,J,Ax,Bx,Cx,Dx) Call Cocomp(Tvals,Yvals,Coefy,J,Ay,By,Cy,Dy) Call Pscint(Y1,Tvals(J),Ax,Bx,Cx,Dx,Ay,By,Cy,Dy) Call Pscint(Y2,Tvals(J+1),Ax,Bx,Cx,Dx,Ay,By,Cy,Dy) Arsum = Arsum + (Y2 - Y1) If (J .EQ. Npoint-1) Iflag = 1 End If 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 Coefficients to the spline segment c Call Cocomp(Tvals,Xvals,Coefx,J,Ax,Bx,Cx,Dx) Call Cocomp(Tvals,Yvals,Coefy,J,Ay,By,Cy,Dy) c c Set the integration limits c If (Tnew(I) .LT. Tvals(1) .OR. Tneww .GT. Tvals(J)) Then T1 = Tneww Arsum = 0.0 Else If (Tnew(I) .GT. Tvals(J+1)) Then T1 = Tvals(J+1) Else T1 = Tvals(J) End If c c Integrate the spline from its initial point to the new X location c Call Pscint(Y1,T1,Ax,Bx,Cx,Dx,Ay,By,Cy,Dy) Call Pscint(Y2,Tnew(I),Ax,Bx,Cx,Dx,Ay,By,Cy,Dy) Y = Y2 - Y1 + Arsum c If (Tnew(I) .LT. Tvals(1) .OR. Tneww .GT. Tvals(J)) Then Call Pscint(Y1,Tnew(1),Ax,Bx,Cx,Dx,Ay,By,Cy,Dy) Call Pscint(Y2,Tvals(J),Ax,Bx,Cx,Dx,Ay,By,Cy,Dy) Arsum = Y2 - Y1 End If Ynew(I) = Y c 30 Continue Return End Subroutine Pscint(Yval,Tval,Ax,Bx,Cx,Dx,Ay,By,Cy,Dy) c c *********************************************************************** c * * c * Subroutine Pscint - To compute the definite integral of * c * parametric splines at a point (Tval) * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Yval = (((((Ax*Ay/2.0*Tval + (2.0*Bx*Ay+3.0*Ax*By)/5.0) * $ Tval + (Ay*Cx+2.0*Bx*By+3.0*Ax*Cy)/4.0) * Tval + $ (Cx*By+2.0*Bx*Cy+3.0*Ax*Dy)/3.0) * Tval + $ (Cx*Cy + 2.0*Bx*Dy)/2.0) * Tval + Cx*Dy) * Tval Return End Subroutine Slintg(Xnew,Ynew,Newpts,Xvals,Yvals,Npoint) c c *********************************************************************** c * * c * Subroutine Slintg - To Integrate a series of straight * c * line segments. * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Dimension Xnew(1),Ynew(1),Xvals(1),Yvals(1) c Inrval = 1 Arsum = 0.0 Iflag = 0 c c Start Loop over the number of points to evaluate: c Do 30 I=1,Newpts 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 c If (Xnew(1) .LE. Xvals(J+1)) Then If (Iflag .EQ. 1) Go To 10 c c Compute the area over an interval with no "new" points c Arsum = Arsum + (Yvals(J) + Yvals(J+1))/2.0 * $ (Xvals(J+1) - Xvals(J)) If (J .EQ. Npoint-1) Iflag = 1 End If 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 Set the integration limits c If (Xnew(I) .LT. Xvals(1) .OR. Xnew(1) .GT. Xvals(J)) Then X1 = Xnew(1) Arsum = 0.0 Else If (Xnew(I) .GT. Xvals(J+1)) Then X1 = Xvals(J+1) Else X1 = Xvals(J) End If X2 = Xnew(I) c c Integrate the segment from its initial point to the new X location c A = (Yvals(J+1)-Yvals(J))/(Xvals(J+1)-Xvals(J)) Y1 = A * X1 - A * Xvals(J) + Yvals(J) Y2 = A * X2 - A * Xvals(J) + Yvals(J) Ynew(I) = (Y2 + Y1)/2.0 * (X2 - X1) + Arsum c If (Xnew(I) .LT. Xvals(1) .OR. Xnew(1) .GT. Xvals(J)) Then Y1 = A * Xnew(1) - A * Xvals(J) + Yvals(J) Y2 = A * Xvals(J) - A * Xvals(J) + Yvals(J) Arsum = (Y2 + Y1)/2.0 * (Xvals(J) - Xnew(1)) End If c 30 Continue Return End Subroutine CSDER(Xnew,Ynew,Newpts,Xvals,Yvals,Npoint,Coefs) c c *********************************************************************** c * * c * Subroutine CSDER - To compute the derivative of the * c * cubic splines * 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 c Inrval = J c c Compute the derivative value c Call Cocomp(Xvals,Yvals,Coefs,J,A,B,C,D) Ynew(I) = (3.0*A * Xnew(I) + 2.0*B) * Xnew(I) + C c 30 Continue Return End Subroutine PSDER(Xnew,Ynew,Tnew,Newpts,Xvals,Yvals,Tvals,Npoint, $ Coefx,Coefy) c c *********************************************************************** c * * c * Subroutine PSDER - to compute the derivative of Parametric * c * Splines * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Dimension Xnew(1),Ynew(1),Xvals(1),Yvals(1),Coefx(1),Coefy(1), $ Tnew(1),Tvals(1) c Inrval = 1 Do 30 I=1,Newpts c c Find the interval along which Tnew(I) resides (interpolation), c or, if Tnew(I) resides outside the values of Tvals (extrapolation), c the nearest interval c Do 10 J = Inrval,Npoint-1 If( Tnew(I) .LT. Tvals(1) .OR. (Tnew(I) .GE. Tvals(J) .AND. $ Tnew(I) .LE. Tvals(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 c Inrval = J c c Compute the Y value c Call Cocomp(Tvals,Xvals,Coefx,J,Ax,Bx,Cx,Dx) Call Cocomp(Tvals,Yvals,Coefy,J,Ay,By,Cy,Dy) T = Tnew(I) Ynew(I) = ((3.0*Ay * T + 2.0*By) * T + Cy) / $ ((3.0*Ax * T + 2.0*Bx) * T + Cx) 30 Continue Return End Subroutine Slder(Xnew,Ynew,Newpts,Xvals,Yvals,Npoint) c c *********************************************************************** c * * c * Subroutine Slder - To compute the derivative along straight * c * line segments. * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Dimension Xnew(1),Ynew(1),Xvals(1),Yvals(1) c Inrval = 1 c c Start Loop over the number of points to evaluate: c Do 30 I=1,Newpts Do 10 J = Inrval,Npoint-1 If( Xnew(I) .LT. Xvals(1) .OR. (Xnew(I) .GE. Xvals(J) .AND. $ Xnew(I) .LT. Xvals(J+1))) Go To 20 c 10 Continue J = Npoint-1 20 Continue c c Increment interval counter so the next search will start in the c same interval c c Inrval = J c c Compute the derivative c If (Xnew(I) .NE. Xvals(J) .OR. J .EQ. 1) Then c c Compute derivative in the middle of a line: c Ynew(I) = (Yvals(J+1) - Yvals(J))/(Xvals(J+1) - Xvals(J)) Else c c Compute derivative where two lines join: c Ynew(I) = (Yvals(J+1) - Yvals(J-1))/(Xvals(J+1) - $ Xvals(J-1)) End If c 30 Continue Return End