Use Chebyshev polynomials in fits

* fit_utils.py: By switching to Chebyshev polynomials T_k(y) instead of
  regular ones, we make it much easier to put bounds on f(y) = Sum(a_k
  T_k(y)).  This is because each T_k(y) is bounded by [-1, 1] so long as
  we remap the domain of y to [-1, 1].
Restrict possible values of component widths

* fit_utils.py: Do not allow the sigmas to fall below 2.0 or above 10.0,
  and also curtail the linear and quadratic terms to amplitude of 1.0
More Chebyshev-ization

* fit-nebula.py:  Adapt the nebular fits to use the new Chebyshev
  method.
Fix bug with proliferating HDUs in FITS files

* fit-nebula.py (main): Try to write the new images to an existing HDU
  of the same name and only create a new HDU if that fails.
Allow manual specification of nebular components

* fit-nebula.py:  Add many new command line parameters to allow greater
  control over the nebular components.  There may now be 1, 2, or 3, and
  the first 2 only get initialised to the 2-delta solution if they are
  not set manually.