Fitting an Optical Image with No PSF

Here we will fit a elliptical Gaussian + constant background model to the DDS image of NGC 3079. After starting XIMGFIT,

ximgfit> data tutorial/ngc3079.fits.gz 100 100 400 400
ximgfit> model cons gauss
ximgfit> setpar 1 2600
ximgfit> setpar 4 12000.
ximgfit> setpar 5 15
ximgfit> setpar 6 96
ximgfit> setpar 7 -13.5
ximgfit> fit

``data'' loads the data from the fits image, extracting the range (100,100) - (400,400). ``model'' sets up the model with default values, and the ``setpar'' commands give better initial guesses for some of the parameters. This should produce output similar to:

7057201.14 -4 2674 2.359 1.072 1e+04 15.63 95.2 -13.5
No new minimim found after 1 tries
7057201.14 -3 2674 2.359 1.072 1e+04 15.63 95.2 -13.5
7056707.04 -4 2674 2.12 1.176 1e+04 15.67 95.5 -13.51
...
6829262.12 0 2675 2.672 1.36 1e+04 16.32 101.7 -13.52
Last change in fit statistic = 977.0837
ximgfit>

Here XIMGFIT is improving the fit, showing the value of $\chi^2$, the fitting matrix ``scale factor'' (ignore this), and the value of the free fit parameters. The ``No new minimim found after 1 tries'' etc., warnings are refering to the fact that the current iteration of perturbing fit values did not reduce $\chi^2$, and accordingly a different scaling parameter (the log of which is shown after $\chi^2$ in each iteration) must be tried (see Numerical Recipes for more information). Since the last change in $\chi^2$ was large, run fit several more times, this time with a large number of iterations by doing ``fit 100''. This may be rather time consuming because the image size is 300 $X$ 300, and each fit iteration requires the evaluation of the fitting function and its partial derivative(s) for each pixel (in this case, $7 \times 300 \times 300 = 630000$ calculations per iteration). After getting a good fit, ``show all'' should produce something like:

Data file: tutorial/ngc3079.fits.gz (100, 100) - (400, 400), size = 301 X 301
PSF: none
Mask file: none (100, 100) - (400, 400), size = 301 X 301
Padding =  none
Rebinning =  none
Plate scale = 1.700000 arcsec/pixel
Fitting statistic: chisq
Model:  Constant Gaussian 
Par Mod   Name            Value (     units)        Description
  1   1  const         2.67e+03 (   unknown)        Constant intensity
  2   2     xc             2.92 (    arcsec)        X offset from image center
  3   2     yc            0.996 (    arcsec)        Y offset from image center
  4   2   norm         1.17e+04 (      comp)        Normalization
  5   2 sigmax             15.3 (    arcsec)        Sigma in the X direction
  6   2 sigmay             96.1 (    arcsec)        Sigma in the Y direction
  7   2  theta            -13.6 (   degrees)        Rotation angle
chisq = 6564727.84, with 90594 degrees of freedom (90601 data bins)

If your terminal can display X-windows, then you should see a ``Tk'' window and can display the data, model and residuals using the ``view'' command:

ximgfit> view data
ximgfit> view model
ximgfit> view resid

which should produce the images shown in Figures 1.

Figure: Result of ``view data'' (top left), ``view model'' (top right), and ``view resid'' (bottom) within XIMGFIT.

Finally, save the current fit using save all n3079_dss_fit.ximg or whatever filename you want. To restore the fit, start XIMGFIT with ximgfit n3079_dss_fit.ximg or, once running XIMGFIT, do @n3079_dss_fit.ximg.

Of course, in this case a Gaussian is only a crude model of the optical continuum, although subtracting this model does a decent job of bringing out the dust lanes. A proper fit of optical images of spiral galaxies should use models such as a deVauculeurs law for the disk and a $e^{-r^{-1/4}}$ for any bulge component.