; GAL
;---------------------------------------------------------------
;! Determine parameters from a velocity field
;# TASK ANALYSIS
;-----------------------------------------------------------------------
;;  Copyright (C) 1995, 1998
;;  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
;-----------------------------------------------------------------------
GAL       LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
GAL Task to determine parameters from a velocity field.
USERID      -32000.0      32000.0  User ID.  0=>current user
                                     32000=>all users
INNAME                             Image name
INCLASS                            Image class
INSEQ            0.0       9999.0  Image seq. #
INDISK           0.0          9.0  Disk drive #
IN2NAME                            Weight image name
IN2CLASS                           Weight image class
IN2SEQ           0.0       9999.0  Weight image seq. #
IN2DISK          0.0          9.0  Weight image disk drive #
OUTNAME                            Output image name (name)
OUTCLASS                           Output image name (class)
OUTSEQ           0.0       9999.0  Output image name (seq. #)
OUTDISK          0.0          9.0  Output disk drive #
BLC              0.0       4096.0  Bottom left corner of image
                                     0=>entire image
TRC              0.0       4096.0  Top right corner of image
                                     0=>entire image
FACTOR           0.0          1.0  Convergence criterion.
                                     0=> 0.001	
FUNCTYPE                           Type of rotation curve: 'BR':
                                   Brandt, 'EX' : exponential,
                                   'CC' : constant, 'SB' : solid
                                   body.
APARM                              First guesses to parameters:
                                   1: x - coordinate centre, 2:
                                   y - coordinate centre, 3: po-
                                   sition angle major axis, 4:
                                   inclination, 5: systemic ve-
                                   locity, 6: maximum velocity,
                                   7: radius maximum velocity
                                   ("), 8: exponent (0=>1).
CPARM                              (1) : residual map; 0 : none,
                                   1 : uses whole field, 2 :
                                   uses selected field, 3: uses
                                   values from INFILE. (2) :
                                   0/1 : plot (no/yes). (3) :
                                   fixed parameters. (4) : 0/1:
                                   use IN2C map for weighting
                                   (no/yes) 5: 0/1: copy old
                                   history (no/yes) 6: 0: write
                                   one line to output;-1/1 write
                                   position to output file.
DPARM                              1 : lower limit radius; 2 :
                                   upper limit radius; 3 : lower
                                   limit cos (pos. angle in ga-
                                   laxy); 4 : upper limit; 5 and
                                   6 : idem for sine; 7 : max.
                                   plotted radius; 8: max. plot-
                                   ted rot. velocity; 9 :0 (def)
                                   => fit parameters, 1=> just
                                   plot using inputs. 10: width
                                   annuli (in pixel separation).
INFILE
                                   Text file containing rotation
                                   curve information.
OUTFILE
                                   Text file to which one line
                                   of output results is written.
PIXSTD                             Estimated rms uncertainty in
                                   the observed radial velocity
                                   (km/s) at one pixel.  0=>10
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.
----------------------------------------------------------------
GAL
Type: Task
Use:  GAL calculates eight parameters specifying the orientation (# 1 to
      5) and the rotation curve (# 6 to 8) of a galaxy, using its
      velocity field as input. If this velocity field consists of M data
      points, the parameters are determined by a nonlinear least squares
      fitting of eight parameters to M data points. If CPARM(1)> 0, a
      residual velocity field (observed minus model) is produced.  If
      CPARM(2)> 0, the observed rotation curve is plotted together with
      the fitted model curve.
Adverbs:
  USERID......User ID of owner of image.  0=>current user
                32000=>all users
  INNAME......Image name(name).  blank=>any
  INCLASS.....Image name(class).  blank=>any
  INSEQ.......Image name(seq. #).  0=>any
  INDISK......Disk drive # of image.  0=>any
  IN2NAME.....Weight image name(name).  blank=>any
  IN2CLASS....Weight image name(class).  blank=>any
  IN2SEQ......Weight image name(seq. #).  0=>any
  IN2DISK.....Disk drive # of weight image.  0=>any
  OUTNAME.....Output image name(name).   blank=>INNAME
  OUTCLASS....Output image name(class).  blank=>'GAL'
  OUTSEQ......Output image name(seq).    0=>next unique
  OUTDISK.....Output image disk drive #. 0=>any
  BLC.........The Bottom Left-hand pixel of the subarray of
              the image to be analysed. The value (0,0)
              means (1,1).
  TRC.........The Top Right-hand pixel of the subarray of
              the image to be analysed. (0,0) => (1,1).
              the top right hand corner of the entire image.
  FACTOR......Criterion to stop least squares fitting.
  FUNCTYPE....Type of rotation curve to be fitted.
                'BR' : Brandt curve, given by :
                                          R / Rmax
                   V / Vmax   =   ---------------------- n 3/2n
                                  (1/3 + 2/3 * (R / Rmax) )
                'EX' : Exponential flat curve :
                                        - ln(100.0) * (R / Rmax)
                   V / Vmax   =   1 - e
                'CC' : Constant curve :    V / Vmax   =   1
                'SB' : Solid Body curve :  V / Vmax   =   R / 60
  APARM.......Initial guesses to the eight parameters to be
	      fitted : (1) : x pixel central position; (2) : y
              pixel central position; (3) : position angle re-
              ceding major axis; (4) : inclination orbital
              plane; (5) : systemic velocity; (6) : Vmax, max.
              rot. velocity; (7) : Rmax, radius(") of maximum
              rot. velocity; (8) : n, measure of extent of
              'flat' part of curve, 0=>1.  (8) is only used if
              FUNCTYPE='BR', and (7) only if FUNCTYPE='BR' or
              'EX'.
  CPARM       (1) : if 1, a residual output map is made using
              the whole field, if 2, this is done for the speci-
              fied part of the field only; (2) : if 1, the ob-
              served and the model rotation curves are plotted.
              (3) : parameters to be held fixed. e.g. if para-
              meters i, j, and k are to be held fixed,  CPARM(3)
              = 2**(i-1) + 2**(j-1) + 2**(k-1). 0 defaults to NO
              fixed parameters. (4) 1: use IN2C map for weights,
              0: don't use weighting image. (5) 1: copy old his-
              tory file; 0: don't copy it. (6) 0: write one line
              with fitted parameters to OUTFILE; 1/-1 write pos-
              ition of major axis to OUTFILE (1 for receding end
              -1 for approaching end).
  DPARM       (1) and (2): lower and upper limit to radius (");
              (3) and (4): limits to cos(azimuth angle); (5)
              and (6): limits to sin(azimuth angle); (7) and
              (8) : max. radius and max velocity to be plotted;
              (9) : no fit, residual field only; (10) : width
              annuli in plotting observed rotation curve.
  INFILE......Text file containing rotation curve information.
              This is used to determine the residual field when
              CPARM(1)=3.
  OUTFILE.....Text file to which one line of output results is
              written. When OUTFILE exists, GAL will append.
              Meant to be used in loop when using concentric
              rings; OUTFILE can then become INFILE of next run.
  PIXSTD......Estimated rms uncertainty in observed radial
              velocity at one pixel (km/s)
  DOTV........> 0 => plot directly on the TV device, otherwise
              make a plot file for later display on one or
              more devices (including the TV if desired).
  GRCHAN......Graphics channel (1 - 7) to use for line drawing.
              0 => 1.
----------------------------------------------------------------
GAL : Task which analyses a velocity field.  In the present im-
plementation it is assumed that the object is a rotating disk.
The user has to specify initial guesses to the parameters to be
fitted to the velocity field, using the adverb APARM.  GAL will
then determine the best least squares' fit to the input velocity
field, and give the final best fit values of the parameters.  Of
course a good initial guess will speed up the procedure, whereas
a very bad one will cause divergence and failure.

Some remarks on specific adverbs :

      IN2NAME, etc.
Specifies a map of which the pixel values are weights to be
applied to the velocities in the velocity field map INNAME. It
is up to the user to supply such a map. Several possibilities
are open: one can use the total HI map, or the square of the
total HI map (using COMB). The latter choice is the correct one
if the noise in the velocity maps can be assumed Gaussian with
respect to the intensity map which normally is a reasonable
assumption. Other possibilities are the width of the profile,
the skewness of the profile, or - if the velocity field was
determined using XGAUS - the goodness of fit.
      BLC , TRC The velocity field is read in as a whole, which
imposes an upper limit to the area that may be used. A 512*512
velocity should be handled easily. When your image is too big, make
sure you are not oversampling too much; more than two points per beam
is unnecessary.
      FACTOR The value of FACTOR determines when to stop the
iteration, the default 0.001 should be adequate for most purposes.
      FUNCTYPE Four schematic rotation curves are available now.  The
Brandt curve (FUNCTYPE = 'BR') attains a maximum value Vmax at R =
Rmax, and declines Keplerian thereafter.  The exponential curve
(FUNCTYPE = 'EX') is flat throughout after an initially linear
increasing part.  The factor ln(100) in the exponent was chosen to
give analogous meanings to Rmax for the two representations, in fact
Rmax is the radius where Vmax attains 99% of its maximum value.
Still, the fitted value of Rmax may well differ depending on the value
of FUNCTYPE.  FUNCTYPE='CC' fits a one parameter rotation curve, V =
constant. This is a good approximation in many cases, but a minimum
radius DPARM(1) should be specified since the approximation breaks
down at low radii. It is particularly useful in applying GAL in
concentric annuli. GAL then should be used in a loop with varying
DPARM(1) and DPARM(2) fitting a constant velocity in each ring.
FUNCTYPE='SB' fits a solid body rotation curve.  In this case the
theoretical velocity field has an infinite number of axes of symmetry,
and the central position becomes undefined.  Furthermore, the
steepness and the inclination are directly related, so one of both
(the inclination) is held fixed at its initial value.  Vmax is defined
as the circular velocity at a 60" radius.  GAL keeps x,y, and i
fixed, although of course you may add other fixed parameters by
means of CPARM(3).
      APARM(1),APARM(2)
The guess to the central position has to be given in pixels; the
resulting fitted value is also displayed in RA and Dec.
      APARM(3)
The position angle is defined measured form north eastward, to
the receding part of the major axis (having the largest veloci-
ties).  An initial guess which is 180 degrees off is bound to
cause divergence.
      APARM(4)
The inclination angle can have values from 0 (face on) to 90
(edge on). Don't use values close to 0 or 90 as initial guesses.
      APARM(5)
A fair initial guess is the mean of the velocities on opposite
ends of the major axis, or values around the minor axis. Specify
velocities in km/s. Note that the units in the map should be in m/s.
      APARM(6)
The height of the rotation curve can be guessed by halving the
difference of two points on opposite ends of the major axis,
and applying an cosec(i) correction.
      APARM(7)
The extent of the rising part may be hard to guess, but usually
the program converges if the initial guess is less off than a
factor 2 or so.
      APARM(8) This is a measure of the flatness of the curve, if
FUNCTYPE = 'BR'.  Very flat curves have APARM(8) well below unity,
while curves with a clear decline in the outer parts have APARM(8)
values like 1.5, 2.0, or even 3.0. If uncertain, start with APARM(8)
= 1. Avoid APARM(8) < 0.1 or > 4.

A successful convergence in the 'BR' case depends critically on
reasonable choices for both APARM(7) and APARM(8).

      CPARM(1)
controls the residual field. If CPARM(1) = 1, 2, or 3, a residual
output map is made. Such a map can be extremely useful to assess
the quality of the fit, and to detect the presence of asymmetries.
In the case CPARM(1)=1, the whole velocity field is used to
generate the residual map, irrespective of the values in DPARM(1)
through DPARM(6). When CPARM(1)=2, only the part of the velocity
field specified in DPARM(1) => DPARM(6) is used. When CPARM(1)=3,
no fitting is done, and the residual field is determined using a
user supplied rotation curve in the textfile INFILE.

      CPARM(2)
If 1, the observed rotation curve is plotted (by integrating the
circular velocities in annuli), as well as the fitted model curve.

      CPARM(3)
is a bit map of the parameters to be held fixed.  If one wants
to fix the position ( both x and y), and Rmax, specify CPARM(3) =
2**0 + 2**1 + 2**6 = 67. CPARM(3)=0 means NO fixed parameters.
Note that some parameters may be fixed by the program anyway
(x,y, and i if FUNCTYPE='SB').

      CPARM(4)
If 0, don't use the image in IN2C as weighting image; if 1, use the
IN2C image (if specified).

      CPARM(5)
If 0, don't copy the old history file; if 1: copy it. The reason
for this is that in the case of very large history files a substan-
tial fraction of the execution time of GAL is spent in copying the
file.
      CPARM(6)
Only important when OUTFILE is specified. A one liner of information
about the fit obtained is appended to OUTFILE. If CPARM(6)=0, this
one line contains the values of all fitted parameters. If CPARM(6)
= 1/-1, the one line of output contains the position (RA, Dec) of
the major axis of the ring under consideration. When +1, GAL will
write the position of the receding major axis, and when -1, of the
approaching major axis. When using GAL many times in rings, the
output file contains a range of positions which trace the major axis
of the galaxy. This file can than be turned in a ST extension file,
and the major axis can be plotted on top of other maps.
      DPARM
The parameters DPARM(1) through DPARM(6) allow the user to spe-
cify a wide variety of subsections of the plane of the galaxy.
DPARM(1) and DPARM(2) specify the minimum and maximum radius to
be used, e.g. to select an annulus, in which case FUNCTYPE='CC'
or 'SB' would be advisable.  DPARM(3) and DPARM(4) specify li-
mits on the cosine of the azimuth angle in the galaxy (as pro-
jected on the sky), e.g. 0,1 selects the receding half of the
galaxy.  DPARM(5) and DPARM(6) specify the sine of this angle,
e.g. 0,1 selects the half plane traversed when the receding
major axis is rotated in an anticlockwise direction to the ap-
proaching major axis.  DPARM(7) specifies the maximum radius in
the plot (0=> auto scaling).  DPARM(8) specifies the maximum
rotation velocity to be plotted (0=> auto scaling).  If DPARM(9)
= 1, no fitting is done, only a plot is generated, based on the
input values in APARM.  DPARM(10) specifies the width of the
annuli used in plotting the actual rotation curve, approxima-
tely in units of a pixel width (if DX=DY : exactly one pixel
width, DX>>DY : DY * SQRT(2), DX<<DY : DX * SQRT(2)).  A value
of 0 defaults to 3.

      INFILE
This is a text file containing rotation curve information
from concentric rings. This file is used when CPARM(1) = 3. No
fitting is done, but a residual map is produced based on the
contents of INFILE. A sample INFILE:

  30 224.5 256.5 120.8 41.5 832 120
  80 224.5 256.5 124.3 43.5 832 210

This INFILE defines two rings, one from 0" to 30", where the
parameter values are given on the first line. The second ring
goes from 30" to 80"; its parameter values are given on the
second line. Such a file typically is the result of several GAL
runs on annuli, with FUNCTYPE = 'CC'.

      OUTFILE
Output results are written to this file, if specified. If
CPARM(6)=0, output parameters are the radii of the ring used,
plus the fitted parameters. To this are added the average dif-
ference between the fitted model and the velocity field (should
be close to zero), and the standard deviation (should be close
to the rms in the ring). If CPARM(6)=+/-1, the location of the
major axis is written (see CPARM(6) explanation).

      PIXSTD
PIXSTD is used to calculate the errors associated with the
final results.  The default is the velocity dispersion usually
observed in unperturbed spirals, but the value may well be
higher as a result of larger random motions and filter width.


Hints
=====

The curve types 'BR', and 'EX' should only be used on the whole
field, and in the initial stage of determining the variables.
Real-world rotation curves deviate from mathematical models. In
the more advanced stages of velocity field analysis, it is
suggested to use 'CC' only, on annuli determined by DPARM(1)
and DPARM(2). The width of these annuli should be a trade-off
between a minimum number of points in the ring to make a good
fit, and a maximum width over which the rotation velocity may
be assumed constant. The whole rotation curve can then be
fitted by a command like:

for i=1:10;dparm(1)=i*30;dparm(2)=(i+1)*30;go gal;wait gal;end

A residual of this fit can be obtained by putting the fitted
values in the textfile INFILE, and running GAL once more with
CPARM(1)=3. Note that in the above example the very central
portion of the velocity field is not used.



