Subroutine LSQSEG(Xnew,Ynew,Nnew,X,Y,Npts,Iorder) c c *********************************************************************** c * * c * Subroutine LSQSEG - To fit points with least squares * c * fits in small segments. * c * * c * Programmed by Steven C. Fisher, DTNSRDC Code 1522 * c * * c *********************************************************************** c Dimension Xnew(1),Ynew(1),X(1),Y(1),Coefs(21) Double Precision Xpoint c c Decode the order of the polynomial, span, and c weight function c Iord = Iorder / 1000 Iorder = Iorder - Iord * 1000 Nspan = Iorder/10 Iw = Iorder - Nspan * 10 c If (Iw .EQ. 0) Iweigt = 1 If (Iw .EQ. 1) Iweigt = 2 If (Iw .EQ. 2) Iweigt = 5 If (Iw .EQ. 3) Iweigt = 25 If (Iw .EQ. 4) Iweigt = 100 c Jpoint = 1 Do 30 I=1,Nnew If (Xnew(I) .LE. X(1)) Then Istart = 1 Ncentr = 1 Else If (Xnew(I) .GE. X(Npts)) Then Istart = Npts Ncentr = Npts Else Do 10 J=Jpoint,Npts If (Xnew(I) .LE. X(J)) Go To 20 10 Continue 20 Continue Jpoint = J Istart = J Ncentr = J End If c c Load up the temporary arrays with data c Istart = Max0(Istart-Nspan/2,1) Istart = Min0(Istart+Nspan-1,Npts) - Nspan + 1 Ncentr = Ncentr - Istart + 1 Call Lssgft(X(Istart),Y(Istart),Nspan,Iord,Coefs,Iweigt, $ Ncentr) Yv = Coefs(1) Xpoint = 1.D0 Do 25 J=2,Iord+1 Xpoint = Xpoint * Xnew(I) Yv = Yv + Xpoint * Coefs(J) 25 Continue Ynew(I) = Yv 30 Continue Return End Subroutine LSSGFT(Xvals,Yvals,Npoint,Ndegre,Coefs,Iweigt,Ncentr) c c ************************************************************************ c * * c * Subroutine LSSGFT - To fit a least squares polynomial to * c * a set of data, with a weight function on the points * 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 * Iweigt - The weight to use on the points * c * Ncentr - The "center" of the weight function * 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(18,19),Power(45) Real*8 Matrix,Power,Xpower,Corr,Xpoint,Wfunct,Slope,Ak,F Integer Ycolum c c Set up the matrix c Ndeg = 1 + Ndegre If (Iweigt .GT. 0) Then F = 1.0D0/Iweigt Else F = -Iweigt End If c 1 Continue c If (Npoint .GT. 2) Then Ak = Max0(Npoint/2,Ncentr-1) Slope = (1.0D0-F) / Ak End If c Ndeg = Ndeg - 1 Nterm = Ndeg + 1 Loop = 2*Nterm Ycolum = Ndeg + 2 Coefs(Nterm+1) = 0.0 c c See if the largest elements will overflow c Amax = 0.0 Do 3 I=1,Npoint Amax = Amax1(Amax,Abs(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. 15.0) Then Corr = 10.0**( 15./Loop - Alog10(Amax)) Else If (Loop*Alog10(Amax) .LT. -15.0) Then Corr = 10.0**(-15./Loop - Alog10(Amax)) End If c 5 Continue c c Initialize arrays used for summing c Do 10 I=1,Loop Power(I) = 0.0D0 10 Continue Do 20 I=1,Nterm Matrix(I,Ycolum) = 0.0D0 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 If (J .EQ. Ncentr) Then Ak = Max0(Npoint-Ncentr,(Npoint-1)/2) Slope = (F-1.0D0) / Ak End If c Wfunct = Slope * (J - Ncentr) + 1.0D0 c Xpower = 1.0D0 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 * Wfunct 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*Wfunct 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,18) c c Load the solution into Coefs c Xpoint = 1.0D0 Do 80 I=1,Nterm If (Dlog10(Xpoint+0.1D0) + Dlog10(Dabs(Matrix(I,Ycolum))+ $ 0.1D0) .GT. 35.0) Go To 1 c Coefs(I) = Xpoint * Matrix(I,Ycolum) c If (Dlog10(Xpoint+0.1D0)+Dlog10(Corr+0.1D0) .GT. 1.0D35) $ Go To 1 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) Go To 1 c Return End