; IMFIT
;---------------------------------------------------------------
;! Fits Gaussians to portions of an image
;# TASK ANALYSIS MODELING
;-----------------------------------------------------------------------
;;  Copyright (C) 1995-1999, 2001, 2003-2005, 2008-2011, 2013-2014
;;  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
;-----------------------------------------------------------------------
IMFIT     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
IMFIT: Task to fit gaussian models to an image
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 IMFIT
                                   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
FSHIFT   @                         RA/Dec shift to put comp 1
         @                         on nearest pixel
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 iteration
                                     0->200*NGAUSS
                                   Solve for model parameters?
DOPRINT        -4.0       132.0    <0 -> Print maps and
                                     solutions on line printer
                                   >0 -> print on terminal
                                   =0 -> maps not printed
FITOUT
                                   Disk file to save fit info
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
STVERS         -1.0                STars output file version
                                      -1 => none, 0 => new
PBPARM                             Primary beam parameters:
                                   (1) level to believe - <= 0
                                   means do not apply a primary
                                   beam (2) > 0 use (3)-(7)
EFACTOR         0.0        20.0    Scale width sigmas in the
                                   deconvolution tests
----------------------------------------------------------------
IMFIT
Type: Task
Use:  IMFIT is a task to fit a portion of an image with up to four (4)
      gaussian components.  It can also fit for a baseline term of up to
      second order.
      The program estimates the error in the fits 2 ways: using the
      quality of the fit and using the image rms and theory.  The rms
      for the latter is determined from the image header keyword
      ACTNOISE (if present and positive) or from fitting the histogram
      of the full (or partial) image plane.  Note that pixels which
      are exactly zero are not used in this fit, allowing balnked
      pixels to be REMAGed to zero.

      The answers are returned in the input guess parameters GMAX,
      GPOS, GWIDTH and the uncertainties are returned in DOMAX, DOPOS,
      and DOWIDTH.
Adverbs:
  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)
  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.
*********   The following are output adverbs:
  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)
  FSHIFT......The values that one could put in RASHIFT and DECSHIFT
              to move the first component to the nearest integer
              pixel.  These values include any previous shifts and are
              corrected for rotation.
*********   The following are input/output adverbs:
  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.
*********   The following are input adverbs:
  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.
              The Clean beam is used as an initial estimate of the
              source size and is used in the deconvolution attempt to
              find the true size from the fit size.  The peak
              intensity printed will be corrected for this effect.
  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 of the
              center of the BLC-TRC box.  If RADIUS < 0, then
              abs(RADIUS) is used as the rms.  Pixels which are
              exactly zero are not used in the fitting for rms.
  NITER.......The maximum number of iterations to use in the
              fitting. 0-> NGAUSS * 200.
  DOPRINT.....<0 -> Plot map, model and residual map and list fit
                    information on the line printer
                    When FITOUT is not blank, DOPRINT=-2 suppresses
                    the page-feed character on page headers and
                    DOPRINT=-3 suppresses page headers and most other
                    header information.
                    DOPRINT=-4 produces output to FITOUT with one line
                    per source component in a special CfA format.
              =0 -> List fit info in message file only
              >0 -> Plot map, model, and residual map and list fit
                    information on the terminal
  FITOUT......Disk file name in which to save the line printer output.
              ' ' => use scratch and print immediately for interactive
              jobs - batch jobs use FITOUT = 'PRTFIL:BATCHjjj.nnn'
              (jjj = job #, nnn = user #).  When FITOUT 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.
  STVERS......The results may also be written to a STars table
              attached to the image.  -1 => do not do this, 0 => write
              a new one, > 0 => add to the specified version.
  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)
                  Otherwise use default parameters for the VLA (or
                  ATCA where appropriate)
              (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
  EFACTOR.....The range of possible deconvolved widths is found by
              trying the deconvolution with Major, Minor, and PA each
              -EFACTOR*sigma, 0, +EFACTOR*sigma from the fit value.
              The highest value found in the 27 tests and the lowest
              value found are reported as well as the values found at
              the nominal (fit) values.  0 -> 1.3.
----------------------------------------------------------------
IMFIT: Task to fit Gaussian models to an image.
DOCUMENTOR: E.B.Fomalont NRAO/VLA
RELATED PROGRAMS : SLFIT,IMEAN,MAXFIT,JMFIT,SAD

                      PURPOSE

     IMFIT 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.  IMFIT 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.
     The fitting algorithm is based on the subroutine LMDER in the
Argonne National Laboratory Subroutine Package and the algorithm uses
a linearized least-square solution to obtain the parameters.
Occasionally, the solution will converge on an obviously unacceptable
fit.  If this occurs when fitting one component to the source, try a
better first guess.  When fitting several components to a blobby
source, the fitted parameters may be absurd.  Careful selection of
fixed parameters will then be 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.

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 N thru E)
          GMAX = Curvature (intensity per pixel**2),
                 Ellipticity of curvature (-1 to 1),
                 Orientation (deg N thru 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 of a circular beam of two pixels FWHP.
Because of poor convergence properties of the algorithm for circular
Gaussian models, the task will introduce a slight ellipticity before
beginning the fitting.  This is not done if either axes is held
constant.

NITER:
     The number of iterations, NITER, is defaulted to 200*NGAUSS if it
is set to 0.   If you are somewhat unsure if your model is reasonable
or is converging to an acceptable solution, especially for fits to
complicated sources with several Gaussians, set NITER=50 and check how
the convergence is going.  The task has several other termination
conditions.  If the solution is unchanging to a level of about 0.1%,
it will terminate.  It some of the fitting parameters are obviously
ridiculous, it will also terminate.

DOWIDTH:
     In fitting complicated sources, it is common to hold some of the
component diameters fixed in order to obtain reasonable convergence.

DOPRINT:
     Set DOPRINT = -1 in most cases.  This produces an automatic hard
copy of the solutions and a digital map of the input image, the first
guess and the residuals. These maps are most useful for determining
the validity of the fit.

DOOUTPUT:
     Catalog the residual map after the fit has been subtracted.  If
all parameters are held fixed, no fitting is done and a residual map
is generated.  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 IMFIT

SETTING UP THE PARAMETERS
     For most simple cases the defaults in IMFIT adequately provide
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
        Note that the zero level information is associated
        with the second 'gaussian' with CTYPE=2
        Additional flags can be specified.
    3.  Fit to one Gaussian, zero-level and slope
        NGAUSS=2;CTYPE=1,3;GO
        Note that the zero level and slope information is
        associated with the second 'gaussian' with
        CTYPE=3.
        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 IMFIT and useful to
obtain the approximate peak and position of a component.

ERRORS OF PARAMETERS:
     An estimate of the each error is determined from theory based on
the actual rms (R) of the image (neglecting signal portions) or the
rms given in header parameter ACTNOISE (if present and positive).
AIPS task IMEAN or verb ACTNOISE may be used to set this header
parameter.  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.  The formulae now used are:
     M = 1                           (Clean beam area > 0.9 fit area)
     M = sqrt (8 * ARbeam/ARimag)    (Clean beam area < 0.1 fit area)
     M = sqrt (0.8 + 0.25*(ARbeam/ARimag-0.1))    (else)
Then
     Delta(P) = M * R
     Delta(W) = W * Delta(P) / P
     Delta(PA) = sqrt(2) * (Smaj+Smin)/(Smaj^2+Smin^2) * Delta(P) / P
     Delta(X) = sqrt[(Delta(Smaj)*sin(PA))^2 + (Delta(Smin)*cos(PA))^2]
                    / sqrt (8 * ln (2.0))
     Delta(Y) = sqrt[(Delta(Smaj)*cos(PA))^2 + (Delta(Smin)*sin(PA))^2]
                    / sqrt (8 * ln (2.0))
     Delta(F) = Delta(P) * ARimag/ARbeam * sqrt(1+2*ARbeam/ARimag)
where P is fit peak, W = Smaj or Smin are fit widths, X and Y are
positions, PA is position angle, ARimag = Smaj*Smin, ARbeam=Bmaj*Bmin.

DECONVOLUTION:
     When fitting to a clean map, IMFIT deconvolves the Clean beam from
the fitted component size.  The nominal deconvolution is obtained by
deconvolving the fit from the Clean beam (corrected for bandwidth
smearing).  A value of 0.0 means that the source is smaller than the
corrected Clean beam in some dimension.  The minimum and maximum
values are obtained by deconvolving the source beam parameters with
all 27 combinations of EFACTOR * (-1, 0, 1) * uncertainties in the
major axis, minor axis, and position angle.  The extremas in these
parameters over all 27 tries are listed.  The default EFACTOR (1.3)
appears to work well with the considerations below.

     An estimate is given concerning whether the source should be
viewed as resolved or unresolved.  The task assumes that the
component is probably unresolved if:
  (a) the deconvolution of the fit answers has the major axis 0
OR
  (b) the fit total flux minus the error in the total flux is
      less than the peak  AND  the minimum deconvolved major axis
      is 0.
The task assumes the conponent is probably resolved if
  (c) the total fit flux minus the error in the total fit flux is
      greater than peak flux  AND  the minimum deconvolved major axis
      is greater than 0.
The task is undecided about resolution if
  (d) the deconvolution of the fit answers has the major axis greater
      than zero.
AND
  (e) the fit total flux minus the error in the total flux is
      less than the peak  BUT  the minimum deconvolved major axis
      is greater than 0.
OR
  (f) the total fit flux minus the error in the total fit flux is
      greater than peak flux  BUT  the minimum deconvolved major axis
      is 0.
Note: the total flux and its error are corrected for primary beam and
the peak is corrected for primary beam and bandwidth smearing (if such
corrections are requested) in the tests described above.  If the
component is unresolved, the best estimate of its total flux is its
peak brightness.  In cases where the task is uncertain, use caution in
deciding if the component is resolved.  Noise seems preferentially to
make sources appear resolved when they are not.  Note too that a
"resolved" source may be clearly unresolved along the minor axis.

