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 + 1
 1    Continue
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 1E15.  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.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
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
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) 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