; JMFIT
;---------------------------------------------------------------
;! fits gaussians to portions of an image
;# TASK ANALYSIS MODELING
;-----------------------------------------------------------------------
;;  Copyright (C) 1995-1999, 2001, 2003-2004
;;  Associated Universities, Inc. Washington DC, USA.
;;
;;  This program is free software; you can redistribute it and/or
;;  modify it under the terms of the GNU General Public License as
;;  published by the Free Software Foundation; either version 2 of
;;  the License, or (at your option) any later version.
;;
;;  This program is distributed in the hope that it will be useful,
;;  but WITHOUT ANY WARRANTY; without even the implied warranty of
;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;  GNU General Public License for more details.
;;
;;  You should have received a copy of the GNU General Public
;;  License along with this program; if not, write to the Free
;;  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,
;;  MA 02139, USA.
;;
;;  Correspondence concerning AIPS should be addressed as follows:
;;         Internet email: aipsmail@nrao.edu.
;;         Postal address: AIPS Project Office
;;                         National Radio Astronomy Observatory
;;                         520 Edgemont Road
;;                         Charlottesville, VA 22903-2475 USA
;-----------------------------------------------------------------------
JMFIT     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
JMFIT: Task to fit Gaussian models to an image by least-squares
USERID     -32000.0     32000.0    User ID. 0 => current user,
                                     32000 => any user.
INNAME                             Input image name (name)
INCLASS                            Input image name (class)
INSEQ           0.0      9999.0    Input image name (seq. #)
INDISK          0.0         9.0    Input image disk drive #
BLC             0.0      4096.0    Bottom Left corner of fit
TRC             0.0      4096.0    Top Right corner of fit
OUTNAME                            Output image name (name)
OUTCLASS                           Output image name (class)
OUTSEQ         -1.0      9999.0    Output image name (seq. #)
OUTDISK         0.0         9.0    Output image disk drive #
NGAUSS          0.0         4.0    Number of components
CTYPE           0.0         5.0    Model types; one for each
                                    component
                                    0->1=Gaussian
                                    2=Zero level; 3=Zero+slope
                                    4=Zero+slope+curve;
                                    5=see HELP JMFIT
                                   Guess of model parameters
GMAX                               Peak of component (JY)
                                     0-> Use maximum value
GPOS          -50.0     16434.0    (X,Y) position (pixels)
                                     0-> Use position of max
GWIDTH       -180.0       180.0    (BMAJ, BMIN, PA) of comp.
                                     (pixels,pixels,deg)
                                     0->Use clean beam
         @                         Fit model parameters
FMAX     @                         Peak of component (JY)
         @                           0-> Use maximum value
FPOS     @    -50.0     16434.0    (X,Y) position (pixels)
         @                           0-> Use position of max
FWIDTH   @   -180.0       180.0    (BMAJ, BMIN, PA) of comp.
         @                           (pixels,pixels,deg)
         @                           0->Use clean beam
DOMAX    $                         Solve for GMAX?  >0 -> yes
         $                           returns the error
DOPOS    $                         Solve for GPOS?  >0 -> yes
         $                           returns the error
DOWIDTH  $                         Solve for GWIDTH? >0-> yes
         $                           returns the error
BWSMEAR         0.0         0.1    Bandwidth smearing corr.
RADIUS                             Radius for finding RMS
NITER           0.0      4000.0    Maximum # of iterations
                                     0->20*NGAUSS
                                   Solve for model parameters?
DOCRT          -1.0       132.0    <=0 -> Print maps and
                                     solutions on line printer
OUTPRINT
                                   Printer disk file to save
DOOUTPUT       -1.0         2.0    >0 -> 1 Catalog residual map
OFFSET         -1.0         1.0    Cutoff level. 0-> None
DOMODEL        -1.0         1.0    > 0 => put solutions in a CC
                                   file with input image
OUTVERS        -1.0                MF table version number
                                      -1 => none, 0 => new
PRTLEV                             Level of messages desired:
                                     0 normal, 1 some extra,
                                     2 iteration results
PBPARM                             Primary beam parameters:
                                   (1) level to believe - <= 0
                                   means do not apply a primary
                                   beam (2) > 0 use (3)-(7)
----------------------------------------------------------------
JMFIT
Type: Task
Use:  JMFIT is a task to fit up to four Gaussian components to a portion
      of an image.  It can also solve for a constant, linear, or
      quadratic two-dimensional ``baseline'' surface.  The program
      estimates the error in the fits using the image rms and theory.
      The rms is determined from the image header keyword ACTNOISE (if
      present and positive) or from fitting the histogram of the full
      (or partial) image plane.

      The answers are returned in the input guess parameters GMAX,
      GPOS, GWIDTH and the uncertainties are returned in DOMAX, DOPOS,
      and DOWIDTH.
Adverbs:
  USERID......The ID of the owner of the images. 0 => current
              user, 32000 => any user.
  INNAME......First image name (name).      Standard defaults.
  INCLASS.....First image name (class).     Standard defaults.
  INSEQ.......First image name (seq. #).    0 => highest.
  INDISK......Disk drive # for the first image.  0 => any.
  BLC.........Bottom left corner of area of image to fit.
  TRC.........Top right corner of area of image to fit.
              Maximum area is 10000 pixels (100x100)
              If the data are a cube (BLC(3) > 1), then the plane number
              if displayed.  For "standard" spectral-line cubes (XYV)
              this will be the spectral channel.  NOTE: there is no
              requirement that the data be in standard coordinates or
              transposition - although am elliptical Gaussian in
              RA-velocity is a curious concept.
  OUTNAME.....Residual map name.            Standard defaults.
  OUTCLASS....Residual map class.           Standard defaults.
  OUTSEQ......Residual map seq. #.          0 => highest unique.
  OUTDISK.....Residual map disk no.      0 => highest with room.
  NGAUSS......The number of components to use in the fitting.
              0->1.  Maximum number is four.
  CTYPE.......Each component type.
              0=>1.  Two-dimensional elliptical Gaussian
              2=Solve for zero level
              3=Solve for zero level and slope
              4=Solve for zero, slope and curvature
              5=Set the six baseline parameters as desired.
                See EXPLAIN for use of GMAX,GPOS and GWIDTH
  GMAX........The peak value guess for each component.
              0=>Use value with largest absolute value in the
              BLC,TRC window
  GPOS........The position (X,Y) guess for components.  The values are
              in pixels in the order (X1,Y1,X2,Y2,X3,Y3,X4,Y4).
              0=>Use pixel location with largest absolute value.
  GWIDTH......The major axis, minor axis and position of major axis
              guess for components.  The values are pixels with
              degrees for position angle and the order is
              (MJ1,MN1,PA1,MJ2,MN2,PA2,...etc)  0->Use clean beam, if
              available; otherwise it will be taken as a circular
              gaussian of diameter 2.
  FMAX........The peak value fit for each component is returned.
  FPOS........The position (X,Y) fit for components is returned.  The
              values are in pixels in the order
              (X1,Y1,X2,Y2,X3,Y3,X4,Y4).
  FWIDTH......The major axis, minor axis and position of major axis
              fit for components are returned.  The values are pixels
              with degrees for position angle and the order is
              (MJ1,MN1,PA1,MJ2,MN2,PA2,...etc)
  DOMAX.......Flags for GMAX: if > 0 fit this parameter, else hold
              fixed.  Returned value is uncertainty in the fit
              parameter.
  DOPOS.......Flags for GPOS: if > 0 fit this parameter, else hold
              fixed.  Returned value is uncertainty in the fit
              parameter.
  DOWIDTH.....Flags for GWIDTH: if > 0 fit this parameter, else hold
              fixed.  Returned value is uncertainty in the fit
              parameter.
  BWSMEAR.....If > 0, the Clean beam will be smeared by a Gaussian in
              the radial direction of FWHM proportional to the radius
              from the pointing position times BWSMEAR.  Set it roughly
              to the channel bandwidth divided by the center frequency.
  RADIUS......If = 0, the rms used for error estimates is taken from
              the image header (keyword ACTNOISE) or found by fitting
              the full image plane.  If RADIUS > 0, then the rms is
              found by fitting only those pixels within RADIUS if the
              center of the BLC-TRC box.  If RADIUS < 0, then
              abs(RADIUS) is used as the rms.
  NITER.......The maximum number of iterations to use in the
              fitting. 0-> NGAUSS * 20.
  DOCRT.......<=0 -> Plot map, model and residual map and list
                fit information on the line printer
              >0 -> List fit info in message file only
  OUTPRINT....Disk file name in which to save the line printer
              output.   ' ' => use scratch and print immediately
              for interactive jobs - batch jobs use OUTPRINT =
              'PRTFIL:BATCHjjj.nnn' (jjj = job #, nnn = user #).
              When OUTPRINT is not blank, multiple outputs are
              concatenated and the file is not actually printed.
  DOOUTPUT....>0 -> Catalog residual map with fitted components
              and write them in a CC file attached to the output
              image.  If DOOUTPUT > 1.5, the components written
              to the CC file are not deconvolved from the beam.
  OFFSET......0-> Include all points in fitting area. Otherwise
              disregard all points less than OFFSET*MAX, where
              MAX is largest value in fitting window.  If MAX
              is less than zero, disregard all point greater
              than OFFSET*MAX.
  DOMODEL.....If true (> 0), put the deconvolved solutions in a
              new CC file attached to the input image.
  OUTVERS.....The results are written into an MF (Model Fit) table file
              with version OUTVERS unless OUTVERS is set < 0.  If
              OUTVERS points at a pre-existing table, the results are
              appended to the file.  OUTVERS = 0 always means to make a
              new MF table.
  PRTLEV......Print level: 0 => normal, 1 => some from fitting,
              2 => all from fitting (to diagnose problems)
  PBPARM......Primary beam parameters:
              (1) Lowest beam value to believe: 0 -> do not do the
                  primary beam correction.  This correction is done
                  to the printed parameters only.  The beam value used
                  is max (PBPARM(1), that computed from (2)-(7)).
              (2) > 0 => Use beam parameters from PBPARM(3)-PBPARM(7)
              (3-7)..For all wavelengths, the beam is described by the
                function:
                   1.0 + X*PBPARM(3)/(10**3) + X*X*PBPARM(4)/(10**7) +
                   X*X*X*PBPARM(5)/(10**10) + X*X*X*X*PBPARM(6)/(10**13)
                   X*X*X*X*X*PBPARM(7)/(10**16)
                where X is (distance from the pointing position in arc
                minutes times the frequency in GHz)**2.  For the VLA,
                these parms are, by default, given by Perley's fits:
                      0.0738 GHz  -0.897  2.71   -0.242
                      0.3275      -0.935  3.23   -0.378
                      1.465       -1.343  6.579  -1.186
                      4.885       -1.372  6.940  -1.309
                      8.435       -1.306  6.253  -1.100
                     14.965       -1.305  6.155  -1.030
                     22.485       -1.417  7.332  -1.352
                     43.315       -1.321  6.185  -0.983
                For the ATCA, these are by default:
                  1.5 GHz  -1.049   4.238  -0.8473  0.09073  -5.004E-3
                  2.35     -0.9942  3.932  -0.7772  0.08239  -4.429E-3
                  5.5      -1.075   4.651  -1.035   0.12274  -6.125E-3
                  8.6      -0.9778  3.875  -0.8068  0.09414  -5.841E-3
                  20.5     -0.9579  3.228  -0.3807  0.0       0.0
                See EXPLAIN PBCOR for details
----------------------------------------------------------------
JMFIT: Task to fit Gaussian models to an image by least squares.
DOCUMENTORS: E.B.Fomalont NRAO/VLA, F. Schwab NRAO, CV
RELATED PROGRAMS : SLFIT,IMEAN,MAXFIT,IMFIT,SAD

                      PURPOSE

     JMFIT fits up to four Gaussian-shaped components to a
selected part of an image.  One of the components can be a
baseline function with a zero level, slope and curvature term.
JMFIT is most commonly used to derive the position, peak and
integrated intensity and angular size of a source which is not
too extended.  An initial guess for the parameters, some of
which are picked as defaults, must be supplied before running
the task.  Solution and error estimates are generated and the
residual image after the fit can be printed on the line-printer.
An arbitrary selection of parameters may be held constant in
the solution.  Default initial guesses chosen by the task may
suffice in the case of a single-component model; otherwise,
sufficiently good guesses for at least the positions must be
provided by the user.
     The fitting algorithm is one which is based on a
algorithm of Davidon.  Occasionally, the solution will converge
on an obviously unacceptable fit.  If this occurs, then try a
better initial guess.  When fitting several components to a
blobby source, the fitted parameters may be absurd.  Careful
selection of fixed parameters then is necessary.

                COMMENTS ABOUT SOME PARAMETERS

BLC, TRC:
     The fitting area should be chosen as small as possible; and
several disconnected components should be fit separately.  The
fitting area is limited to an area of 10000 pixels.  If the data are a
cube (BLC(3) > 1), then the plane number if displayed.  For "standard"
spectral-line cubes (XYV) this will be the spectral channel.

NGAUSS:
     The number of components to fit.  The maximum number is
four, and 0->1.

CTYPE:
     The component types, placed in a scalar array of length 4.
     0->1 Elliptical Gaussian component.
        2 Zero level.
        3 Zero level and slope.
        4 Zero level, slope and curvature.
        5 Insert baseline parameters as follows:
          GMAX = zero level
          GPOS = slope (intensity per pixel),
                 Orientation (deg, measured from N through E)
          GMAX = Curvature (intensity per square pixel),
                 Ellipticity of curvature (-1 to 1),
                 Orientation (deg, measured from N through E)
     DOMAX, DOPOS and DOMAX are used to hold parameters fixed
     for Types 1 and 5 only.

GMAX:
     The initial guess of the model intensity may be supplied.
The units of GMAX are the same as those in the map.  The
default of 0 will place the most extreme value in the fitting
area (negative or positive) in GMAX for the first component. Any
subsequent components with 0 default will be given the value of
0.1 times the extreme value.

GPOS:
     The initial guess of the model position.  The location must
be given by a pair of pixel coordinates.  The default of 0 will
insert the location of the extreme value for all components.
Note that GPOS has the meaning of slope and orientation for a
baseline component.

GWIDTH:
     The defaults for the component widths are generally
reasonable: either the clean beam size or a circular beam of two
pixels FWHP.  To try to get the minimization algorithm off to a
faster start for circular Gaussian models, the task will
introduce a slight ellipticity before beginning the fitting.
This is not done if either axis is to be held constant.

NITER:
     The maximum number of iterations, NITER, is defaulted to 20
times NGAUSS if the input value is 0.   Whenever convergence is
not achieved in NITER iterations, JMFIT so informs the user and
ceases execution.  This is to discourage the user from accepting
a poor solution.  When convergence is not achieved, the user
should increase NITER if it appears that the program may have
been converging on an acceptable solution, but doing so slowly.
Otherwise, one should alter the initial guesses, try fitting to
data within a somewhat smaller or larger window, or alter the
value of OFFSET.

DOWIDTH:
     In fitting complicated sources, it is common to hold some
of the component diameters fixed.

DOCRT:
     Set DOCRT = -1 in most cases.  This produces an automatic
hard copy of the solutions and a digital map of the input image,
the initial guess and the post-fit residuals. These maps are
useful for assessing the quality of the fit.  When DOCRT = -1,
the progress of the (iterative) minimization algorithm is
reported on the line-printer output.  Otherwise, it is displayed
on the terminal.

DOOUTPUT:
     Catalog the post-fit residual map.  If all parameters are
held fixed, no fitting is done and a residual map is generated
(This feature, one which is present in IMFIT, isn't implemented
yet in JMFIT - F.S. 3/7/85).  Only the fitted area is cataloged.
A CC file is also written with this output image listing the
components.  If DOUTPUT > 1.5, the components written are not
deconvolved from the beam.  If 0 < DOOUTPUT <= 1.5, the
components are deconvolved if possible.

OFFSET:
     This adverb permits the exclusion of low valued points when
doing the fit.  If the extremum value (MAX) in the fitted area
is positive, then all points less than OFFSET*MAX are ignored in
the fit.  If the extremum value in the fitted area is negative,
then all points greater than OFFSET*MAX are ignored in the fit.
If OFFSET = 0, then all points are used.

                 COMMENTS ON THE USE OF JMFIT

SETTING UP THE PARAMETERS
     For some simple cases the defaults set by JMFIT are good
enough starting values.  Some examples are as follows: (Always
insert the appropriate input map and always set BLC and TRC to
the smallest area needed for the fit.  The verb TVWINDOW can be
used to set the window.)
    1.  Fit to one Gaussian component.
        Nothing to specify except flags.
    2.  Fit to one Gaussian and zero level
        NGAUSS=2;CTYPE=1,2;GO
        Additional flags can be specified
    3.  Fit to one Gaussian, zero-level and slope
        NGAUSS=2;CTYPE=1,3;GO
        Additional flags can be specified
    4.  Fit to two Gaussians and a zero level
        NGAUSS=3;CTYPE=1,1,2;GO
        Depending on the source complexity it may be
        important to set some of the fitting flags
    5.  Fit one Gaussian with a zero level, slope and
        curvature only in E/W direction
        NGAUSS=2;CTYPE=1,5;DOMAX=1,1;DOPOS=1,1,1,-1
        DOWIDTH=1,1,1,1,-1,-1;GPOS=0,0,0,90;
        GWIDTH=0,0,0,0,0,90

FLUX DENSITY DETERMINATION:
     When attempting to obtain the flux density of a
well-resolved source, the task IMEAN, which integrates the map
values in a specified rectangle, is often more accurate than
fitting the source with several Gaussian components and summing
the integrated flux densities.

PEAK FLUX DENSITY DETERMINATION:

     The verb MAXFIT, a simple fitting of the peak of a
component with a second degree interpolation, is much faster
than JMFIT and useful to obtain the approximate peak and
position of a component.

ERRORS OF PARAMETERS:
     The error estimates should be regarded as tentative. An estimate of
the each error is determined from theory based on the actual rms of the
image (neglecting signal portions) or the rms given in header parameter
ACTNOISE (if present and positive).  AIPS verb ACTNOISE may be used to
set this header parameter.  See the explain file for SAD for details.
Theory gives expressions for the errors in two limiting cases: point
source (the beam area > 0.9*the fitted gaussian area) and expanded
source (the beam area < 0.1*the fitted gaussian area).  The formulae are
taken from J. Condon paper 'Errors in Elliptical Gausian Fit', AA, 1996.
The intermediate case is handled by interpolation between the two limit
cases. It is not clear if this interpolation provides a good estimation
of the error in general relation of beam and fitted gausian sizes and
position angle.

DECONVOLUTION:
     When fitting to a clean map, JMFIT deconvolves the clean
beam from the fitted component size.  The nominal deconvolution
is obtained by deconvolving the fit from the clean beam.  A
value of 0.0 means that the source is smaller than the clean
beam in some dimension.  The minimum and maximum values are
obtained by deconvolving the source beam parameters with all
combinations of 0.7 * error and listing the extreme values.

OUTPUT FROM THE MINIMIZATION SUBROUTINE
     With iterative algorithms, it is important to monitor the
progress of the iteration.  For this reason JMFIT, at each
iteration, prints out the value of the chi-squared function
(the sum of the squares of the residuals) and the Euclidean
norm of the gradient of the chi-squared function (it is
trying to zero the gradient).  chi-squared is labeled F in
the printout, and the gradient norm is labeled simply
GRADIENT.  The residuals are measured in the same units as
the units of the input data, but the units are not printed out
by JMFIT.
     Also, at each iteration, the current estimates of the
solution parameters are printed out in their natural order
(i.e., amplitude, x-position, y-position, major axis fwhm,
minor axis fwhm, and position angle for Gaussian components),
component-by-component.  (Fixed parameters are not included.)
The units are input data amplitude units for amplitudes,
pixel coordinates for positions, pixels for widths, and
radians for angles.  For the usual orientation of the RA and
DEC axes, pi/4 is added to the position angles in these
printouts.
     PRTLEV is provided to suppress some or all of this output.  When
DOCRT = -1 this information is printed on the line-printer.  Otherwise,
it is printed on the terminal.  Thus, if the user is annoyed by the
quantity of output which is written to the terminal, he can set DOCRT =
-1 to have it instead written on the line-printer (and vice-versa, of
course).  If PRTLEV < 2, the output is greatly reduced.

