Subroutine Lsqplt(Xvals,Yvals,Numbr,X,Y,Np,Iorder) Dimension Xvals(1),Yvals(1),X(1),Y(1),Coefs(30) c Call Lsqfit(X,Y,Np,Iorder,Coefs) Call Lseval(Xvals,Yvals,Numbr,Iorder,Coefs) Return End Subroutine LSQFIT(Xvals,Yvals,Npoint,Ndegre,Coefs) c c ************************************************************************ c * * c * Subroutine LSQFIT - To fit a least squares polynomial to * c * a set of data. * c * * c * Input values: * c * * c * Xvals - The array containing the X values to fit * c * Yvals - The array containing the Y values to fit * c * Npoint - The number of points to fit * c * Ndegre - The degree of the polynomial to fit * c * * c * Output Values: * c * * c * Coefs - The Coefficients of the polynomial, with Coefs(1) * c * being the constant term, and Coefs(Ndegre), the * c * highest term * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c ************************************************************************ c Dimension Xvals(1),Yvals(1),Coefs(1) Dimension Matrix(21,22),Power(45) Real*8 Matrix,Power,Xpower,Corr,Xpoint Integer Ycolum c c Set up the matrix c Ndeg = Ndegre c 1 Continue c Nterm = Ndeg + 1 Loop = 2*Nterm Ycolum = Ndeg + 2 c c See if the largest elements will overflow c Amax = 0.0 Do 3 I=1,Npoint Amax = Amax1(Amax,Xvals(I)) 3 Continue c Corr = 1.0 If (Amax .EQ. 0.0) Go To 5 c c Assume the largest X**Loop element we want is on c the order of 1E20. If it is larger, get a scale factor. c If (Loop*Alog10(Amax) .GT. 20.0) Then Corr = 10.0**( 20./Loop - Alog10(Amax)) Else If (Loop*Alog10(Amax) .LT. -20.0) Then Corr = 10.0**(-20./Loop - Alog10(Amax)) End If c 5 Continue c c Initialize arrays used for summing c Do 10 I=1,Loop Power(I) = 0.0 10 Continue Do 20 I=1,Nterm Matrix(I,Ycolum) = 0.0 20 Continue c c Compute here the coefficients to the matrix; all of the c coefficients are of the form: sum of X**(I-1) c Do 50 J=1,Npoint Xpower = 1.0 Xpoint = Xvals(J) * Corr Do 40 I=1,Loop c c Fill Power(I) with the sum of X**(I-1) c Power(I) = Power(I)+Xpower c c Load in the right hand side values into the matrix c (which are the sum of Y*X**(I-1) c If (I .LE. Nterm) Matrix(I,Ycolum) = $ Matrix(I,Ycolum)+Yvals(J)*Xpower Xpower = Xpower * Xpoint 40 Continue 50 Continue c c Now load the sums of X**(I-1) into the matrix c Do 70 I=1,Nterm Do 60 J=1,Nterm Matrix(I,J) = Power(I+J-1) 60 Continue 70 Continue c c Call to the matrix solution routine c Call Matinv(Matrix,Nterm,21) c c Load the solution into Coefs c Xpoint = 1.0 Do 80 I=1,Nterm If (Dlog10(Xpoint+0.1D0) + Dlog10(Dabs(Matrix(I,Ycolum))+ $ 0.1D0) .GT. 35.0) Then Ndeg = Ndeg - 1 Go To 1 End If c Coefs(I) = Xpoint * Matrix(I,Ycolum) c If (Dlog10(Xpoint+0.1D0)+Dlog10(Corr+0.1D0) .GT. 1.0D35) Then Ndeg = Ndeg - 1 Go To 1 End If c Xpoint = Xpoint * Corr c 80 Continue c c If the coefficient for the highest power of X is zero c (or almost), decrement the "power" by one, and recalculate c the least squares fit. c If (Abs(Coefs(Nterm)) .LT. 1E-35) Then Ndeg = Ndeg - 1 Coefs(Nterm) = 0.0 Go To 1 End If c Return End Subroutine LSEVAL(Xvals,Yvals,Npoint,Ndegre,Coefs) c c ************************************************************************ c * * c * Subroutine LSEVAL - To evaluate a polynomial at various X * c * locations. I realize that this subroutine may look * c * inefficient, but it is written this way to eliminate overflow * c * for large values of X at high powers. * c * * c * Input values: * c * * c * Xvals - The X values to evaluate the polynomial at * c * Npoint - The number of points to evaluate * c * Ndegre - The degree of the polynomial (cubic = 3) * c * Coefs - The polynomial coefficients, where Coefs(N) * c * is the coefficent for X**(N-1) * c * * c * Output values: * c * * c * Yvals - The computed Y values corresponding to the X * c * values in Xvals * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c ************************************************************************ c Dimension Xvals(1),Yvals(1),Coefs(1) Double Precision Xpoint c Nterm = Ndegre+1 Do 20 I=1,Npoint Yvals(I) = Coefs(1) Do 10 J=2,Nterm Xpoint = Coefs(J) Do 5 K=2,J Xpoint = Xpoint * Xvals(I) 5 Continue Yvals(I) = Yvals(I) + Xpoint 10 Continue 20 Continue Return End Subroutine MATINV(Matrix,Nsize,Ndim) c c *********************************************************************** c * * c * Subroutine Matinv - To solve a matrix using a Gauss * c * elimination method. The right hand sides of the * c * simultaneous equations are stored in Matrix(I,Nsize+1) * c * The solution to the matrix are then placed there. It * c * should be noted that the rest of the matrix is the * c * identity matrix. * c * * c * Input data: * c * * c * Matrix - The matrix containing the coefficients of the * c * simultaneous equations, with the right hand * c * side of the equation stored in Matrix(I,Nsize+1) * c * Nsize - The size of the matrix * c * Ndim - The first dimension of the Matrix - MUST match * c * that used in the dimension statement of the * c * calling program * c * * c * Output data: * c * * c * Matrix - The solution to the simultaneous equations are * c * stored in Matrix(I,Nsize+1) * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c * * c *********************************************************************** c Dimension Matrix(Ndim,1) Double Precision Matrix,Amax,Save Ncols = Nsize+1 Nloop = Nsize-1 Do 50 I=1,Nloop c c Find the row with the largest I th element magnitude. Then c swap its row with the ith row. c Irow = I Amax = 0.0 Do 10 J=I,Nsize If (Amax .GT. Dabs(Matrix(I,J))) Go To 10 Irow = J Amax = Dabs(Matrix(I,J)) 10 Continue c c Swap the rows, and normalize the new I th row by its element c along the main diagonal (which is the leftmost nonzero c element in that row) c Amax = Matrix(Irow,I) c If (Amax .EQ. 0.0) Amax = 1E-20 c If (Dabs(Amax) .LT. 1D-20) Amax = Dsign(1D-20,Amax) c Do 20 J=I,Ncols Save = Matrix(Irow,J) Matrix(Irow,J) = Matrix(I,J) Matrix(I,J) = Save/Amax 20 Continue c c Now, subtract the Ith row from the other rows c Iplus1 = I+1 Do 40 J=Iplus1,Nsize Save = Matrix(J,I) Do 30 K=I,Ncols Matrix(J,K) = Matrix(J,K)-Save*Matrix(I,K) 30 Continue 40 Continue 50 Continue c c The matrix now consists of 1's in the diagonal, and zeros below c the diagonal. Backfit in the answers in to solve the matrix. c If (Matrix(Nsize,Nsize) .EQ. 0.0) Matrix(Nsize,Nsize) = 1D-20 Matrix(Nsize,Ncols) = Matrix(Nsize,Ncols)/Matrix(Nsize,Nsize) Matrix(Nsize,Nsize) = 1.0 c Do 70 I=2,Nsize J = Nsize-I+1 Sum = Matrix(J,Ncols) Do 60 k = 2,I Nplace = Ncols-K+1 Sum = Sum-Matrix(J,Nplace)*Matrix(Nplace,Ncols) 60 Continue Matrix(J,Ncols) = Sum 70 Continue Return End