Subroutine LSQSPL(Xpt,Ypt,Nxpt,X,Y,Npt,Coefs,Itype) c c *********************************************************************** c * * c * Subroutine LSQSPL - To fit least squares cubic spline * c * segments through your data. The splines satisfy point and * c * slope continuity at the junction points, and use end * c * conditions of zero second derivative. The user provides * c * the X values of the control points at which the segments * c * meet. * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Dimension Xpt(1),Ypt(1),X(1),Y(1),Coefs(1) Dimension Matrix(5,9),Hold(60,4),H(5),H1(5),Holds(60,2) Real Hold Double Precision A,B,C,D,Matrix,Save,H1,H,Holds,V,Xsum,Xpoint, $ X2sum,X3sum,X4sum,Ysum,Xysum,Xnumbr,Xmult, $ X5sum,Xy2sum c Nloop = Nxpt - 1 c Do 60 Iloop = 1,Nloop c c Find the start, end of the segment relative to the c Data to fit c Istart = -1 Iend = Npt Do 10 I=1,Npt If (Xpt(Iloop) .LE. X(I) .AND. Istart .EQ. -1) Istart = I If (Xpt(Iloop+1) .LT. X(I)) Then Iend = I - 1 Go To 15 End If 10 Continue 15 Continue Istart = Max0(1,Istart) Ntot = Iend - Istart + 1 c c Now compute the sums of the squares, etc c Xnumbr = 0.0 Xsum = 0.0 Xysum = 0.0 Ysum = 0.0 X2sum = 0.0 X3sum = 0.0 X4sum = 0.0 X5sum = 0.0 Xy2sum = 0.0 c Do 20 I=1,Ntot Xpoint = X(I+Istart-1) Xmult = 1.0 Xnumbr = Xnumbr + Xmult Xmult = Xmult * Xpoint Xsum = Xsum + Xmult Xysum = Xysum + Xmult*Y(I+Istart-1) Ysum = Ysum + Y(I+Istart-1) Xmult = Xmult * Xpoint Xy2sum = Xy2sum + Xmult*Y(I+Istart-1) X2sum = X2sum + Xmult Xmult = Xmult * Xpoint X3sum = X3sum + Xmult Xmult = Xmult * Xpoint X4sum = X4sum + Xmult Xmult = Xmult * Xpoint X5sum = X5sum + Xmult 20 Continue c c Zero the "matrix" c Do 25 I=1,5 Do 25 J=1,9 Matrix(I,J) = 0.0 25 Continue c c Start to load the matrix c If (Iloop .EQ. 1) Then If (Itype .EQ. 1) Then c c Load in end condition 2nd derivative = 0 c Matrix(1,1) = 6.0 * Xpt(Iloop) Matrix(1,2) = 2.0 Else c c Load in another least squares equation c Matrix(1,1) = X5sum Matrix(1,2) = X4sum Matrix(1,3) = X3sum Matrix(1,4) = X2sum Matrix(1,9) = Xy2sum End If Else c c Load in the last line of the previous matrix solve c Matrix(1,1) = H(1) Matrix(1,2) = H(2) Matrix(1,3) = H(3) Matrix(1,4) = H(4) Matrix(1,9) = H(5) End If c c Load in the least squares conditions c Matrix(2,1) = X3sum Matrix(2,2) = X2sum Matrix(2,3) = Xsum Matrix(2,4) = Xnumbr Matrix(2,9) = Ysum Matrix(3,1) = X4sum Matrix(3,2) = X3sum Matrix(3,3) = X2sum Matrix(3,4) = Xsum Matrix(3,9) = Xysum c If (Iloop .EQ. Nloop) Then If (Itype .EQ. 1) Then c c End condition - set the second derivative to zero c Matrix(4,1) = 6.0 * Xpt(Iloop+1) Matrix(4,2) = 2.0 Else c c Load in another least squares equation c Matrix(4,1) = X5sum Matrix(4,2) = X4sum Matrix(4,3) = X3sum Matrix(4,4) = X2sum Matrix(4,9) = Xy2sum End If Matrix(5,5) = 1.0 Else c c Point continuity condition c K = Iloop + 1 Matrix(4,1) = Xpt(K)*Xpt(K)*Xpt(K) Matrix(4,2) = Xpt(K)*Xpt(K) Matrix(4,3) = Xpt(K) Matrix(4,4) = 1.0 Matrix(4,5) = -Matrix(4,1) Matrix(4,6) = -Matrix(4,2) Matrix(4,7) = -Xpt(K) Matrix(4,8) = -1.0 c c Slope continuity condition c Matrix(5,1) = 3.0 * Xpt(K) * Xpt(K) Matrix(5,2) = 2.0 * Xpt(K) Matrix(5,3) = 1.0 Matrix(5,5) = -Matrix(5,1) Matrix(5,6) = -Matrix(5,2) Matrix(5,7) = -1.0 End If c c Now start to solve the "matrix" c Do 50 I=1,4 c c Pivit the rows; However, do not pivit rows 4 and 5 c up into rows 1, 2, or 3. c Ilp = I Amax = 0.0 Ll = 3 If (I .GT. 3) Ll = 5 c c Find the row with the largest Matrix(*,I) term c Do 30 J=I,Ll If (Abs(Matrix(J,I)) .GT. Amax) Then Amax = Abs(Matrix(J,I)) Ilp = J End If 30 Continue c c Pivit the rows, and nondimensionalize row I c Aval = Matrix(Ilp,I) Do 35 J=I,9 V = Matrix(Ilp,J) Matrix(Ilp,J) = Matrix(I,J) Matrix(I,J) = V / Aval 35 Continue c c Now do a gaussian elimination c Do 45 J=I+1,5 Save = Matrix(J,I) Do 40 L=I,9 Matrix(J,L) = Matrix(J,L) - Save * Matrix(I,L) 40 Continue 45 Continue 50 Continue c c Right now we have c c 1 X X X X 0 0 0 0 | X ---- From previous loop c _______________________ c 1 X X X 0 0 0 0 | X c 0 1 X X 0 0 0 0 | X c 0 0 1 X 0 0 0 0 | X c 0 0 0 1 X X X X | X c 0 0 0 0 X X X X | X c Matrix(2,4) = Matrix(2,4) - Matrix(3,4) * Matrix(2,3) Matrix(2,9) = Matrix(2,9) - Matrix(3,9) * Matrix(2,3) Matrix(1,4) = Matrix(1,4) - Matrix(2,4) * Matrix(1,2) - $ Matrix(3,4) * Matrix(1,3) Matrix(1,9) = Matrix(1,9) - Matrix(2,9) * Matrix(1,2) - $ Matrix(3,9) * Matrix(1,3) c c Eliminate terms from previous matrix solve c If (Iloop .GT. 1) Then Holds(Iloop-1,1) = H1(4) - H1(2)*Matrix(2,4) - $ H1(1)*Matrix(1,4) - H1(3)*Matrix(3,4) Holds(Iloop-1,2) = H1(5) - H1(2)*Matrix(2,9) - $ H1(1)*Matrix(1,9) - Matrix(3,9)*H1(3) End If c c Now we have: c c c 1 0 0 0 X 0 0 0 0 | X ---- From previous loop c _______________________ c 1 0 0 X 0 0 0 0 | X | c 0 1 0 X 0 0 0 0 | X | Current matrix c 0 0 1 X 0 0 0 0 | X | c 0 0 0 1 X X X X | X | c 0 0 0 0 X X X X | X ---- Insert back into next matrix top c H1(1) = Matrix(4,5) H1(2) = Matrix(4,6) H1(3) = Matrix(4,7) H1(4) = Matrix(4,8) H1(5) = Matrix(4,9) c c Save the partially solved matrix here c Hold(Iloop,1) = Matrix(1,4) Hold(Iloop,2) = Matrix(1,9) Hold(Iloop,3) = Matrix(2,4) Hold(Iloop,4) = Matrix(2,9) Ypt(Iloop) = Matrix(3,4) ! Temporary storage Coefs(Iloop) = Matrix(3,9) ! Temporary storage c c c Save the remaining portion of the matrix in H - c Will fit it back in to the matrix top row. c H(1) = Matrix(5,5) H(2) = Matrix(5,6) H(3) = Matrix(5,7) H(4) = Matrix(5,8) H(5) = Matrix(5,9) c c Go back for the next segment c 60 Continue c c Now we will backfit to get the results we want c D = H1(5) Do 80 I = Nloop,1,-1 A = - D * Hold(I,1) + Hold(I,2) B = - D * Hold(I,3) + Hold(I,4) C = - D * Ypt(I) + Coefs(I) c c Save the results as Y values at the segment endpoints c and the 2nd derivative at the segment endpoints c Xv = Xpt(I+1) Ypt(I+1) = ((((A * Xv) + B) * Xv) + C) * Xv + D Coefs(I+1) = (3.0*A*Xv + 2.0*B)*Xv + C c If (I .GT. 1) Then D = - D * Holds(I-1,1) + Holds(I-1,2) End If 80 Continue c Xv = Xpt(1) Ypt(1) = ((((A * Xv) + B) * Xv) + C) * Xv + D Coefs(1) = (3.0*A*Xv + 2.0*B)*Xv + C c Return End Subroutine Lssgev(Xnew,Ynew,Nnew,X,Y,Npt,Coefs) c c *********************************************************************** c * * c * Subroutine Lssgev - To evaluate intermediate points * c * on a cubic spline defined by point and 1st derivative * c * continuity only. * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Dimension Xnew(1),Ynew(1),X(1),Y(1),Coefs(1) c Istart = 1 Iflag = 1 J = 1 c Do 40 I=1,Nnew If (Xnew(I) .LE. X(1)) Go To 20 c c Find the interval that the X value lies on c Do 10 J=Istart,Npt-1 If (Xnew(I) .GE. X(J) .AND. Xnew(I) .LE. X(J+1)) Go To 20 10 Continue J = Npt-1 20 Continue c c Compute the coefficents of the cubic equation c If (J .GT. Istart .OR. Iflag .EQ. 1) Then Iflag = 0 Call Lsccmp(X,Y,Coefs,J,A,B,C,D) End If c c Compute the new Y value c Ynew(I) = D + Xnew(I) * (C + Xnew(I) * (B + Xnew(I) * A)) c If (I .LT. Nnew) Then If (Xnew(I+1) .GE. Xnew(I)) Then Istart = J Else Istart = 1 End If End If c 40 Continue Return End Subroutine LSCINT(Xnew,Ynew,Newpts,Xvals,Yvals,Npoint,Coefs) c c *********************************************************************** c * * c * Subroutine LSCINT - To Integrate Least Squares 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 Lsccmp(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 Lsccmp(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 LSCDER(Xnew,Ynew,Newpts,Xvals,Yvals,Npoint,Coefs) c c *********************************************************************** c * * c * Subroutine LSCDER - To compute the derivative of the * c * least squares 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 If (I .LT. Newpts) Then If (Xnew(I+1) .GE. Xnew(I)) Then Inrval = J Else Inrval = 1 End If End If c c Compute the derivative value c Call Lsccmp(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 Lsccmp(X,Y,Coefs,J,A,B,C,D) c Dimension X(1),Y(1),Coefs(1) c E3 = X(J)*X(J)*X(J) - X(J+1)*X(J+1)*X(J+1) E2 = X(J)*X(J) - X(J+1)*X(J+1) E1 = X(J) - X(J+1) A11 = 3.0 * E2 A12 = 2.0 * E1 A21 = E3 - 3.0 * E1 * X(J)*X(J) A22 = E2 - 2.0 * E1 * X(J) A = ((Coefs(J) - Coefs(J+1))*A22 - ((Y(J)-Y(J+1)) - $ Coefs(J)*E1)*A12) / (A11*A22 - A21*A12) B = (Coefs(J) - Coefs(J+1) - 3.0 * A * E2) / (2.0 * E1) C = Coefs(J) - X(J) * (2.0 * B + 3.0 * A * X(J)) D = Y(J) - X(J) * (C + X(J) * (B + X(J) * A)) c Return End