; TVSAD
;---------------------------------------------------------------
;! Finds and fits Gaussians to portions of an image with interaction
;# TASK ANALYSIS MODELING INTERACTIVE TV
;-----------------------------------------------------------------------
;;  Copyright (C) 2014, 2016
;;  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
;-----------------------------------------------------------------------
TVSAD     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
TVSAD:   Interactive task fits 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 #
INVERS         -1.0     46655.0    Input image MF ext. version
                                      <= 0 => new
BLC             0.0      4096.0    Bottom Left corner of fit
TRC             0.0      4096.0    Top Right corner of fit
DORESID        -1.0         2.0    >0 => 1 Catalog residual map
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     40000.0    Max. Number of components
CPARM           0.0      9999.0    Flux cutoff levels in
                                   descending order
ICUT            0.0      9999.0    Retry level
BWSMEAR        -0.1         0.1    Bandwidth smearing corr.
SORT                               Sort Order of output ' '=RA
DOCRT          -4.0       132.0    >0 => print fits on CRT
                                   > 72 => terminal width
                                   =0 => print suppressed
                                   <0 => line printer
FITOUT
                                   Disk file to save fit info
OUTVERS        -1.0     46655.0    CC file version # for
                                   solutions: -1 => none, 0 new
STVERS         -1.0                STars output file version
                                      -1 => none, 0 => new
DOALL          -1.0         1.0    >0 => fit multiple peaks
DOWIDTH        -1.0         2.0    >0 => fit widths
GAIN            0.0         1.0    Amp-dependent part of retry
                                   and warning levels
DPARM                              Reject components with
                                   (1) peak < D(1)
                                   (2) flux < D(2)
                                   (3) rms > D(3) + GAIN * flux
                                       (in quadrature)
                                   (4) widths > D(4) cells
                                   (5) peaks > D(5) cells
                                          outside interior area
                                   (6) peaks > D(6) cells
                                          outside image
                                   (7) total residual flux in
                                          the fitting box
                                          [no default]
                                   (9) Use IN2NAME etc to set
                                   image of rms
                                    >= 2 => CPARM is in S/N
                                        (see HELP)
                                   (10) Expand fit box pixels
IN2NAME                            RMS image name (name)
IN2CLASS                           RMS image name (class)
IN2SEQ          0.0      9999.0    RMS image name (seq. #)
IN2DISK         0.0         9.0    RMS image disk drive #
PRTLEV          0.0      1023.0    Debug print level
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
----------------------------------------------------------------
TVSAD
Type: Task
Use:  TVSAD and SAD (Search And Destroy) are tasks to find a number of
      potential sources in a (sub)image, fit Gaussian components to
      these, and optionally print the results, store them in a CC
      file, and create a residual file.  TVSAD is an interactive task
      which enables you to look at each island before and after it is
      fit to determine if one or more than one (or no) Gaussian should
      be fit to that portion of the image.  The searching may be
      restarted with different parameters to remove additional
      sources.  Only one plane of a cube is fit at one time.  However,
      the MF file keeps track of the plane number and may contain data
      for more than one plane.  With additional tasks (to be written
      yet), one should be able to use this feature to do spectral
      fitting of the objects found spatially by TVSAD and/or SAD.

      The printer display produced by TVSAD and SAD shows the results
      of the fitting, estimates of the uncertainties in the answers
      (based on the actual noise and theory rather than the quality of
      the fitting), and attempts to deconvolve the Clean beam from the
      fit sources.  The display contains a mark character near the
      beginning of the row.  If the scale for a particular peak or
      flux has to be reduced by 1000 an * will appear as the mark and
      if it must be reduced by 1000000 an ! will appear as the mark.
      When the residual image near a source exceeds a limit controlled
      by the cutoff (CPARM) and GAIN parameters, an H, L, or S will be
      shown as the mark character if the numbers are not scaled.  H is
      for a high point in the residual, L for a rather negative point,
      and S for a large rms.  If the width of the display allows, the
      max residual is displayed with the deconvolution results and an
      * is placed near excessive residuals.

      If a primary beam correction is applied and the source lies
      outside the primary beam limit (PBPARM(1)), the uncertainty in
      the peak flux is surrounded by ! characters rather than
      parentheses.

       Details of the interactive operation and options are described
       in the EXPLAIN file.  See also AIPS Memo 119, "TVSAD:
       interactive search and destroy", January 2016 for a more
       detailed description of the use of this task.
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.
  INVERS......Version number of a MF table extension file attached to
              the input image and copied to the ouput residual image (if
              any).  <= 0 => create a new one, > 0 => use a pre-existing
              one and subtract the components in it before searching for
              more.
  BLC.........Bottom left corner of area of image to fit.
  TRC.........Top right corner of area of image to fit.
  DORESID.....>0 -> Create a residual map of subimage
  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 maximum number of components to search for
              0->10.   Maximum number is 40000.
  CPARM.......Search for potential sources down to these levels in
              units of the input image or in units of S/N (when
              DPARM(9) >= 2.0 and IN2NAME et al are specified). No
              defaults - list in descending order, A value of 0 => no
              more levels.  Thus CPARM = 8, 4, 0.666, 0  => 3 rounds
              of searching.  It has been found useful to fit the
              higest levels first, then look somewhat lower, and then
              even lower.  Low level signal around a strong source can
              sometimes confuse the fitting.  CPARM all 0 -> 3,0 in
              S/N or 3*rms,0 in image units.
  ICUT........If the peak residual in an island, after fitting a single
              Gaussian exceeds ICUT in image or S/N units (as CPARM)
              adjusted in quadrature by GAIN*flux, then the program
              will attempt a 2-Gaussian fit to the island.   0 =>
              CPARM(i) for each search level i.  Note that TVSAD can
              fit up to 8 Gaussians in any one island, but it has to
              be directed manulayy to do so.
  BWSMEAR.....If not 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 abs(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; set BWSMEAR < 0 to suppress this
              correction.
  SORT........Sort order of output.
              ' ' -> 'R'-> Right Ascension
              'D' -> Declination
              'X' -> X pixel (~RA if no rotation)
              'Y' -> Y pixel (~Dec if no rotation)
              'S' -> Flux
  DOCRT.......False (< 0) use the line printer if FITOUT = ' '
                   else write named FITOUT file only.
                   -2 causes the line feed characters in FITOUT to
                      be suppressed
                   -3 causes the documentation lines to be suppressed
              Null  (= 0) print suppressed.
              True  (> 0) use the terminal interactively.  The task will
                   use the actual terminal width as a display limit
                   unless 72 < DOCRT < width.  In that case, the display
                   limit will be DOCRT characters.
  FITOUT......Disk file name in which to save the fit (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.
  OUTVERS.....= 0 => Create a new CC file on input image.
              < 0 => do not write any CC file
              > 0 => append components to pre-existing CC file
              SAD writes the deconvolved widths (0 on failure of
              deconvolution) into the CC files, if the input image is a
              Clean image.  If DOWIDTH <= 0, then it writes point
              components  into the CC file.
  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.
  DOALL.......>0 -> If an island has multiple (up to 4) peaks, fit
              multiple gaussians.  Otherwise fit 1 gaussian to each
              island.
  DOWIDTH.....<= 0     -> Fit all sources with Clean beam
              >0 & <=1 -> Let fit find size of source starting with the
                          Clean beam
              >1       -> Let program guess the sizes and then fit them
                          (this can be somewhat unstable)
  GAIN........Retry, cutoff and warning levels are computed by
                 SQRT ( level**2 + (GAIN*flux)**2 )
              where the level is the no-signal retry (ICUT), cutoff
              (DPARM(2)) or warning (CPARM cutoff) level and flux is
              the peak of the component in question.  No default.
  DPARM.......Components are fit within islands that are rectangular in
              shape.  It may be possible to get bad solutions, so the
              program will reject components with:
              (1) peak <= DPARM(1).  No default; 0 is okay if you
                  mean it.  Units are image units unless DPARM(9)>=2,
                  in which case they are in S/N.
              (2) flux <= DPARM(2).  No default; 0 okay
              (3) rms in the box > sqrt (DPARM(3)**2 + (GAIN*flux)**2).
                  0 => 1000.  Units are image units unless DPARM(9)>=2,
                  in which case they are in S/N.
              (4) widths > DPARM(4) cells.   0 => 1000.
              (5) peaks > DPARM(5) cells outside the island (fitting
                  area) but interior to the image.  No default;
                  < 0 => inside the fitting area required.
              (6) peaks > DPARM(6) cells outside image.  No default.
                  This allows SAD to fit components legitimately on the
                  edge of the whole image while rejecting spurious ones
                  extrapolated to positions outside the fitting
                  rectangle but fully in the image.
              (7) total residual flux in the fitting box >
                  SQRT (DPARM(7) + (GAIN*flux)**2).  No default.
     ***      (9) > 0 => use IN2NAME, IN2CLASS, etc as an image of the
                  RMS so that it may vary with position.  = 1 => use
                  the image only to evaluate the uncertainties in the
                  fit components.  >= 2 => use the image also to
                  evaluate the noise at each pixel as data are read.
                  In that case CPARM, ICUT, DPARM(1), and DPARM(3) are
                  in S/N units.
     ***      (10) Expand (or contract) the box used to fit the
                   component.
              Whem a component is rejected, TVSAD will resturn to the
              display of the island with its initial guess(es).  You
              may then try some other initial guess or reject the
              island entirely.
  IN2NAME.....RMS image name (name).        NO DEFAULT.
              The task RMSD may be used to prepare such an image.
              Note that the RMS image must be at the same coordinates
              more or less as the INNAME image but they do not have to
              allign.  The RMS is found by looking up the coordinate
              of the coponent in the RMS image.
              If DPARM(9) >= 2, the RMS image is ASSUMED to have the
              same size and coordinates as the input image.  Bad
              things will happen, perhaps not in obvious ways, if this
              is not true.  In this case, the noise image is used
              throughout rather than simply in evaluating the
              uncertainties in the component parameters.
  IN2CLASS....RMS image name (class).       Standard defaults.
  IN2SEQ......RMS image name (seq. #).      0 => highest.
  IN2DISK.....RMS drive # for the first image.  0 => any.
  PRTLEV......To find out which islands are rejected by the above for
              what reasons, set PTRLEV to the sum of the desired
              reasons: 1 peak, 2 flux, 4 rms, 8 widths, 16 X position
              outside interior window, 32 Y position outside interior
              window, 64 X position outside image, 128 Y position outside
              image., and 256 for residual flux.  511 does all of them.
              If you don't want to see the low peak ones, set PRTLEV
              to an even number.  Any that are rejected for low peak
              will not be displayed even if they are also rejected for
              a selected reason.
  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.
                See explain 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.
----------------------------------------------------------------
TVSAD: Task to find potential sources and fit Gaussian models to an
       image by least squares with TV interactivity
DOCUMENTORS: W. Jaffe, Leiden
RELATED PROGRAMS : SAD, MFPRT, SLFIT, IMEAN, MAXFIT, IMFIT, JMFIT

TV interactivity in TVSAD is essentially the only difference between
SAD and TVSAD.  The task finds all islands of emission above CPARM(i)
and then loops to fit each island in the order of peak flux.  TVSAD
makes an initial guess of the Gaussian components.  It then invites
you to look at an interpolated image of the island and surrounding
area.  You may adjust the island boundaries and change the initial
guess beofre asking TVSAD to do the fit.  You can also tell it to omit
fitting this island and to go on to the next.  If the fit fails for
any reason, it will loop back to the TV display, turning that display
back on if you had chosen to do the fitting non-interactively.  If
the fit succeeds, it will display the island image with the fit
components removed.  You may tell the task to go on to the next
component or to re-try this island using the initial TV display
routine again.  When all islands above CPARM(i) are fit or omitted,
TVSAD will repeat the process with CPARM(i+1) if it is also greater
than zero.

The TV menu in the main interactive routine has two columns.  The left
column offers options:

---------------
| OFF TRANS   |    reset color and black-and-white enhancements
| OFF TVZOOM  |    reset TV zoom
| TVFIDDLE    |    usual zoom and enhancement verb
| CURVALUE    |    display image intensities under TV cursor
| ZOOM IN     |    increase the interpolation factor by one
| ZOOM OUT    |    decrease the interpolation factor by one
| REBOX       |    change the island boundaries
| ENTER GUESS |    enter a guess for up to 4 Gaussians with TV
| REDO GAUSS 1|    enter a guess only for Gaussian number 1
| REDO GAUSS 2|    enter a guess only for Gaussian number 2
| REDO GAUSS 3|    enter a guess only for Gaussian number 3
| REDO GAUSS 4|    enter a guess only for Gaussian number 4
| REDO GAUSS 5|    enter a guess only for Gaussian number 5
| REDO GAUSS 6|    enter a guess only for Gaussian number 6
| REDO GAUSS 7|    enter a guess only for Gaussian number 7
| REDO GAUSS 8|    enter a guess only for Gaussian number 8
---------------

Only a suitable number of REDO options are offered.

The right-hand columns has options:

---------------
| DO FIT      |    try the fit for this island
| TVOFF       |    try the fit for this island and turn off
                   interactivity until there is a failure
| NEXT ISLAND |    forget this island, go on to the next
| QUIT        |    store all results, do reports, and exit task
---------------

The post-fit residual display offers the options in the left-hand
column

---------------
| OFF TRANS   |    reset color and black-and-white enhancements
| OFF TVZOOM  |    reset TV zoom
| TVFIDDLE    |    usual zoom and enhancement verb
| CURVALUE    |    display image intensities under TV cursor
| ZOOM IN     |    increase the interpolation factor by one
| ZOOM OUT    |    decrease the interpolation factor by one
---------------

and in the right-hand column

---------------
| GOOD        |    accept this result and go on to the next island
| RE-TRY      |    return to the main TV interaction to try this
                   island over again
| NEXT ISLAND |    forget this island, go on to the next
| QUIT        |    store all results, do reports, and exit task
---------------


      -------------- Explain file for SAD --------------------

                      PURPOSE

SAD attempts to find all sources in a subimage whose peak is brighter
than a given level.  It searches the subimage specified by BLC and TRC
for all points above this level and merges such points in contiguous
"islands".  For each island, initial estimates of the strength and size
are generated.  Then the gaussian fitting algorithm used in JMFIT is
called to determine the least square guassian fit.  If DOALL is < 0.,
only one gaussian is fit per island, with initial estimates generated
from 2nd moments.  If DOALL is TRUE and multiple peaks (above a cutoff
= CPARM(i)) are found within the island, then up to 4 gaussians are
fit.  If the peak residual of a single Gaussian fit exceeds ICUT then
a more complex 2-Gaussian fit ai attempted.  The results of the fit
are printed on the screen or the line printer and written into a
ModelFit (MF) extension table.  If DORESID is TRUE, the gaussians are
subtracted from the input subimage, and a residual imageis cataloged,
using OUTDISK, OUTNAME etc. If DOWIDTH(1) is negative, a gaussian with
fixed size = CLEAN beam will be fit.

                WARNING!!!

The list generated by SAD is not COMPLETE in any statistical sense.  You
will be warned about islands where the fitting algorithm failed.  Note
that sizes etc. of sources near the noise are very unreliable.  For
something approaching completeness Walter Jaffe recommends:

    Start with a well Cleaned map
    Set DOWIDTH = FALSE and only fit point sources.
    Do a series of search on progressively more tapered maps and
       compare results to determine flxes of extended sources.
    Make residual maps and inspect them to see where SAD blew it.

Suggestion from Eric Greisen:
    DOWIDTH false always produces poor results even for only slightly
extended objects.  For decent S/N, the width fitting is quite reliable.

                COMMENTS ABOUT SOME PARAMETERS

BLC, TRC:
     These specify the area to be searched for sources.  If an
residual file is requested, it covers the entire input image.

NGAUSS:
     The maximum number of sources to find.  Currently limited
to 40000.

CPARM:
     Search for source peaks down to this level.  If specified as 0.,
SAD will determine the real (i.e. source-free) R.M.S. in the full
image area and take CPARM(1) to be 3.0 times this value.  If you make
some CPARM(i) too small, you'll get lots of nonsense.  The task uses
CPARM(1), then restarts using CPARM(2), etc so long as the CPARM(i) >
0.  This lets you extract the strongest points first without low-level
interference and then go after weaker things.

DOCRT:
     DOCRT = 1. will cause a list of fluxes, positions, and extension
parameters, with their estimated errors, to be displayed on your
terminal.

     DOCRT = -1. will cause the list to be printed on the line printer.

     DOCRT = 0. causes listing to be supressed.

DORESID
     Catalog the post-fit residual map.

DOWIDTH
     If you think most sources are unresolved, use DOWIDTH(1) = -1.0.
This will fit simple point sources to each peak (and do very badly on
those few that are slightly resolved).

PRIMARY BEAM CORRECTION

     TVSAD corrects an image for the primary beam attenuation of
the antennas.  The function used to model the primary beam for normal
VLA frequencies

            F(x) =  1.0
                   + parm(3) * 10E-3  * x
                   + parm(4) * 10E-7  * x*x
                   + parm(5) * 10E-10 * x*x*x
                   + parm(6) * 10E-13 * x*x*x*x
                   + parm(7) * 10E-16 * x*x*x*x*x

where x is proportional to the square of the distance from the
pointing position in units of [arcmin * freq (GHz)]**2, and F(x)
is the multiplicative factor to divide into the image intensity at the
distance parameter x.  For other antennas, the user may read
in appropraite constants in PBPARM(3) through PBPARM(7).  The
flag, PBPARM(2) must be set to a positive number to invoke this
option and PBPARM(3) must not be zero.
     This correction scales with frequency and has a cutoff
beyond which the map values are set to an undefined pixel value GIVEN
in PBPARM(1).  At the VLA frequencies the default cutoff is
                 1.485 GHz     29.8  arcmin
                 4.885 GHz      9.13 arcmin
                15     GHz      2.95 arcmin
                22.5   GHz      1.97 arcmin
and occurs at a primary beam sensitivity of 2.3% of the value at
the beam center.  Corrections factors < 1 are forced to be 1.
The estimated error of the algorithm is about 0.02 in (1/F(x))
and thus leads to very large errors for x>1500, or at areas
outside of the primary response of 20%.  The cutoff level
may be specified with DPARM(1).

Default values of PBPARM for the VLA are 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
For the Karl G Jansky VLA ("EVLA"), the defaults are frequency
dependent.  If the observing frequency is between two tabulated
frequencies, then the beam is computed for each of the tabulated
frequencies and then interpolated to the observing frequency.  The
values used are far too numerous to give here, see EVLA Memo 195,
"Jansky Very Large Array Primary Beam Characteristics" by Rick Perley,
revision dated June 2016.  Obtain it from
http://library.nrao.edu/evla.shtml

                 COMMENTS ON THE USE OF SAD

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.  Of
course for a well resolved source a gaussian fit is only a crude
approximation in any case.  TVSTAT allows you to use an irregular area
for this determination.

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, but more
reliable than those from formal standard error estimates.  The latter
are meaningless since the objects are normally well fit using only 1-4
Clean beam areas.  This means that the number of parameters fit exceeds
the number of truly independent samples (neglecting the samples of 0
which were part of the selection of the island but not part of the
formal fit of course).  The errors are produced using an estimate of the
true (signal free) rms R.  It is taken from the image header keyword
ACTNOISE, if present and positive, or by fitting the histogram of the
image with a Gaussian.  Note that this fit excludes both blanked and
pure zero pixels, thereby allowing blanked pixels to be REMAGed to
zero when needed.

     One method - NO LONGER USED - employs John Ball's article in
Method in Computational Physics, 14, 177, 1975.  That gives the noise
in the peak of the fit Gaussian as:
        Delta(P) = R * (Clean-beam area) / (fit area)
the noise in the widths as:
        Delta(W) = Delta(P) / P * W
the noise in the center positions as:
        Delta(X) = Delta(W) / 2
The noise in the total flux is found by the usual combinatorial formulae
and the noise in the position angle is estimated from Condon's article
        Delta(PA) = Delta(p) / P * sqrt (2 WxWn / (WxWx - WnWn))
where Wx is the maximum width and Wn the minimum width.

     The estimates from Ball have been replace by more complicated
formulae from Condon.  The following remarks now apply.  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 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, SAD 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.

ON THE MESSAGE OUTPUT FORMAT
     For readability and compactness the message formats are chosen to
fit in the width of the "printer" (real printer or CRT screen width
given by DOCRT).  The UNITS of the printed flux are scaled so that the
weakest found flux lies in the range 1.0 to 999.0.  The UNIT will be Jy
or milliJy or whatever (indicated in the header to the printout).  If
another source has a flux greater than 999.999 in this UNIT, the flux
will be divided by 1000., before printing, and an asterisk "*" will be
printed before the flux to indicate this.  A flux greater than
999,999.999 will be divided by 1.e6 and be prefixed by a "#".

     If non-point sources are fit (DOWIDTH > 0), and a CLEAN beam size
can be found in the header, the derived source sizes and orientations
are deconvolved, with the results printed in a second listing.  The
numbers printed are the deconvolution of the best fit major axis,
minor axis and position angle.  (arcsec, arcsec, degrees), followed by
the EFACTOR sigma lower limits on these quantities and the EFACTOR
sigma upper limits.  An additional column lists whether the source
source is probably resolved (R), probably unresolved (U), or somehwre
in between (?).

ON THE MEANING OF BMAX, BMIN, AND BPA

     Formally the position angle of anything, such as the major axis of
a gaussian source, is defined as the angle of the thing East from North.
But the direction North only coincides with a column of the pixel array
when you are on a meridian through the phase center.  Thus in general
two different positional angles can be defined, either relative to local
North, or relative to the local Y-direction of the pixel matrix.
Similar small differences can occur in the definition of an astronomical
major axis versus a pixel-based major axis.  Note for example that the
formal astronomical position angle of the point spread function of an
unresolved source changes with position in the field, due to the
curvature of the North/South coordinate system.

Sooooooooo, following the convention in JMFIT, the extention parameters
bpa, bmin, and bmax printed in the message file are the astronomical
values, relative to local North etc., while the parameters put into the
CC file are the pixel based values.  This latter so that operations like
UVSUB will work correctly.  This has the consequence, for example, that
for point sources the message file values differ from the values for the
BEAM, whereas the CC file values should agree with the BEAM.  The MF
file contains everything and the kitchen sink, so all parameters are
present there.
