; SPCOR
;---------------------------------------------------------------
;! Task to correct an image for a primary beam and spectral index
;# TASK IMAGING ANALYSIS
;-----------------------------------------------------------------------
;;  Copyright (C) 2011, 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
;-----------------------------------------------------------------------
SPCOR     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SPCOR:  Applies primary beam & spectral index corrections
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 #
IN2NAME                            Spectral index image name
IN2CLASS                           Spectral index image class
IN2SEQ          0.0      9999.0    Spectral index image seq. #
IN2DISK         0.0         9.0    Spectral index disk drive #
IN3NAME                            Spectral curvature image name
IN3CLASS                           Spectral curvature class
IN3SEQ          0.0      9999.0    Spectral curvature image seq
IN3DISK         0.0         9.0    Spectral curvature disk drive
BLC             0.0      4096.0    Bottom Left Corner
TRC             0.0      4096.0    Top Right Corner
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 disk drive #
PBPARM                             Beam parameters:
                                   (1) Cutoff 0 -> 0.023
                                       >= 1 no primary beam
                                   (2) > 0 -> Use (3)-(7)
                                   (3)-(7) Beam shape
X                                  Frequency to reference GHz
                                   0 -> use image reference
COORDINA                           Pointing position RA DEC as
                                      HH MM SS.SS, DD MM SS.S
                                   ** used if not all 0 **
RADIUS           0.7               radius over which to average
                                   in spectral index image(s)
DOBLANK         -1.0        1.0    > 0 => blank when the SPIX
                                   image is blanked, else just
                                   do no correction
DOINVERS        -1.0        1.0    > 0 => multiply image by beam
                                   and undo the spectral index
                                   correction
----------------------------------------------------------------
SPCOR
Type: Task
Use:  SPCOR will apply a primary beam correction (divide by the primary
      beam gain factor) to images and apply spectral index corrections
      based on images of spectral index and spectral index curvature.
      It can also do the inverse of the corrections.  The image is
      required to have RA, DEC and FREQ (or FQID) axes in order to
      compute the correction.  However, they may be in any order and
      "cubes" are fully supported.  Note that the spectral index
      images do not have to be on the same geometry as the input
      image.  The output image is corrected to be as if it were
      observed at the frequency given by X - thus the image at each
      point is multipled by Cor(RA,Dec,Freq)/Cor(RA,Dec,X).
Adverbs:
  INNAME......Input name of image(name).     Standard defaults.
  INCLASS.....Input name of image(class).    Standard defaults.
  INSEQ.......Input name of image(seq. #).   0 => highest.
  INDISK......Disk drive # of image.         0 => any.
  IN2NAME.....Spectral index name of image(name).     ' ' => none
  IN2CLASS....Spectral index name of image(class).    ' ' => none
  IN2SEQ......Spectral index name of image(seq. #).   0 => highest.
  IN2DISK.....Spectral index disk drive # of image.   0 => any.
  IN3NAME.....Spectral curvature name of image(name).     ' ' => none
  IN3CLASS....Spectral curvature name of image(class).    ' ' => none
  IN3SEQ......Spectral curvature name of image(seq. #).   0 => highest.
  IN3DISK.....Spectral curvature disk drive # of image.   0 => any.
  BLC.........The bottom left-hand pixel of the input image which
              becomes the bottom left corner of the corrected
              subimage.  The value (0,0) means (1,1).
  TRC.........The top right-hand pixel of the input image which
              becomes the top right corner of the corrected subimage.
              The value (0,0) means take the top right hand corner of
              the image.
  OUTNAME.....Output name of image(name).   Standard defaults.
  OUTCLASS....Output name of image(class).  Standard defaults.
  OUTSEQ......Output name of image(seq. #). 0 => highest unique.
  OUTDISK.....Disk drive # of Output image. 0 => highest w room.
  PBPARM......Primary beam parameters:
              (1) Lowest beam value to believe: <= 0.001 -> 0.023
                  Image pixels outside this range are blanked.
                  >= 1 => do not do a primary beam correction.
              (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
  COORDINA....The RA and DEC coordinates for the pointing position are
              found as:
              RA = abs(CO(1)) + abs(CO(2))/60 + abs(CO(3))/3600
                 if any of CO(1), CO(2), CO(3) < 0: Xpos = -Xpos
              DEC = abs(CO(4)) + abs(CO(5))/60 + abs(CO(6))/3600
                 if any of CO(4), CO(5), CO(6) < 0: Ypos = -Ypos
              RA is in hours, minutes, seconds of time
              DEC is in degrees, minutes, seconds of arc
              WARNING: if any of these 6 are not zero, then COORDINA
              defines the pointing direction for the image.  Otherwise
              the pointing position in the header (or if that is 0,
              the source coordinate in the header) is used.
  RADIUS .....Radius in pixels over which the spectral index images
              are smoothed.
  DOBLANK.....> 0 => blank any pixels for which the spectral index
                     image or the curvature (if any) are blanked
              <=0 => do no spectral index correction at pixels for
                     which there is a blanked pixewl in the spectral
                     index images.
  DOINVERS....> 0 Multiply the image by the beam pattern.
              <= 0 Divide the image by the beam pattern.
----------------------------------------------------------------
SPCOR: Task to apply the primary beam correction
RELATED PROGRAMS:  PRTIM

                          PURPOSE

     SPCOR 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


                 RICK PERLEY'S (OLD) REPORT

	Polynomial Coefficients from LSq Fit to VLA Primary
	Beam raster scans.

	Functional form fitted:

		1 + G1.X^2 + G2.X^4 + G3.X^6

	where X = r.F,

	and 	r = radius in arcminutes
		F = frequency in GHz.

	Fits were made to 3% cutoff in power for 24 antennas.
Poor fits, and discrepant fits were discarded, and the most
consistent subset of antennas had their fitted coefficients
averaged to produce the following 'best' coefficients.


Freq.		G1		G2		G3
----------------------------------------------------------
1.285           -1.329E-3       6.445E-7        -1.146E-10  *
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 (old)    -1.350          6.526           -1.090      *
22.485 (new)    -1.417          7.332           -1.352
43.315          -1.321          6.185           -0.983
-----------------------------------------------------------

	The estimated errors (from the scatter in the fitted
coefficients) are generally very small:

	G1: .003 at all bands except Q (.014)
	G2: .03 to .07 at all bands except Q (.15)
	G3: .01 to .02 at all bands except Q (.04)

	R. Perley  21/Nov/00

* The 1.285 and 22.485 old feed values are not used.

                        COMMENTS

COORDINA:
     If COORDINA is all zero, the pointing position is taken to be
the pointing position given in the header or, if that is 0,
the position of the reference pixel (equal to the map making
position which is equal to the phase tracking position if the
map was made with no shift).  Any other pointing position can be
specified in COORDINA, an array of dimension 6 in the format

            COORDINA= HH MM SS.SS DD MM SS.SS

where the RA is HHh MMm SS.SSs and the DEC is DDd MM' SS.SS".
Set some or all of COORDINA to non-zero to have it used by the
program.

ACCURACY:
     The accuracy of the correction is about 1% within the 50%
primary beam sensitivity and about 2% beyond this region.  Thus
the errors of the correction F(x) become quite large near the
edges of the primary beam response.

WHEN TO USE SPCOR:
     Do NOT use SPCOR on the dirty map before running APCLN or
any of the other deconvolution processes.  SPCOR is usually
run near the end of data processing.  It may now be run on
data "cubes" but it does require that frequency, right
ascension, and declination axes appear in the header.

     The speed of the task is dependent on the axis order and
the accuracy of the position computations.  In general, with
linear computations, the RA-DEC-FREQ order of axes is faster.
With non-linear position computations (needed for 1-degree
fields and larger), the FREQ-RA-DEC order appears to be
significantly faster than RA-DEC-FREQ.  Any order will give
correct results, however, and non-linear is always rather
slower than linear.
