NAME:
   mysvdfit

 PURPOSE:
   Perform a general least squares fit with optional error estimates.

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

 CATEGORY:
   Function fitting

 CALLING SEQUENCE:
   Result = SVDFIT(X, Y, M)

 INPUTS:
   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.

 OPTIONAL INPUTS:
   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.
      
 OUTPUTS:
   SVDFIT returns a vector of M coefficients.

 OPTIONAL OUTPUT PARAMETERS:
   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.

 COMMON BLOCKS:
   None.

 SIDE EFFECTS:
   None.

 MODIFICATION HISTORY:
   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.