; ----------------------------------------------------------------------
;   GJ_polyfit...  (Polynomial fitting with Gauss-Jordan elimination)...
;
;   Input: npt      : Number of points (npt must be > order+1).
;          x[]      : x values.
;          y[]      : y values.
;          wt[]     : weight values (usually 1/sig^2).
;                       (Entries with wt=0 are excluded from the data set.)
;          order    : Polynomial order (0 to 99).
;          forcexoff: Force routine to use the given "xoff".
;                       0 = no force, xoff determined by this routine.
;                       1 = force, use user xoff value for offset.
;  In/Out: xoff     : X value offset. Subtract this value from x values
;                        before using coeffecients.
;  Output: ucoef[]  : User polynomial coeffecients (from fit # >= order+1).
;
;    Returns "1" if successful, otherwise returns "0".
;
;> use poly_fit
;> take some care to avoid 1/sqrt(negative)

function GJ_polyfit, npt, x, y, wt, order, forcexoff, xoff, ucoef

  xx = x[0:npt-1]
  yy = y[0:npt-1]
  ww = wt[0:npt-1]
  indx = where(ww gt 0., count)
  if count eq 0 then return, 0
  xx = xx[indx]
  yy = yy[indx]
  measure_errors = 1/sqrt(ww[indx])
  if forcexoff eq 1 then xx = xx - xoff else xoff = 0.
  ucoef = POLY_FIT( xx, yy, order, /DOUBLE, MEASURE_ERRORS=measure_errors)

  return, 1

end