; PBEAM
;---------------------------------------------------------------
;! Fits the analytic function to the measured values of the beam
;# Task Calibration
;-----------------------------------------------------------------------
;;  Copyright (C) 1995-1997, 2000, 2005, 2008, 2012, 2015-2016, 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
;-----------------------------------------------------------------------
PBEAM     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
PBEAM     Fits an analytic function to raster-scanned beam data.
OUTNAME                            Save beam images if not blank
OUTSEQ             0.0     9999.0  Output files name (seq. #)
OUTDISK                            Output files disk unit #
INFILE                             First input data file
IN2FILE                            Second input data file
OUTTEXT                            Output file to record the
                                   fitted parameters
VPARM                              1. Fractional power cutoff
                                      for fit
                                      0 => 0.05 (5% of peak)
                                   2. Antenna diameter, meters
                                      0 => 25
                                   3. number of iterations for
                                      nonlinear least square fit
                                      0 => 4
                                   4. Degree of the fit polynom.
                                      0 => 3
                                   5. exclude data with
                                      |phase| > APARM(5) degrees
                                      0 => The cutoff = 60 deg.
                                      -1 =>do not exclude points
                                           on the basis of phase
                                   6. 0 => combine the two input
                                           files (Stokes 'I')
                                      1 => Fit model to each
                                           input file separately
                                           and write out
                                           differences (squint)
                                   7. 0 => do not change the
                                           input files
                                      1 => switch the first two
                                           columns of the input
                                           files
                                   8. > 0 -> Complex normalize
                                      else normalize by peak amp
                                   9. width factor to convolve
                                      to grid 0 -> 1.4,
                                      smaller is smoother, wider
                                   10. > 0 Do not normalize the
                                      fit parameters
                                   11. > 0 => data are power not
                                      voltage
BPARM                              These are used to simulate
                                   data for testing the program.
                                   1. Coefficient of the first
                                      degree term (G1)
                                      0 => -0.1
                                   2. Coefficient of the second
                                      degree term (G2)
                                      0 => 0.0034
                                   3. Coefficient of the third
                                      degree term (G3)
                                      0 => -4.60E-05
                                   4. XSHIFT - horizontal beam
                                      offset
                                      0 => 0.1 arc minutes
                                   5. YSHIFT - vertical beam
                                      offset
                                      0 => 0.05 arc minutes
                                   6. eccentricity
                                      0 => 0.1
                                   7. Error of eccent. guess
                                      0 => 0
                                   8. Major axis position angle,
                                      deg.
                                      0 => 60
                                   9. Error of position
                                      angle guess, degrees
                                      0 => 5
DOPLOT                             1,3,5,7 -> make contour plots
                                   2,3,6,7 -> make IRING plots
                                   4,5,6,7 -> make IRING
                                              difference plots
OPTYPE                             '    ' => observed data used,
                                   eccentricity and Pos. angle
                                   are not fitted
                                   'ELLI' => observed data used,
                                   eccentricity and Pos. angle
                                   are  fitted.
                                   'SIMU' => simulated data are
                                   used; eccentricity and Pos.
                                   angle are  fitted.
PRTLEV                             0 => print the polynomial
                                        coefficients and errors.
                                   1 => print observed data and
                                        fitted values as well.
                                        X,Y coordinates use
                                   2 => print observed data and
                                        fitted values as well.
                                        RO, PHI coordinates use
                                        Header is not printed
DOTV             -1.0       1.0    > 0 Do plot on the TV, else
                                   make a plot file
GRCHAN            0.0       8.0    Graphics channel 0 => 1.
BADDISK           0.0     9999.0   Disks to avoid for scratch
----------------------------------------------------------------
PBEAM
Type: Task
Use:  Fits an analytic function to raster scan beam data.  It can also
      make contour plots of the data, the fit model, and the residual
      (data-model) and can output interpolated images of them.  In
      addition, it can make plots of the data and model in rings about
      the fit center (ring-by-ring and sum interior to ring) and plots
      of the same rings in the residual image.
Adverbs:
  OUTNAME....Output image file name (name).  If not blank, output
             images will be made of the data, model, and residual.
             If DOTV <= 0, OUTNAME must be set and the initial plots
             of the image will be attached to the image files.
             If OUTNAME is blank, DOTV will be forced to true.
  OUTSEQ.....Output UV file name (seq. #).    0 => highest.
  OUTDISK....Disk drive # of output UV file.  0 => any.
  INFILE.....Input file with the raster scanned beam data for one
             polarization.  The accepted format is that used by UVHOL.
             Lines starting with '#!' are headers, and contain a
             variety of information.  The lines beginning '#! RefAnt'
             are particularly inportant and should not be altered from
             those provided by UVHOL.  The give reference antenna,
             antenna number, Stokes and frequency for the following
             beam data. The data lines have a free format giving X, Y,
             Amplitude, and Phase.
  IN2FILE....Input file with beam data for the 2nd stokes (optional).
  OUTTEXT....Output file to record the fitted parameters and data.
  VPARM..... 1. Minimum normalized power level used.  0 => 0.05 (5%)
             2. Diameter of antenna, meters. 0 => 25
             3. number of iterations for nonlinear least square.
                0 => 4
             4. Degree of the fitted polynomial. 0 => 3
             5. Exclude points with abs(phase) > VPARM(5) degrees
                0 => The cutoff is equal 60 degrees
                .LE. 0 => do not exclude points on the basis of phase
             6. 0 => combine the two input files and fit to Stokes 'I'
                1 => Fit the model to both input files separately,
                     and compute the beam squint parameters. The INFILE
                     data define the reference position.
             7. PBEAM expects the first column of the input files to
                be the horizontal offset and the second to be the
                vertical offset.  If your data file has the opposite
                convention, set VPARM(7)=1 to switch the columns
                before fitting.
             8. The incoming data are normalized by the value of the
                peak amplitude.  This can be done in a complex fashion
                using the amp and phase at the peak iff VPARM(8) > 0.
             9. The images are constructed by smoothing to a grid with
                a circular Gaussian function.  The width of that
                function is controlled by VPARM(9) with default 1.4.
                Larger values mean that the data samples do not smooth
                together and so by 2.0 the image is rather "boxy"
                around each sample.  Smaller values smooth the image
                further but also fatten the beam.  This does not
                affect the fit results, only the displays.
             10. The beam fit parameters are usually normalized so the
                G0 parameter has value 1.0.  To omit this
                normalization, set VPARM(10) > 0.
             11. When the UVHOL output is for one moving and one fixed
                antenna, then the values are effectively voltages and
                must be squared to fit the beam.  However, when both
                antennas are moving, then the values are visibilities
                or power and the data should not be squared.  Set
                VPARM(11) > 0 to indicate this case.
  BPARM......Beam simulation parameters.
             1. Coefficient of G1, the first degree term; 0 => -0.1
             2. Coefficient of G2, the second degree term;  0 => 0.0034
             3. Coefficient of G3, the third degree term; 0 => -4.60E-05
             4. XSHIFT; 0 => 0.1 (Beam offset in azimuth)
             5. YSHIFT; 0 => 0.05 (Beam offset in elevation)
             6. Eccentricity; 0 => 0.1
             7. Error of the eccentricity's guess
             8. Position angle, degrees; 0 => 60
             9. Error of the position angle guess, degrees; 0 => 5
             	G0 is fixed at 1.0.
  DOPLOT.....Controls whether plots are done:
             DOPLOT = 1,3,5,7 -> make contour plots of each image
                                 (data, model, residual)
             DOPLOT = 2,3,6,7 -> make IRING-like plots comparing data
                                 and model with average beam in ring
                                 and sum of beam values inside ring
             DOPLOT = 4,5,6,7 -> make plots of the IRING-like average
                                 in each ring and sum interior from
                                 the difference image and of the
                                 difference IRING divided by the
                                 model IRING.
  OPTYPE.....'    ' => the given number of the polynom's coefficients,
                       and the X-, Y-shifts are fitted to the data
             'ELLI' => the given number of the polynom's coefficients,
                       the X-, Y-shifts, eccentricity and position
                       angle  are fitted to the observed data
             'SIMU' => the given number of the polynom's coefficients,
                       the X-, Y-shifts, eccentricity and position
                       angle  are fitted to the simulated data
             The position angle is defined as the angle between the
             the major axis and the Y (elevation) axis.  A positive
             angle is defined as the rotation which turns the positive
             X axis into the positive Y axis.
  PRTLEV.....0 => print only the coefficients and errors of the fit.
             1 => print observed and fitted data
                  X,Y coordinates use
             2 => print observed and fitted data
                  RO, PHI coordinates use
                  Header is not printed
  DOTV.......> 0 Do plot on the TV.
  GRCHAN.....Graphics channel 0 => use channels 1,2,3,4 which allows
             color to separate things being plotted.  Labels go in 1
             (usually yellow), contours and the IRING connecting line
             go in 2 (usually green), data IRING X's go in 4 (usually
             cyan), and model IRING plus signs go in 3 (usually pink).
  BADDISK....A list of disks on which scratch files are not to
             be placed.  This will not affect the output file.
------------------------------------------------------------------------
PBEAM:  Fits for the coefficients of an analytical model to
        raster-scanned data produced by UVHOL.
Documentors: L.R. Kogan, R.A. Perley

The program fits a polynomial of the form:

		GO + G1*ARG + G2*(ARG^2) + G3*(ARG^3) + ...

	where ARG = (X1*EC)^2 + Y1^2
              X1 =  (X-X0)*COS(PHI) + (Y-Y0)*SIN(PHI)
              Y1 = -(X-X0)*SIN(PHI) + (Y-Y0)*COS(PHI)
              X0 and Y0 are offsets from the fiducial center.
              EC = b/a is eccentricity of the ellipse
              PHI is position angle of the major axis

In the case of circular beam (EC and PHI are not fitted) the fit
polynomial looks simpler:


		GO + G1*ARG + G2*(ARG^2) + G3*(ARG^3) + ...

	where ARG = (X-X0)^2 + (Y-Y0)^2
              X0 and Y0 are offsets from the fiducial center.

   The order of the polynomial is specified by the user via VPARM(4).
Non linear least square method is used for the fitting.
The number of iterations in the non linear least square fitting
is specified by the user via VPARM(3).
The default for VPARM(3)=4 and VPARM(4)=3 works fine in most cases.
	
   The task reads from external files produced, for example, by
the AIPS task UVHOL.  UVHOL processes VLA holography data produced
by the VLA's online holography programs.
   PBEAM assumes, but does not require, calibrated data.
Normally, this means the maximum amplitude is 1.0, and the phase is
zero.  PBEAM searches for the maximum amplitude in the raster data,
and normalizes all measurements by this value.  (This is necessary
because large pointing errors will throw off the amplitude
calibration).
The phase is used to discriminate against data taken in the first sidelobe
This test can be turned off by setting VPARM(5)=1.
   After normalizing the amplitude, PBEAM squares the amplitudes
and applies a cutoff given by VPARM(1).  High sidelobe data are rejected
via the phase test, if desired.   The requested model is then applied to
the surviving data.

   If the model is to be applied to only a single polarization,
IN2FILE must be left blank.
   Both circular polarizations may be processed simultaneously
if they are specified via INFILE and IN2FILE.  Although it is not
necessary for these data to be on the same grid, the only points which
will be utilized in the analysis will be those which are taken at
very similar grid points.

There are two distinct modes available when both polarizations are
provided:
	If VPARM(6) = 0, the input files are vector averaged (normally
producing a Stokes' I value) and the model applied.  The X and
Y offsets found then represent the pointing error of the antenna.
	If VPARM(6) = 1, the two input files are separately analyzed
(producing offsets for the two polarizations separately).  The program
then takes the difference in each axis, and writes out the beam squint
amplitude and position angle.
	The sign of the offsets is such that a positive XSHIFT means
the antenna beam is located at a greater azimuth than the reference
position.  A positive YSHIFT means the antenna beam is at a higher
elevation than the reference position.  The squint angle (PHISHIFT)
is the p.a. of the IN2FILE beam's maximum w.r.t. the INFILE beam
maximum, with the horizontal axis defining the zero, and a CCW
direction defined as positive.  The units of XSHIFT, YSHIFT, and
ROSHIFT (squint separation) are in arc minutes.  Error estimates,
derived from the errors in the R and L beam centers, of the squint
parameters are produced and printed.

	The program first divides the data by a rough estimate of the
beam maximum.  The analytic model described above is then fitted.
The coefficients are then divided by G0 (the fitted beam maximum) and
printed. Error estimates are also given, after normalization by the
value of G0.
	The post-fit RMS is printed.
	VPARM(5) is a phase cutoff, intended to prohibit high sidelobe
values from entering the analysis.  Only use this if the data are
calibrated (so the main-beam values are near 0 degrees).
	The BPARMs are used for generating fake beam data, for testing
the program.  This option is enabled when OPTYPE = 'SIMU'.
	The program will produce PL files, showing the input data and
the fitted curves.
	The program will fit for beam ellipticity by setting
OPTYPE = 'ELLI'.  It is advised to use a high power cutoff (say, 40%)
to avoid skewing the results by the quadrupod structures.

PBEAM can work with several reference antennas but only with one
investigared antenna. So the input file should include data only
for one antenna.

Three images are written when OUTNAME is not blank - the data, the
model, and the residual (data-model).  They are accompanied by plot
files if DOTV <= 0.  If OUTNAME is blank or DOTV > 0, the plots are
displayed on the TV.  The plots are contour images with LEVS from
-0.95 to +0.95 in steps of 0.1 for the data and model images and from
-0.145 to +0.145 in steps of 0.01 for the residual image.  The plots
also put little plus signs where the data samples actually occurred.
