Perform a general least squares fit with optional error estimates.

   This version uses SVD.  A user-supplied function or a built-in
   polynomial is fit to the data.

   Function fitting

   Result = SVDFIT(X, Y, M)

   X:   A vector representing the independent variable.  If this an array,
           the columns are taken to be the precomputed independant vectors
           and no actual function is computed here.

   Y:   Dependent variable vector.  This vector should be same length 
      as X.

   M:   The number of coefficients in the fitting function.  For 
      polynomials, M is equal to the degree of the polynomial + 1.

   Weight:   A vector of weights for Y(i).  This vector should be the same
      length as X and Y.

      If this parameter is ommitted, 1 is assumed.  The error for 
      each term is weighted by Weight(i) when computing the fit.  
      Frequently, Weight(i) = 1./Sigma(i) where Sigma is the 
      measurement error or standard deviation of Y(i).

   Funct:   A string that contains the name of an optional user-supplied 
      basis function with M coefficients. If omitted, polynomials
      are used.

      The function is called:
         R = FUNCT(X,M)
      where X is an N element vector, and the function value is an 
      (N, M) array of the N inputs, evaluated with the M basis 
      functions.  M is analogous to the degree of the polynomial +1 
      if the basis function is polynomials.  For example, see the 
      function COSINES, in the IDL User Library, which returns a 
      basis function of:
         R(i,j) = cos(j*x(i)).
      For more examples, see Numerical Recipes, page 519.

      The basis function for polynomials, is R(i,j) = x(i)^j.
   SVDFIT returns a vector of M coefficients.

   NOTE:  In order for an optional keyword output parameter
   to be returned, it must be defined before calling SVDFIT.
   The value or structure doesn't matter.  For example:

      YF = 1            ;Define output variable yf.
      C = SVDFIT(X, Y, M, YFIT = YF)    ;Do SVD, fitted Y vector is now
                  ;returned in variable YF.

   YFIT:   Vector of calculated Y's.

   CHISQ:   Sum of squared errors multiplied by weights if weights
      are specified.

   COVAR:   Covariance matrix of the coefficients.

    VARIANCE:   Sigma squared in estimate of each coeff(M).

    SINGULAR:   The number of singular values returned.  This value should
      be 0.  If not, the basis functions do not accurately
      characterize the data.



   Adapted from SVDFIT, from the book Numerical Recipes, Press,
   et. al., Page 518.
   minor error corrected April, 1992 (J.Murthy)
   93/10/12, Marc W. Buie, Lowell Observatory.  Added option to make this
             work similar to "regress".
   97/03/20, MWB, Changed to use SVDC and SVSOL and everything is now in
               double precision.
   2005/02/02, MWB, Changed thresh from 10^-9 to 2x10^-12
   2005/06/21, MWB, added error trap keyword
   2009/07/14, MWB, fixed a bug in the original covar calculation
   2013/09/14, AMZ, fixed a bug for using default polynomials.