; UVMOD
;---------------------------------------------------------------
;! Modify UV database by adding a model incl spectral index
;# TASK UV ANALYSIS MODELING SPECTRAL
;-----------------------------------------------------------------------
;;  Copyright (C) 1995, 1997, 2001, 2005, 2007-2010, 2014-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
;-----------------------------------------------------------------------
UVMOD     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
UVMOD     Task which inserts a model into 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 #
SOURCES                            Source name
QUAL            -10.0              Calibrator qualifier -1=>all
CALCODE                            Calibrator code '    '=>all
STOKES                             Stokes of output
TIMERANG                           Time range to use
SELBAND                            Bandwidth to select (kHz)
SELFREQ                            Frequency to select (MHz)
FREQID                             Freq. ID to select.
SUBARRAY          0.0     1000.0   Sub-array, 0=>all
BIF                                Low IF number to do
EIF                                Highest IF number to do
BCHAN             0.0              First channel included
ECHAN             0.0              last channel included
DOCALIB          -1.0      101.0   > 0 calibrate data & weights
                                   > 99 do NOT calibrate weights
GAINUSE                            CL (or SN) table to apply
DOPOL            -1.0       10.0   If >0.5 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.5 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.
DOACOR                             Include autocorrelations?
OUTNAME                            Output UV file name (name)
OUTCLASS                           Output UV file name (class)
OUTSEQ           -1.0     9999.0   Output UV file name (seq. #)
OUTDISK           0.0        9.0   Output UV file disk unit #
NGAUSS            0.0     9999.0   Number of sources
CTYPE             0.0        5.0   Component type: 0 point, 1
                                   Gaussian, others see help
FMAX                               Peak value of sources
FPOS                               (X,Y) position in asec
FWIDTH                             (Major,Minor,PA) in
                                   asec,asec,deg
INLIST                             List of sources up to 9999
ZEROSP                             Zero spacing flux for (1) I,
                                   (2) Q, (3) U and (4) V in Jy.
FLUX                               Noise level in Jy/Weight or
                                   Jy if FACTOR=0.
DPARM                              Spectral indices to use
FACTOR                             Multiplication factor.
WTUV                               If >0 then weights are
                                   reset to a value of 1.
DOHIST                             > 0 => list sources in
                                   history file
FQCENTER                           >= 0 -> center frequency axis
BADDISK                            Disks to avoid for scratch
----------------------------------------------------------------
UVMOD
Type: Task
Use: Modification of existing UV data by the addition of models.
Adverbs:
  INNAME.....Input image name (name).       Standard defaults.
  INCLASS....Input image name (class).      Standard defaults.
  INSEQ......Input image name (seq. #).     0 => highest.
  INDISK.....Disk drive # of input image.   0 => any.
  SOURCES....Source to be copied.   '  '=> all; if any starts with a
             '-' then all except ANY source named.
  QUAL.......Qualifier of source to be copied. -1 => all.
  CALCODE....Calibrator code of sources to copy. ' '=> all.
  STOKES.....Specifies which STOKES parameters are written in the
             output data set:  ' ' => 'FULL'
               'I','Q','U','V', 'IV', 'IQU', 'IQUV'
               'RR','LL', 'RL', 'LR', 'RRLL', 'RLLR', 'RLRL'
               'XX','YY', 'XY', 'YX', 'XXYY', 'XYYX', 'XYXY'
             'HALF', 'CROS', and 'FULL' have sensible interpretations
             depending on the Stokes present in the data.  The last in
             each of the 3 rows above == 'FULL'.  Note that many
             combinations of polarizations in the input and values
             above are not supported.
  TIMERANG...Time range of the data to be copied. In order: Start day,
             hour, min. sec, end day, hour, min. sec. Days relative to
             ref. date.
  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. For data which contain multiple
             bandwidths/frequencies the task will insist that some form
             of selection be made by frequency or bandwidth.
  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.
  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
             override that of FREQID.  However, setting SELBAND and
             SELFREQ may result in an ambiguity.  In that case, the task
             will request that you use FREQID.
  SUBARRAY...Sub-array number to copy. 0=>all.
  BIF........First IF to include.  0 -> 1.
  EIF........Last IF to include.   0 -> max.
  BCHAN......First channel to copy. 0=>all.
  ECHAN......Highest channel to copy. 0=>all higher than BCHAN
  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....version number of the CL table to apply to multi-source
             files or the SN table for single source files.
             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 apply. <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.
             IMAGR uses DOBAND as the nearest integer; 0.1 is therefore
             "false".
  BPVER......Specifies the version of the BP table to be applied
                0 => highest numbered table.
               <0 => no bandpass correction to be applied.
  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).
  DOACOR.....> 0 => include autocorrelations as well as cross
             correlation data.
  OUTNAME....Output image name (name).      Standard defaults.
  OUTCLASS...Output image name (class).     Standard defaults.
  OUTSEQ.....Output image name (seq. #).    0 => highest unique
  OUTDISK....Disk drive number of output image.  0 =>
             highest number with sufficient space.
  NGAUSS.....Number of sources 1 - 4 use CTYPE, FMAX, FPOS, and
             FWIDTH; > 4, use INLIST file.  Limit 9999.
             Note, if INLIST not blank and NGAUSS > 1, the INLIST file
             will be used BUT the value of NGAUSS will limit how many
             components are read from INLIST.  Set it large when you
             use INLIST if you want the full contents of INLIST.
  CTYPE......Type of each source: 1 Gaussian, 2 solid disk, 3 solid
             rectangle, 4 optically thin sphere, 5 exponential
             otherwise point.
  FMAX.......I polarization max of each source
  FPOS.......Offset value of the X,Y centroid of the sources in
             arcsec. Positive values mean increasing R.A. and DEC.
             These are offsets in the coordinates with R.A. scaled by
             cos(DEC) and are applied with the W term taken into
             account (as of 1997-04-16).  See Explain at the end.
  FWIDTH.....Model major and minor axis in arcsec and PA in degrees.
             Full width or full width to half maximum (GAUS).
  INLIST.....Set this blank unless you want an input components list.
             Text file containing one line per source, giving
                I, DX, DY, Maj, Min, PA, type, spix, Dspix, Q, U, V
             blank separated free format and trailing zeros may be
             omitted.  Used if NGAUSS > 1 and INLIST not blanks.  The
             resulting NGAUS will be min (NGAUSS, #lines in INLIST).
             Limit 9999.  Note that I, Q, U, V and the spectral
             indexes are assumed to apply at the header (reference)
             frequency.  Any line in INLIST beginning with a $ or a #
             is taken as a comment and ignored.  Negative values for I
             are allowed, but I=0 lines are ignored.
  ZEROSP.....Zero spacing flux of (2) Q, (3) U and (4) V in Jy used
             with the first source only and only if INLIST blank or
             NGAUSS=1.  For more polarization in models, use the
             INLIST option.  Be sure to put in the zeros (or other
             values) for the parameters that precede the Q,U,V values
             on each line in INLIST.
  FLUX.......Noise level to be added (in Jy. per Weight or Jy if
             FACTOR=0).  This means, when FACTOR is not zero, that
             the noise added is FLUX/sqrt(weight) so that FLUX scales
             the data rms, assumed to be 1/sqrt(weight).
  DPARM......(1-5) Spectral index to use for components 1,2,3,4,noise
             (6-10) Spectral index curvature for comp 1,2,3,4,noise
             A non-blank INLIST gives these if NGAUS > 1.  DPARM(5)
             and (10) still used for noise.  See explain at the end.
             Curvature numbers are based on base 10 logarithms with a
             reference frequency of 1.0 GHz.
  FACTOR.....Factor by which original data are multiplied before they
             are added to the model.  If FACTOR = 0 then only the model
             will be left.
  WTUV.......If WTUV > 0 then all Weights will be set to a
             value of 1 in the output file.
  DOHIST.....List sources in history file if NGAUS <= 4 and/or DOHIST
             > 0.
  FQCENTER...>  0 => Change frequency axis reference pixel to
                     Nchan / 2 + 1
             else => do not change reference pixel
  BADDISK....The disk numbers to avoid for scratch files (sorting
             tables mostly).
----------------------------------------------------------------
UVMOD: Task which modifies UVDATA by scaling the existing data,
       and adding a specified model (See also IMMOD).
DOCUMENTATION: Eric R. Nelson NRAO/VLA/UNM
DATE OF DOCUMENTATION: 14 JUNE 1983
RELATED PROGRAMS: IMMOD, UVMAP, APCLN, COMB

                  PURPOSE

     UVMOD modifies an already existing UV data file by the addition of
one of several model types.  The original data may be scaled by a
multiplicative factor, including negative values and zero, before they
are added to the model.  Random noise may also be added to the UV data.
The program could be useful in investigating the affects of CLEAN on a
specific geometry, or for removing models from data, i.e. planetary
disks.

     The six available models to choose from are 0) point source, 1)
Gaussian, 2) solid disk  3) solid rectangle, 4) optically thin sphere,
and 5) exponential.  These models are first Fourier transformed and
then added to the UV data.  The resulting functions are:

          0) Point -> A constant visibility amplitude is added to the
                      data.  The GWIDTH adverbs have no affect on this
                      model.  Used with CTYPE <= 0 or > 5.

          1) Gaussian -> The function EXP(-3.559707*R**2) is added to
                      the UV data.  The function R is given by:
                         R = Sqrt(UU**2 + VV**2) where
                         UU = BMAJ*(V*COS(BPA)+U*SIN(BPA))
                         VV = BMIN*(U*COS(BPA)-V*SIN(BPA))

          2) Disk ->  The function J1(R)/R is added to the UV data,
                      where J1 is the Bessel function of order 1 and R
                      is the same as above.

          3) Rectangle -> The function SINC(UU)*SINC(VV) is added to
                      the UV data, where SINC(X) = SIN(X)/X.  UU and
                      VV are defined as above for the Gaussian.

          4) Sphere -> The function (SIN(A)/A - COS(A)) / (A*A) is
                      added to the UV data where
                         A = BMAJ * Sqrt (U*U + V*V)
                         A = max (A, 2 pi / 100)
                      The GWIDTH adverbs have no affect on this model.

          5) Exponential -> The function
                           2 Pi / (1 + a * a * R * R) ** 3/2
                      is added to the UV data where R is defined in 2
                      above and a is Pi/ln(2).

Note that all functions are scaled by the total flux and a complex
vector representing the phase of the model before being added to the
scaled input visibility data.  The spectral index is applied to make
the peak flux
         log(F/F0) = spix * log(nu/nu0) + Dspix * log^2 (nu/nu0)
where F0 is the model peak flux at the header reference frequency
nu0.  log functions are base 10.


BMAJ, BMIN, BPA :

     In this discussion, BMAJ = GWIDTH(1,i) for the i'th component,
BMIN = GWIDTH(2,i) for the i'th component, and BPA = GWIDTH(3,i) for
the i'th component,  The dimensions of the resulting functions are
determined by BMAJ, BMIN and BPA (position angle).  For the Gaussian
the first two values are the FWHM of the two axis.  For the Disk and
the Rectangle, the first two values are the absolute dimensions of the
two available axis.  If BMAJ and BMIN are both zero then all the
models reduce to the point model.

ZEROSP:

     This adverb, array indices 2-4, allows you to add polarization to
the first model component only.

FACTOR :

     The FACTOR term allows one to add a scaled version of the
original data to the model.  FACTOR is simply multiplied by the
original data which is then added to the model.  If FACTOR = 0, then
only the model will remain in the final UV data base.

Coordinate considerations:

     FPOS is translated to an RA, Dec following the formula (assuming
no rotations):
          Dec = Dec0 + FPOS(2)
          RA  = RA0 + FPOS(1) / cos (Dec0)
These are then turned into l,m,n for phase = ul + vm + wn as (for -SIN
geometry):
          l = cos (Dec) * sin (Ra-Ra0)
          m = sin (Dec) * cos (Dec0) -
              cos (Dec) * sin (Dec0) * cos (Ra-Ra0)
          n = sin (Dec) * sin (Dec0) +
              cos (Dec) * cos (Dec0) * cos (Ra-Ra0)
Suitable formulae are used for -NCP geometry as well.

Spectral index parameters are entered as:
          x = log10 (freq / 1 GHz)
          F = F0 * exp ((Spix + Dspix * x) * x)
where F is the flux and F0 is the flux at 1.0 GHz.  The model fluxes
at the reference frequency are converted to fluxes at 1 GHz.
