; UVFIT
;---------------------------------------------------------------
;! Fits source models to uv data.
;# TASK UV ANALYSIS CALIBRATION MODELING
;-----------------------------------------------------------------------
;;  Copyright (C) 1995-1997, 2000, 2003-2004, 2006, 2010-2012, 2016,
;;  Copyright (C) 2018
;;  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
;-----------------------------------------------------------------------
UVFIT     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
UVFIT     Task to fit gaussian or spheres to UV data.
INNAME                             Input UV file name (name)
INCLASS                            Input UV file name (class)
INSEQ             0.0     9999.0   Input UV file name (seq. #)
INDISK            0.0        9.0   Input UV file disk unit #
BCHAN             0.0     2048.0   Lowest channel number 0=>all
ECHAN             0.0     2048.0   Highest channel number
BIF               0.0      100.0   Lowest IF number 0=>all
EIF               0.0      100.0   Highest IF number 0=>all
STOKES                             Stokes type to fit.
SELBAND                            Bandwidth to select (kHz)
SELFREQ                            Frequency to select (MHz)
FREQID                             Freq. ID to select.
SOURCES                            Source name (1 only)
ANTENNAS                           Antennas to include. 0=all
SUBARRAY          0.0     1000.0   Subarray, 0 => 1.
TIMERANG                           Time range to use.
DOCALIB          -1.0      101.0   > 0 calibrate data & weights
                                   > 99 do NOT calibrate weights
GAINUSE                            CL/SN table to apply.
DOPOL            -1.0       10.0   If >0 correct polarization.
PDVER                              PD table to apply (DOPOL>0)
BLVER                              BL table to apply.
FLAGVER                            Flag table version
DOBAND           -1.0       10.0   If >0 apply bandpass cal.
                                   Method used depends on value
                                   of DOBAND (see HELP file).
BPVER                              Bandpass table version
SMOOTH                             Smoothing function. See
                                   HELP SMOOTH for details.
UVRANGE                            Range of baseline lengths
                                   in kilo wavelengths.
SOLMODE                            Soln. mode: 'A&P','A','RE',
                                   ' '
OPCODE                             'GAUS' or 'SPHE'
NGAUSS          0.0         4.0    Number of components
                                   Guess of model parameters
GMAX                                Peak of component (JY)
GPOS                                (X,Y) position (asec)
GWIDTH                              (BMAJ, BMIN, PA) of comp.
                                    (asec, asec, deg)
                                    0->Use small value.
GAINERR                            Antenna gains 0=>1
NITER           0.0     10000.0    Maximum # of iterations
EDROP           0.0                Tolerance level
                                   Solve for model parameters?
DOMAX          -1.0         1.0      Solve for GMAX?
DOPOS          -1.0         1.0      Solve for GPOS?
DOWIDTH        -2.1         1.0      Solve for GWIDTH?
INLIST                             List of sources up to 9999
IMSIZE          0.0                Min, max allowed component
                                   size (asec.)  IMPORTANT
PRTLEV          0.0         2.0    Print level,0=>none, 1,2 some
DOCAT          -1.0         3.0    0=>NO CC file. 1=>write CC.
INVER                              Output CC file version no.
                                   0=>create new
FITOUT                             Not ' ' => write results to
                                   this text file
BADDISK                            Disks to avoid for scratch
----------------------------------------------------------------
UVFIT
Task:  This task will fit a model consisting of elliptical Gaussians
       or optically thin, uniform spheres to a uv data file.  Output
       can optionally be written into a CC (components) table
       associated with the input uv data base.  Multiple channels or
       IFs can be averaged before fitting.  A circular Gaussian can be
       used by using an appropriate value for DOWIDTH.
          The current version will handle up to 2500000 visibilities.
       If there is more data in the file it should be averaged before
       running UVFIT.  There is also a limit of 390 parameters to be
       solved for (60 components, 30 antennas).  Errors are the formal
       errors from the fitting process and are usually unrealistic.

       The task UVGIT is a version of UVFIT that uses a different
       fitting method (that of IMFIT) rather than the one used by
       JMFIT.  UVGIT sometimes converges to correct answers where
       UVFIT does not (and can simply return the initial guess).
Adverbs:
  INNAME.....Input UV file name (name).      Standard defaults.
  INCLASS....Input UV file name (class).     Standard defaults.
  INSEQ......Input UV file name (seq. #).    0 => highest.
  INDISK.....Disk drive # of input UV file.  0 => any.
  BCHAN......First channel to use. Data from BCHAN to ECHAN are
             averaged, 0=>all.
  ECHAN......Highest channel to use. 0=>all higher than BCHAN
  BIF........First IF to use. Data from IF BIF to EIF are
             averaged 0=>all.
  EIF........Highest IF to use. 0=>all higher than BIF
  STOKES.....The desired Stokes type of the data:
             'I','V','Q','U','RR','LL','RL','LR'
  SELBAND....Bandwidth of data to be selected. If more than
	     one IF is present SELBAND is the width of the
	     first IF required. Units = kHz, 0=> all
  SELFREQ....Frequency of data to be selected. If more than
	     one IF is present SELFREQ is the frequency of the
	     first IF required. Units = MHz, 0=> all
  FREQID.....Frequency identifier to select (you may determine
	     which is applicable from the OPTYPE='SCAN' listing
	     produced by LISTR. If either SELBAND or SELFREQ
	     are set their values overide that of FREQID,
	     however setting SELBAND and SELFREQ may occasionally
	     result in an ambiguity, in which case the task will
	     request that you use FREQID.
  SOURCES....Selected source from multi source file, only 1
             allowed.
  ANTENNAS...A list of the antennas to include in the solutions
             If any number is negative then all antennas listed
             are NOT to be used to determine
             solutions and all others are. All 0 => use all.
  SUBARRAY...Subarray number to use. 0 => 1.
  TIMERANG...Time range of the data to be used. In order:
             Start day, hour, min. sec,
             end day, hour, min. sec. Days relative to ref.
             date.
  DOCALIB....If true (>0), calibrate the data using information in the
             specified Cal (CL) table for multi-source or SN table for
             single-source data.  Also calibrate the weights unless
             DOCALIB > 99 (use this for old non-physical weights).
  GAINUSE....(multisource) version number of the CL table to
             apply to the data.   SN version GAINUSE will be
             applied for a single source file.  0 => highest.
  DOPOL......If > 0 then correct data for instrumental polarization as
             represented in the AN or PD table.  This correction is
             only useful if PCAL has been run or feed polarization
             parameters have been otherwise obtained.  See HELP DOPOL
             for available correction modes: 1 is normal, 2 and 3 are
             for VLBI.  1-3 use a PD table if available; 6, 7, 8 are
             the same but use the AN (continuum solution) even if a PD
             table is present.
  PDVER......PD table to apply if PCAL was run with SPECTRAL true and
             0 < DOPOL < 6.  <= 0 => highest.
  BLVER......Version number of the baseline based calibration
             (BL) table to appply. <0 => apply no BL table,
             0 => highest.
  FLAGVER....Specifies the version of the flagging table to be
             applied.  0 => highest numbered table.  <0 => no flagging
             to be applied.
  DOBAND.....If true (>0) then correct the data for the shape of the
             antenna bandpasses using the BP table specified by BPVER.
             The correction has five modes:
             (a) if DOBAND=1 all entries for an antenna in the table
             are averaged together before correcting the data.
             (b) if DOBAND=2 the entry nearest in time (including
             solution weights) is used to correct the data.
             (c) if DOBAND=3 the table entries are interpolated in
             time (using solution weights) and the data are then
             corrected.
             (d) if DOBAND=4 the entry nearest in time (ignoring
             solution weights) is used to correct the data.
             (e) if DOBAND=5 the table entries are interpolated in
             time (ignoring solution weights) and the data are then
             corrected.
  BPVER......Specifies the version of the BP table to be applied.
             <0 => no bandpass correction done.
  SMOOTH.....Specifies the type of spectral smoothing to be applied to
             a uv database . The default is not to apply any smoothing.
             The elements of SMOOTH are as follows:
             SMOOTH(1) = type of smoothing to apply: 0 => no smoothing
               To smooth before applying bandpass calibration
                 1 => Hanning, 2 => Gaussian, 3 => Boxcar, 4 => Sinc
               To smooth after applying bandpass calibration
                 5 => Hanning, 6 => Gaussian, 7 => Boxcar, 8 => Sinc
             SMOOTH(2) = the "diameter" of the function, i.e. width
               between first nulls of Hanning triangle and sinc
               function, FWHM of Gaussian, width of Boxcar. Defaults
               (if < 0.1) are 4, 2, 2 and 3 channels for SMOOTH(1) =
               1 - 4 and 5 - 8, resp.
             SMOOTH(3) = the diameter over which the convolving
               function has value - in channels.  Defaults: 1,3,1,4
               times SMOOTH(2) used when input SMOOTH(3) < net
               SMOOTH(2).
  UVRANGE....Range of allowed baseline lengths in
             kilowavelengths.
  SOLMODE....Solution mode:
              'A&P ' => amplitude and phase.
              'A   ' => Amplitude only
              'RE  ' => Real part only
              '    ' => 'A&P '
  OPCODE.....The model type. 'GAUS' = gaussian, 'SPHE'=>
             Optically thin homogenous sphere.
  NGAUSS......The number of components to use in the fitting.  0->1.
              Maximum number is four using adverbs.  Ignored if INLIST
              is not blank.
  GMAX........The peak value for each component. MUST be given.
              Ignored if INLIST is not blank.
  GPOS........The position (X,Y) for each component.  The values are
              in arcsec. No default.
              Ignored if INLIST is not blank.
  GWIDTH......The major axis, minor axis and position of major axis
              for each component.  The values are in asec with no
              default. For uniform spheres only the first width
              parameter is used which is the radius in asec, the other
              values should be zeroed.  Ignored if INLIST is not blank.
  GAINERR.....The antenna gains to use for each antenna.
              a positive value means fit that parameter with
              the given initial guess.  Zero means use 1.0 and
              don't fit.  Negative values mean use the absolute
              value as the antenna gain but don't fit.
  NITER.......The maximum number of iterations to use in the
              fitting. 0-> NGAUSS * 100.
  EDROP.......A tolerance level for cutoff.  0->10**(-12)
              A value >1 relaxes the convergence criterum,
              A value <1 (but >0) increases.  One or more
              factors of ten may be appropriate.
           Flags for holding any of the parameters constant
            1 -> let parameter vary.  (> 0 -> 1)
           -1 -> keep parameter fixed.  (<= 0 -> -1)
  DOMAX.......Flags for GMAX
              Ignored if INLIST is not blank.
  DOPOS.......Flags for GPOS.  If SOLMODE='A' is used then the
              position of at least one component MUST be fixed.
              Ignored if INLIST is not blank.
  DOWIDTH.....Flags for GWIDTH, if DOWIDTH(2,n) = -2 then use a
              circular gaussian of width GWIDTH(1,n).  For uniform
              spheres ONLY DOWIDTH(1,n) should be set true (+1).
              Ignored if INLIST is not blank.
  INLIST......Set this blank unless you want an input components list.
              Text file containing one line per source, giving
                channel, gmax(i), gpos(1,i), gpos(2,i), gwidth(1,i),
                gwidth(2,i), gwidth(3,i), domax(i), dopos(1,i),
                dopos(2,i), dowidth(1,i), dowidth(2,i), dowidth(3,i)
              blank separated free format and trailing zeros may be
              omitted.   A comment field may follow dowidth(3,i) which
              will be appended to the output line in the FITOUT file.
              Input lines in INLIST are limited to 512 characters and
              the comment field may be no longer than 256 characters.
              If the INLIST adverb is not blank, this file is used to
              set the NGAUSS, initial guess parameters, and whether
              they are fit or held fixed.  There is a limit of 60
              components for any one channel and lines in the file
              beginning with # are comments and those with channel not
              matching BCHAN are ignored.
  IMSIZE......Minumum and maximum size allowed for components
              in arc seconds.  MUST be given.
  PRTLEV......Print level for messages from fitting routine:
              0 => none, 1 => some, 2 => lots.
              (mostly inscrutable stuff)
  DOCAT.......0-> Don't write a CC extension file with the
              results. 1-> Write the derived gaussian
              components to a CC file.  The source name is added
              as a table keyword; use PRTAB to examine.
  INVER.......Output CC file version number desired.
              0=>Create a new extension file.
  FITOUT......If not blank, the name of a text file to receive a
              summary of the fit results.  A one line title plus one
              line per component are printed.  The text is appended to
              pre-existing files.  The text includes a line giving
              column titles and one line per component giving in
              order: channel, component number, RA position and error,
              Dec position and error, Peak flux and error, major axis
              and error, minor axis and error, position angle and
              error, postfit RMS, Chi-squared, start time and end time
              of the interval.  Units for time are days, for flux and
              RMS Jy, for position offsets and sizes arc-seconds or
              milli-arcseconds, and for position angles degrees.
  BADDISK....The disk numbers to avoid for scratch files (sorting
             tables mostly).
----------------------------------------------------------------
UVFIT:  Task to fit a model brightness distribution to uv data.
DOCUMENTOR:  Chris Flatters, NRAO

                             PURPOSE

UVFIT fits a model brightness distribution comprising up to
four discrete components to uv data by adjusting the component
parameters to minimize the chi-squared statistic for the ideal
visibilities derived from the model and the observed data.  The
model components may be elliptical Gaussian brightness
distributions or optically-thin spheres according to the
setting of the OPCODE adverb; UVFIT can not handle a mixture of
Gaussian and spherical components.  The optimized model may be
saved as a clean component (CC) file by setting DOCAT and
INVERS.

You may also use UVFIT to derive constant telescope gain
corrections.  This is useful for calibrating VLBI data where
there may be significant differences in the a priori flux
scaling at each antenna.

                         DATA SELECTION

UVFIT will operate on both single-source and multi-source data
files.  If multi-source data is used then a single source must
be selected using the SOURCES adverb.

Multichannel and multi-IF data may be used as input. All
selected channels and IFs will be averaged together.  VLBI data
should therefore be fringe-fitted before running UVFIT and any
bandpass corrections that are available for spectral- line data
should be applied.

UVFIT will not use more than 2500000 visibilities.  If the selection
criteria include more than this number of visibilities then UVFIT will
use the first 2500000 visibilities that meet the criteria.  If there
are more than 2500000 selected visibilities, it is advisable to
average the data so that uv coverage is not lost.

A priori amplitude calibration should, of course, be applied to
the data.  The STOKES adverb should normally be restricted to
'I', 'RR' or 'LL'.

           FITTING A MODEL BRIGHTNESS DISTRIBUTION

UVFIT works by minimizing the chi-squared statistic (the sum of the
squared deviations of the data from the ideal visibilities).  As is
usual case for minimization techniques, UVFIT will find a local
minimum of chi-squared with respect to all of the model parameters
that are allowed to vary.  In order for this to be a global minimum,
UVFIT must be given an initial guess for the model parameters that is
close to the location of the global minimum.

You supply the initial guess by giving the number of parameters
that you want to use and their parameters.

Gaussian components are specified by giving a peak brightness
(GMAX), a position in arcseconds relative to the phase center
(GPOS), a major axis size in arcseconds, a minor axis size and
a position angle for the major axis in degrees (GWIDTH) for
each component.  If the major axis or minor axis size is given
as zero a small number will be used instead.  Note that
component widths are constrained to be between IMSIZE(1) and
IMSIZE(2) during the fit so you should make sure that IMSIZE(1)
is less than your smallest GWIDTH and that IMSIZE(2) is greater
than larger than your biggest GWIDTH.  Always check IMSIZE
before running UVFIT since its usage in other programs tends to
favour IMSIZE(1) being equal to IMSIZE(2) which will force
UVFIT to use circular Gaussian components all of the same size.

The same parameters are used for spherical components with the
exception that only the major axis size is specified using
GWIDTH and the other two components of GWIDTH are ignored.

You should then choose whether you wish to fit to amplitude and
phase or amplitude only.  If you have VLBI data that has not
yet been self-calibrated in phase them you should use
amplitudes only (SOLMODE='A') otherwise you should normally use
both amplitudes and phases (SOLMODE='A&P').

The next step is to choose which parameters you will allow to
vary during the fit.  Do this by setting the appropriate flags
in DOMAX, DOPOS and DOWIDTH to +1 for the parameters that will
vary and -1 for those that will be held constant.  If you are
using Gaussian components you can force any of the to be
circular by setting DOWIDTH(2, i) to 2 (to constrain the i-th
component).  Note that if you are fitting to amplitudes only
then the position of at least one component may be fixed.

You may then start UVFIT to find the optimized parameters.

I recommend that you run UVFIT with PRTLEV set to 1.  This will
give a small amount of information about every iteration.  If
the iteration number goes over about thirty or if you don't see
the chi-squared statistic (labelled F in the iteration summary)
go down significantly in the first few iterations then UVFIT is
going nowhere and you might as well abort it and start again
with a different model (small numbers of iterations can also be
a problem as discussed below).  PRTLEV 0 doesn't give any
information about the iterations and PRTLEV 2 gives much more
information than is needed.

                   EVALUATING THE RESULTS

UVFIT will report the post-fit chi-squared statistic and rms.
Ideally, the rms should be close to the average thermal noise
for a data point after averaging over the selected channels and
IFs.  If the rms is much higher than this you obviously have a
poor fit to the data.  To use the post-fit chi-square you need
to determine the number of degrees of freedom, d = k * v - p
where k = 1 if you are fitting amplitudes only or k = 2 if you
are fitting amplitude and phase, v is the number of
visibilities used and p is the number of parameters allowed to
vary.  If the reduced chi-squared statistic (obtained by
dividing the post-fit chi-squared by d) is about 1 or less then
you have a good fit; it it is much greater than one then you
have junk.

More quantitatively, you can use the chi-squared statistic to
get the probability of obtaining the observed data given the
model in question (assuming errors have a Gaussian
distribution).  You will probably need a statistics program to
calculate this since d is usually much larger than the number
of degrees of freedom for which the chi-square distribution is
tabulated in reference books.  This may not seem particularly
useful on its own but given two models, A and B, with no reason
to prefer one over the other and given the probabilities of
obtaining the observed data given model A, P(d|A) and that of
obtaining the observed data given model B, P(d|B), the
probabilities of either model being correct given the observed
data are in the ratio P(A|d) / P(B|d) = P(d|A) / P(d|B) (this
is a trivial application of Bayes' theorem).  More generally,
the higher the probability of obtaining the data given a
particular model, the more you should believe the model.

You should note, however, that the calculation of the
chi-squared statistic assumes that the data are weighted by the
inverse of the noise variance; if this is not true then the
value of the chi-squared statistic is useless (although lower
chi-square numbers are still preferred over larger ones).

                 NOTES ON THE ERROR BOUNDS

UVFIT presents nominal error bounds on all of the parameters
that were allowed to vary during the fit.  These errors are
actually the square roots of the diagonal elements of the
covariance matrix (the matrix obtained by inverting the Hessian
matrix -- the matrix where element ij is the partial derivative
of chi-squared with respect to the i-th and j-th model
parameters) for the final model parameters.  Provided that the
data weights reflect the expected noise in the data these
elements give the formal 1-sigma error for each parameter with
all other components at their fitted value.  This means that
they will underestimate the true errors since the error ranges
of the parameters are all interconnected; it is, however
difficult to derive any better estimates.

A further complication arises from the algorithm used by UVFIT.
UVFIT uses a member of the quasi-Newton family of minimization
algorithms.  Quasi-Newton algorithms build up the Hessian
matrix through successive approximations.  This has the
consequence that the error estimates will not be accurate
unless UVFIT goes through several iterations.  If UVFIT
converges in less than about 5 iterations you may want to
change one or more of the parameters (to force it to do more
iterations) and run it again to get better error estimates.

                   FITTING ANTENNA GAINS

You may also introduce antenna gain factors as parameters to be
fit.  If GAINERR(i) is zero the gain factor for antenna i is
fixed at unity (the default); you may fix the gain factor at
another value, C, by setting GAINERR(i) to -C or you may allow
UVFIT to fit the gain factor by setting GAINERR(i) to a
positive value (which is taken as the starting guess).  You can
use the gain factors output by UVFIT to correct antenna gains
using CLCOR (set CLCORPRM(1) to the factor given by UVFIT for
the antenna in question, set the remaining elements of CLCORPRM
to zero and set OPCODE to 'GAIN').  This is useful for VLBI
data where different antennae may have significant systematic
errors in gain measurements.

Unless you know the flux density of the source, in which case
the model fluxes should be fixed by setting all elements of
DOMAX to -1, the gain of at least one antenna in the array
should be fixed otherwise the flux scale will drift.
