; UVSUB
;---------------------------------------------------------------
;! Subtracts/divides a model from/into a uv data base
;# TASK UV IMAGING MODELING CALIBRATION
;-----------------------------------------------------------------------
;;  Copyright (C) 1995-1997, 1999, 2001-2003, 2007-2009, 2014, 2019
;;  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
;-----------------------------------------------------------------------
UVSUB     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
UVSUB     Task to subtract CLEAN components from a uv data base.
                                   (or, to divide observed visi-
                                   bility by model visibility)
INNAME                             Input UV file name (name)
INCLASS                            Input UV file name (class)
INSEQ             0.0     9999.0   Input UV file name (seq. #)
INDISK                             Input UV file disk unit #
NMAPS             0.0     4096.0   No. maps to use for model.
CHANNEL          -1.0     9999.0   Spectral channel (0=>all)
                                   Use 0 for continuum
BIF               0.0     9999.0   First IF (0=>1)
EIF               0.0     9999.0   Highest IF (0=>BIF to last)
IN2NAME                            Cleaned map name (name)
IN2CLASS                           Cleaned map name (class)
IN2SEQ            0.0     9999.0   Cleaned map name (seq. #)
IN2DISK                            Cleaned map disk unit #
INVERS           -1.0    46655.0   CC file version #.
OUTNAME                            Output UV file name (name)
OUTCLASS                           Output UV file name (class)
OUTSEQ           -1.0     9999.0   Output UV file name (seq. #)
OUTDISK                            Output UV file disk unit #.
BCOMP                              First CLEAN comp to sub.
                                   1 per field.
NCOMP                              Last CLEAN comp to sub.
                                   to use (0 => all)
FLUX                               Lowest CC component used.
CMETHOD                            Modeling method:
                                   'DFT','GRID','    '
CMODEL                             Model type: 'COMP','IMAG'
                                   'SUBI' (see HELP re images)
FACTOR                             Factor times CLEAN fluxes.
                                   0->1.0  Subtract
                                     -1.0  Add
OPCODE                             'DIV ' => divide visibility
                                     observation by model vis.
                                   'DIV4' => divide visibility
                                     observation by I model vis.
                                     All 4 Stokes divided
                                   'MODL' => replace visibility
                                     with model visibility
                                   anything else => subtract
SMODEL                             Source model, 1=flux,2=x,3=y
                                   See HELP SMODEL for models.
BADDISK                            Disks to avoid for scratch
----------------------------------------------------------------
UVSUB
Task:  Subtracts/divides a model from/into a uv data base.  The model
       may be a specific model, a set of CLEAN components files, or a
       set of images.  "CLEAN" models may be points, Gaussians or
       uniform, optically thin spheres.  The task will also subtract the
       model and then re-read the result replacing the difference with
       the model value.
            Model images made with both values of IMAGR's DO3DIMAG
       option are handled correctly, as are multi-scale images.  Set
       NMAPS = NFIELD * NGAUSS.
            UVSUB works only on single-source files.

       NOTE: this task does NOT apply flagging or calibration tables
       to the input UV data.  Run SPLIT first if that operation is
       desired.
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.
  NMAPS......Number of image files to use for model.  For multi-scale
             models, set NMAPS = NFIELD * NGAUSS to include the Clean
             components of the extended resolutions.  If more than one
             file is to be used, the NAME, CLASS, DISK and SEQ of the
             subsequent image files will be the same as the first file
             except that the LAST 3 or 4 characters of the CLASS will
             be an increasing sequence above that in IN2CLASS.  Thus,
             if INCLASS='ICL005', classes 'ICL005' through 'ICLnnn'
             or 'ICnnnn', where nnn = 5 + NMAPS - 1 will be used.  Old
             names (in which the 4'th character is not a number) are
             also supported: the last two characters are '01' through
             'E7' for fields 2 through 512.  In old names, the highest
             field number allowed is 512; in new names it is 4096.
  CHANNEL....Frequency channel, 0 => all (use 0 for continuum)
             If > 0, then subtract/divide only this channel but copy
             all channels with the others unchanged.  NOTE WELL:
             EIF will be set to BIF if CHANNEL >= 1.
  BIF........First IF to process. 0=>1
  EIF........Highest IF to process 0=> do BIF to highest.
             Note: not all data sets will have IFs.  See note under
             CHANNEL directly above too.  If the IF axis precedes the
             FREQ axis in the header, then all IFs must be done if all
             spectral channels are to be done.
  IN2NAME....Model map name (name).      Standard defaults.
  IN2CLASS...Model map name (class).     Standard defaults.
  IN2SEQ.....Model map name (seq. #).    0 => highest.
  IN2DISK....Disk drive # of model map.  0 => any.
  INVER......CC file ver. number.          0 => highest.
  OUTNAME....Output UV file name (name).   Standard defaults.
  OUTCLASS...Output UV file name (class).  Standard defaults.
  OUTSEQ.....Output UV file name (seq. #). 0 => highest unique.
  OUTDISK....Disk drive # of output UV file.  0 => highest with
             space
  BCOMP......The first clean component to process. One value is
             specified for each field used.
  NCOMP......Number of Clean components to use for the model, one
             value per field.  If all values are zero, then all
             components in all fields are used.  If any value is not
             zero, then abs(NCOMP(i)) (or fewer depending on FLUX and
             negativity) components are used for field i, even if
             NCOMP(i) is zero.  If any of the NCOMP is less than 0,
             then components are only used in each field i up to
             abs(NCOMP(i)), FLUX, or the first negative whichever
             comes first.  If abs(NCOMP(i)) is greater than the number
             of components in field i, the actual number is used.  For
             example
                   NCOMP = -1,0
             says to use one component from field one unless it is
             negative or < FLUX and no components from any other
             field.  This would usually not be desirable.
                   NCOMP = -1000000
             says to use all components from each field up to the
             first negative in that field.
                   NCOMP = -200 100 23 0 300 5
             says to use no more than 200 components from field 1, 100
             from field 2, 23 from field 3, 300 from field 5, 5 from
             field 6 and none from any other field.  Fewer are used if
             a negative is encountered or the components go below
             FLUX.
  FLUX.......Only components > FLUX in absolute value are used in the
             model.
  CMETHOD....This determines the method used to compute the
             model visibility values.
             'DFT' uses the direct Fourier transform, this
             method is the most accurate.
             'GRID' does a gridded-FFT interpolation model
             computation.
             '    ' allows the program to use the fastest
             method.
             NOTE: data in any sort order may be used by the
             'DFT' method but only 'XY' sorted data may be used
             by the 'GRID' method.
             NOTE: CMETHOD='GRID' does not work correctly for RL
             and LR data; DO NOT USE CMETHOD='GRID' for RL, LR!
  CMODEL.....This indicates the type of input model; 'COMP' means that
             the input model consists of Clean components, 'IMAG'
             indicates that the input model consists of images.
             'SUBI' means that the model consists of a sub-image of
             the original IMAGR output.  If CMODEL is '   ' Clean
             components will be used if present and the image if not.
             SUBI should work for sub-images made with DO3DIM true and
             sib-images of the central facet made with DO3DIM false,
             but probably will not work well for shifted facets with
             DO3DIM false.  Use BLANK rather than SUBIM in such cases.
             Note also that Clean images do not make good models.  The
             Clean components have been convolved with the Gaussian
             Clean beam making their Fourier transform be rather
             tapered compared to the original uv data.  CCRES can make
             a points image for you.
  FACTOR.....This value will be multiplied times the CLEAN component
             flux densities before subtraction.  The default 0->1.0, so
             the clean component model will be subtracted from the UV
             data.  FACTOR=-1 will add the clean component model to the
             UV data.  FACTOR will be adjusted by UVSUB for image
             models that are in JY/BEAM.  FACTOR is used with all
             OPCODEs.
  OPCODE.....OPCODE='DIV ' => divide observed visibility by model
                visibility.
             OPCODE='DIV4' => divide observed visibility (all 4
                Stokes) by Ipol model visibility.
             OPCODE='MODL' => replace the visibility with the model
                visibility.
             Any other setting of OPCODE causes the task to subtract the
             model visibility from the observed visibility (the normal
             mode of operation).
  SMODEL.....A single component model to be used instead of a CLEAN
             components model; if abs (SMODEL) > 0 then use of this
             model is requested.
                SMODEL(1) = flux density (Jy)
                SMODEL(2) = X offset in sky (arcsec)
                SMODEL(3) = Y offset in sky (arcsec)
                SMODEL(4) = Model type:
                  0 => point model
                  1 => elliptical Gaussian and
                       SMODEL(5) = major axis size (arcsec)
                       SMODEL(6) = minor axis size (arcsec)
                       SMODEL(7) = P. A. of major axis (degrees)
                  2 => uniform sphere and
                       SMODEL(5) = radius (arcsec)
  BADDISK....The disk numbers to avoid for scratch files.
----------------------------------------------------------------
UVSUB: task to subtract a clean model from a uv data base
DOCUMENTOR: B.E.Turner  NRAO/Charlottesville
RELATED PROGRAMS: APCLN, PRTCC, UVMAP, UVPLT, CLIP, UVFND

                        PURPOSE

     The task UVSUB can be used to subtract (or add) or divide
the Fourier transform of a specified model file(s) from/to/into
a visibility data set.  The input model may consist of the
CLEAN components associated with one of more input files;
one or more input images; or, a specified point model.
In addition, there are two model methods. The DFT method
(CLEAN components and point model only) takes a direct
Fourier transform to determine the model visibility values.
The gridded-FFT method does an FFT on either an image
consisting of the CLEAN components or another specified image
and interpolates the model visibility values. There are at
least three principal circumstances for which this is required.
(The division option can be used to solve for "correlator
offsets." See VLA Scientific Memo. No. 152.)

SORT ORDER:
     Data in any sort order may be used with clean components
using the 'DFT' model computation method.  All other cases
require 'XY' sorted data.


EDITING A VISIBILITY DATA SET FOR BAD DATA:
     Suppose the field contains a moderately strong source
source, which produces signal on many baselines comparable or
stronger than the suspected bad data.  One then maps and cleans
the source (down to a level which does not include the effects
of the bad data), then subtracts the Fourier transform of the CC
(Clean Components )file with UVSUB, leaving a vis data set in
which the bad data should now be dominant and easily isolated.
One then can identify it explicitly by use of UVFND, and flag
it.  Or, one can judge its level (above noise) by UVPLT, then
use CLIP to flag all vis data above a desired level.  Finally,
one adds the Fourier transform of the CC file back into the
edited vis data set, using UVSUB with FACTOR = -1.0, to produce
an edited final data set.

SUBTRACTING "CONTAMINATING" SOURCES FROM A FIELD:
(Note: the following can also be done using IMAGR.)
     Often, one or more strong sources may limit the dynamic
range of a radio map.  If the sources lie within the inner
quarter of the radio map, the source and its sidelobes can often
be removed by cleaning using APCLN and instrumental errors can
be removed by self-calibration using CALIB.
     If these sources are sufficiently far from the phase center
then two problems occur: 1) unreasonably large maps would be
needed to place the sources within the inner quarter of the map
where they could be adequately cleaned and 2) distortions are
caused by the sky curvature and the non-coplanar baselines
inherent to VLA synthesis.
     These sources and their effects can be removed making
relatively small maps centered around each bright source using
UVMAP with the UVSHIFT option.  This task handles the phase
shifts accurately so that the non-coplanar distortions are
removed.  The source is then cleaned and the subsequent CC file
is then used with UVSUB.  UVSUB accurately handles the
appropriate phase shifts.  This subtraction can be done, in
turn, on all of the bright sources and several iterations may be
required if the sidelobes interact among the sources.
Generally, if 95% of the strong sources are removed, a large
gain is made in lessening these distortions.
     Other distortions associated with bright sources far from
the phase center of the observations (bandwidth smearing,
time-averaging smearing, inaccurate pointing of the antennas)
are not corrected by this technique; only problems associated
with the non-coplanar baseline effect and inaccuracies of
cleaning a source which is not near the middle of a map.

SUBTRACTION OF CONTINUUM FROM SPECTRAL LINE DATA
     The continuum can be subtracted from spectral line channels
either in map-plane or in u-v data.  If the source size S, the
synthesized beam width s, the center frequency F, and the
frequency range f satisfy the following constraint

              (S/s) * (f/F) < 0.1

then little scaling occurs from channel to channel and the dirty
continuum map can be subtracted from the channel maps.  The
continuum map is generally produced from a sum of the channel
maps which have no significant line radiation or absorption.
Cleaning and self-calibration of the differential channel maps
may be useful if there is sufficient signal to noise.  For the
purposes of self-calibration, the sum of the u-v data of the
channels with significant line emission or the original
continuum u-v data set can be used to increase the dynamic
range.  The self-calibration solution can then be applied to
each channel.
     If the above constraint is not satisfied, then it is more
accurate to subtract the continuum model directly from the u-v
data for each channel.  The continuum model is obtained by
summing the channel maps which contain little significant line
emission and then cleaning the map to obtain the CC file.  UVSUB
is then used to subtract this continuum model from each of the
continuum u-v data sets.  As above, cleaning and
self-calibration may be needed for the residual channel data is
the signal to noise warrents it.
     Because of the differential scaling between channels which
may be summed to generate a continuum map, several continuum
maps, each composed of relatively closely-spaced frequency
channels, may have to be made and the cleaned separately.  Each
CC file will then have to be subtracted from the channels using
UVSUB.  Factor will then have to be adjusted so that the sum of
the FACTOR's used is 1.0 with each FACTOR akin to a weight of
that particular continuum map (i.e. the number of channels).

MULTIPLE IF DATA SETS

     An IF consists of a set of one or more equally spaced
frequency channels; multiple IFs with arbitrary frequency
spacings are allowed.  For data sets with more than one IF the
range of IF to be processed can be specified by BIF and EIF.
If multiple IF are to be processed, all frequency channels in
each IF must be processed.


CALIBRATION:
     UVSUB can be used in conjunction with DBCON, UVSRT, CALIB,
and IMAGR to calibrate uv data.  Each calibrator source is
processed separately using CALIB etc. until an acceptable model
is derived.  Then the model is divided into the RAW calibrator
data.  DBCON can be used to concatinate the divided data from
all calibrators using the options DOPOS=FALSE, and DOARRAY=TRUE.
At this point the data are all in the form of the response of a
1 Jy point source at the phase center.  After sorting this
combined data set to 'TB' order (if necessary) run CALIB with
APARM(2) = 1 which will cause it to skip the model subtraction
step and produce a gain table which, when applied by most tasks
to calibrate the data for a program source.


NMAPS:

     Multiple input files may be specified although all must be
of the same type (i.e. CLean components files or images).
Subsequent model files must have the same NAME, CLASS, DISK and
SEQ as the first file (IN2NAME, IN2CLASS, IN2DISK, IN2SEQ)
except that the LAST three characters of the CLASS should be '001'
thru '999' for the second through the 999 file.  This is the
convention used by IMAGR to name multiple fields.

CMETHOD:

     There are two model computation methods available. The
'DFT' method is more accurate but under some circumstances the
'GRID' method is considerable faster. Only the 'GRID' method
may be used for image models.

TIMING:
     UVSUB on a VAX 11-780 with an array processor runs slowly
With an otherwise empty VAX, the run time for CLEAN component
is:

DFT method
  T(real) = 6.0E-6 * No. vis * No. CLEAN components   seconds

Gridded-FFT method
  T(real) = No. Fields * (No. frequencies) *
           ((No. vis * 0.4E-3)
            + SQRT (4*NX*NY) * 4.0E-2
            + No. Clean components * 1.0E-3   seconds)

where the size of the image is NX x NY.  The time for the model
calculations using an image model is similar the that for the
Gridded-FFT CLEAN components method for 0 clean components with
the addition of the time for the FFT which is approximately:
       No. fields * No. frequencies * (SQRT(NX*NY)**1.3) * 180

