%-----------------------------------------------------------------------
%;  Copyright (C) 2014-2015, 2017, 2020, 2025
%;  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
%-----------------------------------------------------------------------
\documentclass[twoside]{article}
\usepackage{palatino}
\renewcommand{\ttdefault}{cmtt}
% Highlight new text.
\usepackage{color}
\usepackage{alltt}
\usepackage{graphicx,xspace,wrapfig}
\usepackage{pstricks}  % added by Greisen
\definecolor{hicol}{rgb}{0.7,0.1,0.1}
\definecolor{mecol}{rgb}{0.2,0.2,0.8}
\definecolor{excol}{rgb}{0.1,0.4,0.1}
\newcommand{\Hi}[1]{\textcolor{hicol}{#1}}
%\newcommand{\Hi}[1]{\textcolor{black}{#1}}
\newcommand{\Me}[1]{\textcolor{mecol}{#1}}
%\newcommand{\Me}[1]{\textcolor{black}{#1}}
\newcommand{\Ex}[1]{\textcolor{excol}{#1}}
%\newcommand{\Ex}[1]{\textcolor{black}{#1}}
\newcommand{\No}[1]{\textcolor{black}{#1}}
\newcommand{\hicol}{\color{hicol}}
%\newcommand{\hicol}{\color{black}}
\newcommand{\mecol}{\color{mecol}}
%\newcommand{\mecol}{\color{black}}
\newcommand{\excol}{\color{excol}}
%\newcommand{\excol}{\color{black}}
\newcommand{\hblack}{\color{black}}
%
\newcommand{\AIPS}{{$\cal AIPS\/$}}
\newcommand{\eg}{{\it e.g.},}
\newcommand{\ie}{{\it i.e.},}
\newcommand{\etal}{{\it et al.}}
\newcommand{\tablerowgapbefore}{-1ex}
\newcommand{\tablerowgapafter}{1ex}
\newcommand{\keyw}[1]{\hbox{{\tt #1}}}
\newcommand{\sub}[1]{_\mathrm{#1}}
\newcommand{\degr}{^{\circ}}
\newcommand{\vv}{v}
%\newcommand{\vv}{\varv}
\newcommand{\eq}{\hbox{\hspace{0.6em}=\hspace{0.6em}}}
\newcommand{\newfig}[2]{\includegraphics[width=#1]{data.fig#2}}
%\newcommand{\putfig}[1]{\includegraphics{data.fig#1.eps}}
\newcommand{\putfig}[1]{\includegraphics{#1.eps}}
\newcommand{\whatmem}{\AIPS\ Memo \memnum}
\newcommand{\boxit}[3]{\vbox{\hrule height#1\hbox{\vrule width#1\kern#2%
\vbox{\kern#2{#3}\kern#2}\kern#2\vrule width#1}\hrule height#1}}
%
\newcommand{\memnum}{118 \Hi{Revised}}
\newcommand{\memtit}{Modeling Spectral Cubes in \AIPS}
\title{
   \vskip -35pt
   \fbox{{\large\whatmem}} \\
   \vskip 28pt
%   \vskip 10pt
%   \fbox{{\Huge
%              \Ex{D R A F T}}}
%   \vskip 10pt
   \memtit \\}
\author{Eric W. Greisen}
%
\parskip 4mm
\linewidth 6.5in                     % was 6.5
\textwidth 6.5in                     % text width excluding margin 6.5
\textheight 9.0 in                  % was 8.81
\marginparsep 0in
\oddsidemargin .25in                 % EWG from -.25
\evensidemargin -.25in
\topmargin -0.4in
%\topmargin 0.2in
\headsep 0.25in
\headheight 0.25in
\parindent 0in
\newcommand{\normalstyle}{\baselineskip 4mm \parskip 2mm \normalsize}
\newcommand{\tablestyle}{\baselineskip 2mm \parskip 1mm \small }
%
%
\begin{document}

\pagestyle{myheadings}
\thispagestyle{empty}

\newcommand{\Rheading}{\whatmem \hfill \memtit \hfill Page~~}
\newcommand{\Lheading}{~~Page \hfill \memtit \hfill \whatmem}
\markboth{\Lheading}{\Rheading}
%
\vskip -.5cm
\pretolerance 10000
\listparindent 0cm
\labelsep 0cm
%
%

\vskip -30pt
\maketitle

\normalstyle
\begin{abstract}
  \AIPS\ has done Gaussian fitting along the $x$-axis of image cubes
  with task \keyw{XGAUS} since the 1980s.  That task \Hi{was overhauled
    in 2013} to be much easier to use and much more capable.  In like
  fashion, tasks \keyw{ZEMAN} and \keyw{RMFIT} were developed.  The
  former fits the standard leakage and scaling terms for Stokes V cubes,
  including a new option to do this for each of the Gaussians found by
  \keyw{XGAUS}\@.  The latter fits polarization models to Stokes Q and U
  cubes, using the output of Faraday Rotation Measure Synthesis
  (\AIPS\ task \keyw{FARS}) to assist with initial guesses.  The models
  can contain multiple components each with a polarization flux, angle,
  rotation measure, and rotation measure ``thickness.''  The present
  memo will describe the functions of these tasks in some detail with
  numerous graphical examples.  A revision of this memo discusses
  changes made for the {\tt 31DEC15} release and two new tasks to plot
  spectra with model fits and a number of tasks which make visibility
  and image model files.  The November 2017 revision includes new
  display options in {\tt RMFIT}\@.  \Hi{The April and July 2025
    revisions add display and edit options and correct display and
    pixrange setting.}
  \footnote{Color is used to indicate changes made to this document:
    the all revisions prior to 2025 \Hi{and the April and July 2025
      revisions}.}
\end{abstract}

\renewcommand{\floatpagefraction}{0.75}
\typeout{bottomnumber = \arabic{bottomnumber} \bottomfraction}
\typeout{topnumber = \arabic{topnumber} \topfraction}
\typeout{totalnumber = \arabic{totalnumber} \textfraction\ \floatpagefraction}

\section{Introduction}

{\tt XGAUS} was written, and re-written, a long time ago in the hopes
that converting a spectral cube, with many spectral channels at each
spatial location, into one or more components each with a set of
images of peak, center, width, and integral would simplify the data
presentation.  One could even hope that the separate Gaussians fit
might have some separate physical reality.  The task could fit a
single Gaussian to each spatial position without much human guidance,
but detailed human interactivity was required for more complex
spectra.  Use of the task involved an endurance contest between the
human and the power company, requiring both luck and skill on the part
of the human. The task had no way to back up and fix fits generated at
previous spatial positions and either wrote out all of the fit
parameters in images or discarded all work done on the cube.  Clearly,
this was not a desirable situation and would-be users tended to prefer
fitting routines in other packages such as {\tt GYPSY}\@.

There has always been some interest in fitting rotation-measures to
polarization images taken at different frequencies.  Until the advent
of the Jansky VLA however, the images available for fitting rotation
measures were at a modest number of scattered frequencies and could
generally be fit by rather simple algorithms.  This is no longer the
case and better fitting algorithms were required.  Initially there was
considerable interest in a method known as Faraday Rotation Measure
Synthesis.\footnote{Brentjens, M. A. and de Bruyn, A. G. 2005, {\it
Astron.\ \& Ap.}, 1217.}  The method was implemented in \AIPS\ by
Leonid Kogan as task \keyw{FARS} with supplementary tasks and
additions by the author.  However, it was found that the ``synthesized
beam width'' in this method was too wide to separate most
scientifically interesting cases, even if the data were taken to
frequencies as low as 1 GHz.  Data from lower frequencies would be
needed to produce better component discrimination, but most radio
sources are completely depolarized at these low frequencies.  A more
direct fitting method is therefore required.

The author undertook a re-write of \keyw{XGAUS} \Hi{in 2013} in order to
make it a much more useful and user-friendly task and to provide a
software framework for a new task to do polarization fitting.  The
suggestion was then made to create a task to fit standard
Zeeman-splitting parameters as well.  The following Memo attempts to
describe in detail the function of all three of the tasks that have
resulted.  This memo has been revised for the {\tt 31DEC15} version of
these tasks and does not attempt to describe differences with the {\tt
  31DEC13} and {\tt 31DEC14} versions, in which these tasks first
appeared.  \Hi{The 2025 modifications have added useful new
  capabilities.}  Two tasks, {\tt AGAUS} and {\tt ZAMAN}, have appeared
in the {\tt 31DEC17} release of \AIPS\ to fit Gaussians and Zeeman
splitting in absorption spectra.  They are described in \AIPS\ Memo 122,
\Hi{which was also revised in 2025}.

Two tasks have been written to prepare plot files from the fits done by
these tasks.  There are also tasks in \AIPS\ which modify $uv$-data
files to add model components both spatially and spectrally.  In
addition there are tasks which add polarization and spectral models to
images.  These will be described in a final section.

\section{Gaussian fitting: {\tt XGAUS}}

Gaussian fitting is the process of determining the peak brightness,
center position, and width of a number of one-dimensional Gaussians,
plus a spectral baseline, which best matches a one-dimensional slice
of input data.  {\tt XGAUS} will do this for every voxel in a
three-dimensional image cube, taking each row as the one-dimensional
slice.  It then produces images over the second and third input axes
of the fit parameters.  Normal usage has the first axis of the cube as
a spectral axis, either in frequency or velocity units, and the second
and third axes are celestial coordinates.  This usage will be assumed
in the later discussion, but it is not required.  Mathematically, the
function fit is
\begin{equation}
  T(x,y,z) = A + Bx + \sum_1^{\tt NGAUSS} T_i e^{-\alpha (x - x_i)^2 /
     \sigma_i^2}
\end{equation}
where $x$ is the coordinate along the first axis of the image cube,
$A(y,z) + B(y,z)x$ is the spectral baseline, $T_i(y,z)$ are the peak
values of each of the {\tt NGAUSS} components, $x_i(y,z)$ are the
center channels of the components, $\sigma_i(y,z)$  are the full
widths at half maximum of the components and $\alpha = 4 \ln(2)$ to
impart this meaning to the $\sigma_i$.

\subsection{Inputs}

The usual {\tt INNAME} {\it et al.} adverbs define the cube to be fit
and the {\tt OUTNAME} {\it et al.} adverbs define the (eventual)
output name.  The {\tt INVERS} adverb controls which {\tt XG} table is
used by the task, with $\leq 0$ meaning a new table.  {\tt BLC} and
{\tt TRC} define the pixel ranges to be used in the current execution,
where {\tt BLC(1)} and {\tt TRC(1)} control the spectral channels that
will be fit and {\tt BLC(2)}, {\tt TRC(2)}, {\tt BLC(3)}, and {\tt
TRC(3)} control the area in celestial coordinates over which the
fitting is done.  New {\tt XG} table files are created for the
entire input image cube, but are filled with peak brightnesses limited
by {\tt BLC(1)} through {\tt TRC(1)}\@.  Therefore, it would be wise
to use as much of the first axis as is reliable when creating a new
table.  The second and third values of these adverbs then limit the
area over which fitting is done during this execution of {\tt
  XGAUS}\@.  Adverbs {\tt YINC} and {\tt ZINC} control the stride
taken in the first pass through the cube; a second pass will then fit
all voxels not fit in the first pass.  {\tt FLUX} controls which
spectra will be fit during this execution; all positions with 3
consecutive channels averaging above {\tt FLUX} will be fit.  The
initial guess for the linear baseline is now always zero.  Adverb
  {\tt ORDER} controls the order of the baseline with $<0$ meaning
  none, 0 meaning a constant, and $\ge 1$ meaning a constant plus a
  slope.  {\tt DOOUTPUT} controls what files are written --- this may
be changed interactively \Hi{(except for one-dimensional images)} so
leave it zero at this point.  Set \Hi{{\tt DOTV}$ > 0$} to use TV menus
to prompt you.  Even when fitting only 1 Gaussian component, it is best
to watch what is happening so you should never set this adverb false.
{\tt DORESID} controls whether residuals are plotted on the fit spectra;
such plots often provide clues when more Gaussians are needed for the
best fits.  In general it is best to set {\tt PIXRANGE} to zero to see
the full range of image values.  \Hi{In {\tt 31DEC25}, you may change
  {\tt PIXRANGE} interactively, but only when the task reaches the edit
  stage.  Cutting off the highest values in order to see weaker features
  better may} make it harder to set initial guesses however since you
can't point at the peak of a clipped component.  Set {\tt LTYPE} to your
favorite type of labeling, set {\tt PIXVAL} to zero to see all positions
being fit, and leave {\tt NITER} zero since 100 is more than enough
iterations. Set {\tt NGAUSS} to the number of Gaussians to be fit in
this execution.  It may be changed if you re-start on a pre-existing
{\tt XG} table since those tables contain room for the maximum number
of Gaussians allowed (8 at this writing).  {\tt RMSLIMIT} is an upper
limit for the rms of a fit before the fit is viewed as ``failed''
which causes the TV and interaction to be turned back on after you
have turned it off.  You should get a good idea of an appropriate
value from your initial uses of {\tt XGAUS} or from your knowledge of
the noise in your data cube.

The task begins by creating an {\tt XG} table and populating each row
with the largest average brightness over three consecutive channels in
the corresponding row in the range {\tt BLC(1)} through {\tt
TRC(1)}\@.  Then it reads the table every {\tt YINC} rows and {\tt
ZINC} planes and, for those with an average brightness greater than
{\tt FLUX}, attempts a fit.  Your interaction with this fit will be
described below.  After the first pass, the task loops over every row
and plane fitting those positions which have enough brightness and
which have not already been fit.  Finally, after all pixels above {\tt
  FLUX} have been fit, the task goes into an ``edit'' mode.  It
constructs images of each fit parameter and the integral flux of the
Gaussian and of the uncertainties in these parameters.  You may view
these images, select positions explicitly or by their parameter values
or rms and revisit the fits of those positions.  This stage will be
described in detail below.

At any time you may exit the task and then re-start it using the
same {\tt XG} table.  Good reasons for doing this include fitting
smaller regions with each pass using the appropriate number of
Gaussians for that region.  Doing small regions which will have
similar parameter values helps a great deal with the initial guessing
done by the task (mostly using the previous solution).  You might also
fit the cube initially with a high value of {\tt FLUX} and then
re-start with a lower value to extend the areas fit.

\subsection{Fitting}

\begin{figure}
\begin{center}
\resizebox{5.3in}{!}{\putfig{XGAUS.a}}
\centerline{\hphantom{MM}}
\resizebox{5.3in}{!}{\putfig{XGAUS.b}}
\caption{\Hi{Top: First model spectrum to fit, initial ``guess'' is all
    zero.  Bottom: better initial guess is requested.  Note change of
    color.}}
\label{fig:XGAUS.init}
\end{center}
\end{figure}

The fitting process starts with a plot of the spectrum with the data
and axis labels in graphics channel one (usually yellow) and the
initial guess as X's in graphics channel two (usually green).  Then
you are offered a menu of options, as shown in the present figures, on
the TV (\Hi{{\tt DOTV}$> 0$}).  The first spectrum to be fit (from a
model image cube for clarity) is illustrated in  the top of
Figure~\ref{fig:XGAUS.init} showing that the first guess for more than
one Gaussian is not useful.  The menu that appears at this point is\\

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{XGAUS.c}}
\caption{\Hi{First model spectrum to fit: better initial ``guess''
  has been entered.}}
\label{fig:XGAUS.guessed}
\end{center}
\end{figure}

\begin{center}
  \hicol
\begin{tabular}{|l|l|}\hline
   {\tt DO FIT}   & Attempt a fit \\
   {\tt RE-GUESS} & Interactively enter a new initial guess \\
   {\tt BAD}      & Mark pixel bad and go on \\
   {\tt QUIT}     & Exit {\tt XGAUS} now \\ \hline
\end{tabular}
\end{center}
\hblack

You select a menu option by moving the cursor to the desired option
with the mouse and registering that move with the TV by clicking the
left mouse button.  The selected menu item will change color as shown
in the figure.  If you press TV ``button'' {\tt D} at this point
(actually keyboard character D), helpful information about the
selected item will appear in your terminal window.  If you press one
of TV ``buttons'' {\tt A}, {\tt B}, or {\tt C} (actually keyboard
characters A, B, C), the selected function will be performed.  The
option to {\tt QUIT} causes the task to quit at this point.  You may
re-start later.  The option {\tt BAD} will mark this position as failed
and go on to the next position.  The option {\tt DO FIT} will cause the
task to attempt the non-linear fit with the current initial guess.  The
selected option in Figure~\ref{fig:XGAUS.init} is {\tt RE-GUESS} which
causes the task to change the color of the plot as illustrated in the
bottom of Figure~\ref{fig:XGAUS.init}.  {\tt XGAUS} will then prompt you
first to {\it ``Position cursor at center \&\ height of Gaussian
  component 1''}.  Move the cursor to the peak of component 1 and press
any TV button.  This selects the peak value and center of component 1.
Then the task prompts you to {\it ``Position cursor at half-width of
  Gaussian component 1''}.  Move the cursor horizontally to the
approximate position of the half-power point of component 1 and press
any TV button.  The horizontal position of the cursor then sets the
initial guess of the full width of the component.  These prompts are
then repeated for components two through {\tt NGAUSS}\@.  If you do not
want to fit a particular component at this position, move the cursor
outside the rectangular border line (\ie\ outside the data area of the
plot, but still inside the TV) before pressing the TV button for that
component; it will be omitted.

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{XGAUS.d}}
\caption{\Hi{First model spectrum to fit: good fit obtained.}}
\label{fig:XGAUS.good}
\end{center}
\end{figure}

After the improved initial guess has been entered, the plot is changed
to show the new guess, as illustrated in Figure~\ref{fig:XGAUS.guessed}.
This guess is good so the {\tt DO FIT} option is selected.  After TV
``buttons'' {\tt A}, {\tt B}, or {\tt C} are pressed the task attempts
the fit with the current initial guess.  Then the plot is changed with
the addition of the fit function in graphics plane 4 (usually cyan)
and the residual in graphics plane 7 (usually a dirty yellow).  This is
illustrated in Figure~\ref{fig:XGAUS.good}.  A different menu appears
at this point containing\\
\begin{center}
  \hicol
\begin{tabular}{|l|l|}\hline
   {\tt GOOD}     & Mark solution good and go on \\
   {\tt DO FIT}   & Attempt a new fit \\
   {\tt RE-GUESS} & Interactively enter a new initial guess \\
   {\tt TVOFF}    & Mark solution good and turn off interactivity \\
   {\tt HAND}     & Enter a solution from the terminal \\
   {\tt BAD}      & Mark pixel bad and go on \\
   {\tt 1}        & Interactively enter a 1-Gaussian initial guess \\
   {\tt 2}        & Interactively enter a 2-Gaussian initial guess \\
   {\tt 3}        & Interactively enter a 3-Gaussian initial guess \\
   {\tt QUIT}     & Exit {\tt XGAUS} now \\ \hline
\end{tabular}
\end{center}
\hblack

The option to {\tt QUIT} causes the task to quit at this point.  You may
re-start later.  The option {\tt BAD} will mark this position as failed
and go on to the next position.   The option {\tt RE-GUESS} will loop
back to prompt you for a new guess and repeat the fit.  Options {\tt 1},
{\tt 2}, $\ldots$, {\tt NGAUSS} will loop back to request an initial
guess with the selected number of Gaussians.  ({\tt NGAUSS=2} in the
current example.)  Option {\tt HAND} will prompt you to enter using the
terminal the Gaussian parameters for each component.  Enter on one line
for each component, the peak value of the Gaussian (in image units), the
center (in pixels with respect to the reference pixel), and width (in
pixels).  Appropriate ranges of parameters in these units can be seen
from the display of the current fit values.  You may also enter flags to
cause one or more parameter values to be fixed should you re-fit the
current spectrum.  The flags are entered after the 3 parameter values,
flags $\leq 0$ mean to hold the corresponding parameter fixed and
omitted flags are taken as 1.  {\tt XGAUS} will then repeat the display
in Figure~\ref{fig:XGAUS.good} to see if you made a good guess.
Immediately after a {\tt HAND} operation only, the option {\tt DO FIT}
is offered to go back with the hand-entered values as the initial guess
for a new fit.

Option {\tt TVOFF} allows you to turn off interactivity, allowing the
task to run using its own initial guesses until it finds a completely
unreasonable solution or one with an rms greater than {\tt RMSLIMIT}\@.
When that happens, you are shown the offending fit parameters \Hi{and
the task resumes interactivity.  Frequently this will drop you
immediately into entering a new guess as in the bottom of
Figure~\ref{fig:XGAUS.init} where the color change should make this
apparent.  Point to the Gaussian peak and half power following the
instructions on the terminal.  Sometimes it simply shows something like
the top part of Figure~\ref{fig:XGAUS.init} and you must select {\tt
  RE-GUESS} to try to fix things.}  Option {\tt GOOD} tells the task
that you are (reasonably) happy and that it should go on to the next
position.

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{XGAUShia}}
\caption{\Hi{First 1-Gaussian HI spectrum of NGC6503 to fit: good fit
    obtained automatically.}}
\label{fig:XGAUS.HIgood}
\end{center}
\end{figure}

When fitting only a single Gaussian, {\tt XGAUS} makes an initial
guess based mostly on finding a real peak in the spectrum.  This is
quite reliable, so turning off the TV interaction may save a great
deal of effort, although there will possibly be bad positions to be
fixed up in the next stage of this task.  An example of an HI cube, on
its first position, is shown in Figure~\ref{fig:XGAUS.HIgood}.
%\vfill\eject

\subsection{Editing and output}

Eventually all positions selected by {\tt FLUX} and {\tt BLC} and {\tt
  TRC} will have been fit.  At this point, the task computes images of
the fit parameters plus the ``flux'' (area under the Gaussian) and
their uncertainties.  \Hi{This ``edit'' stage is omitted if the input
image is one dimensional, but, in that case, channel flagging is
available in the initial fitting of the spectrum.}  It then offers a
lengthy menu of options which will allow you to view these images and
revisit positions that seem to have produced incorrect fits.  If {\tt
  NGAUSS} $> 1$, options to swap portions of image $n$ with
corresponding portions of image $m$ are also offered.  This ``edit''
menu is illustrated in Figure~\ref{fig:XGAUS.edit}.  Note that the size
of the signal portion of the model image is quite small.  When an
appropriate sub-image is selected for display, as was done for this
figure, {\tt XGAUS} replicates pixels in both directions to make the
image large enough to see.

There are three kinds of editing implemented here.  In the first, the
user establishes the parameter extrema which should be viewed as
acceptable.  The extrema currently set are shown in the title lines.
Then, {\tt XGAUS} may be told to flag all solutions not meeting these
criteria, or, more profitably perhaps, to revisit those positions to
see if a better fit can be obtained.  The other editing methods are
similar, but act on a list of pixel positions.  These may be entered
by typing in values, by clicking on suspect pixel positions in the
{\tt CURVALUE} function described below, \Hi{or by the {\tt BLOTCH LIST}
function also described below}.  The contents of the list may be viewed,
the solutions at the positions may be flagged, or they may be revisited
to attempt for a better solution, or the solutions at the listed
positions may be swapped between components $n$ and $m$.  The menu will
offer only appropriate swaps between components, thus 1 and 2 in our
{\tt NGAUSS=2} example.  Swapping may be required if {\tt XGAUS} gets
confused as to which component you want to call number 1 and which
number 2.  Such confusion appears in Figure~\ref{fig:XGAUS.edit}.  After
the flagging, revisiting, or swapping, the list is cleared.    The first
column of the menu shows the following:
%\vfill\eject

\begin{figure}
\begin{center}
\resizebox{5.1in}{!}{\putfig{XGAUS.e}}
\caption{\Hi{Model source, edit view}}
\label{fig:XGAUS.edit}
\end{center}
\end{figure}

\begin{center}
\begin{tabular}{|l|l|}\hline
 {\tt EXIT           } & Exit {\tt XGAUS}, writing output images if
                         {\tt DOOUTPUT} is now $> 0.$\\
 {\tt SET MIN S/N    } & Set minimum amplitude S/N(s) for okay
                         solutions\\
 {\tt SET MAX RES    } & Set maximum residual for okay solutions\\
 {\tt SET PEAK RANGE } & Set peak value range(s) for okay solutions\\
 {\tt SET CENTER RANGE } & Set X center range(s) for okay solutions\\
 {\tt SET WIDTH RANGE} & Set width range(s) for okay solutions\\
 {\tt SET MAX ERR WID} & Set maximum error(s) in width for okay
                         solutions\\
 {\tt REDO ALL       } & Re-do all solutions which are not okay
                         following the above criteria\\
 {\tt FLAG ALL       } & Mark bad all solutions which are not okay\\
 \Hi{{\tt SET PIXRANGE }} & \Hi{Enter new spectra plot range on terminal}\\
 {\tt OFF ZOOM       } & Turn of TV zoom\\
 {\tt OFF TRANSFER   } & Turn off black \&\ white and color TV
                         enhancements\\
 {\tt RESET WINDOW   } & Display full view of current image \\
 {\tt LABEL WEDGE?   } & Turn on/off labeling of step wedge \\
 {\tt SET DOOUTPUT   } & Increment {\tt DOOUTPUT} in loop 0-3 --- with
                         1 and 3 causing residual\\
                       & images and 2 and 3 causing parameter images
                         to be written on {\tt EXIT}\\
 {\tt ADD TO LIST    } & Type in pixel coordinates to add to list\\
 {\tt SHOW LIST      } & Display coordinates in list\\
 {\tt REDO LIST      } & Re-do solutions for all pixels in list\\
 {\tt FLAG LIST      } & Flag solutions for all pixels in list\\
 {\tt SWAP LIST 1-2  } & Swap solutions for components 1 and 2 for all
                         pixels in list\\ \hline
\end{tabular}
\end{center}
\vfill\eject

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{XGAUS.rc1}}
\caption{\Hi{Spectrum to be re-done after doing a {\tt RE-GUESS}.
    {\tt FLAG+FIT} selected for illustrative purposes.}}
\label{fig:XGAUS.redo}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\resizebox{5.08in}{!}{\putfig{XGAUS.rc2}}
\centerline{\hphantom{MM}}
\resizebox{5.08in}{!}{\putfig{XGAUS.rc3}}
\caption{\Hi{Top:Spectrum to be re-done in flag menu.  Some channels are
    flagged for illustrative (not scientific) purposes and the pink
    lines show another range of channels which could be flagged.
    Bottom: Acceptance menu after {\tt FLAG+FIT}.  Note that the good
    initial guess constrains the fitting despite the questionable
    channel flagging.}}
\label{fig:XGAUS.flag}
\end{center}
\end{figure}

\hicol
In {\tt 31DEC25}, {\tt XGAUS} offers an additional capability in the
edit stage.\footnote{\Hi{This option is available at the beginning when
    the input image is one-dimensional.}}  When doing {\tt REDO ALL} and
{\tt REDO LIST}, the menu shown in Figure~\ref{fig:XGAUS.init} and
Figure~\ref{fig:XGAUS.guessed} has an additional option {\tt FLAG+FIT}
as shown in Figure~\ref{fig:XGAUS.redo}.  When this option is selected a
new menu is shown as in Figure~\ref{fig:XGAUS.flag}:
\begin{center}
\begin{tabular}{|l|l|}\hline
  {\tt UNDO ALL}   & Undo all channel flagging \\
  {\tt FLAG CHANS} & Interactively flag a range of channels \\
  {\tt UNDO CHANS} & Interactively undo flags in a range of channels \\
  {\tt RE-GUESS}   & Interactively enter a new guess for the
                     Gaussians \\
  {\tt DO THE FIT} & Leave the flag routine and do the fit \\ \hline
\end{tabular}
\end{center}

Flagged channels are plotted in a different color than unflagged
channels (usually reddish and yellow, respectively).  {\tt UNDO ALL}
simply removes all flags in the current spectrum.  {\tt FLAG CHANS} and
{\tt UNDO CHANS} bring up two vertical lines (usually pink in color).
The cursor position controls the position of the lines, button A
switches between the two, button B sets a channel range and continues,
button C sets a range and exits, and button D exits without setting a
new range.  {\tt RE-GUESS} takes you to the screen to set a new guess
(Figure~\ref{fig:XGAUS.init}) and then returns to this flag menu.
Finally {\tt DO THE FIT} attempts the Gaussian fit and returns to the
menu to accept the result.  That screen is illustrated in the bottom of
Figure~\ref{fig:XGAUS.flag}.  Note the {\tt FLAG+REDO} option added to
this menu in the edit phase.  That option takes you back through the
sequence just illustrated.
\hblack
\vfill\eject

The second (and potentially third and even more) columns of the main
edit menu contain {\tt NGAUSS} sets of functions

\begin{center}
\begin{tabular}{|l|l|}\hline
 {\tt SHOW IMAGE A1 } & Enter image interaction with peak value of
                     component 1\\
 {\tt SHOW IMAGE C1 } & Enter image interaction with center pixel of
                     component 1\\
 {\tt SHOW IMAGE W1 } & Enter image interaction with width of
                     component 1\\
 {\tt SHOW IMAGE F1 } & Enter image interaction with "flux" of
                     component 1\\
 {\tt SHOW IMAGE EA1} & Enter image interaction with uncertainty in
                     peak value of component 1\\
 {\tt SHOW IMAGE EC1} & Enter image interaction with uncertainty in
                     center pixel of component 1\\
 {\tt SHOW IMAGE EW1} & Enter image interaction with uncertainty in
                     width of component 1\\
 {\tt SHOW IMAGE EF1} & Enter image interaction with uncertainty in
                     "flux" of component 1\\ \hline
\end{tabular}
\end{center}

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{XGAUShib}}
\caption{HI galaxy: peak value coloredby {\tt OFMCOLOR}; image
  interaction}
\label{fig:XGAUS.HIpeak}
\end{center}
\end{figure}

On very crowded menus, the word {\tt SHOW} may be omitted.  When you
select one of these functions most of the following operations will
appear in yet another menu.  This menu is illustrated in
Figure~\ref{fig:XGAUS.HIpeak}.  These options mostly invoke familiar
functions from \AIPS\ to control the {\tt FUNCTYPE} used in loading the
image to the display, to enhance the image intensities, to color the
enhanced image intensities, and to zoom the display.  Only one of the
{\tt LOAD AS} options will appear, with the next one in the sequence
offered when the current one has been invoked.  The {\tt SET WINDOW}
option allows you to select a sub-image to view in greater detail, while
{\tt RESET WINDOW} returns to viewing the full image.  The {\tt SWAP
  $n$-$m$} options will appear as needed when {\tt NGAUSS} $ > 1$ and
are used to change your selection of component number $n$ and component
$m$ in case {\tt XGAUS} got confused.  The {\tt NEXT WINDOW} option
appears when needed to display an image too large to fit on the display
screen.

\vfill\eject

\begin{center}
\begin{tabular}{|l|l|}\hline
 {\tt RETURN     } & Return to the above menus, image stays displayed\\
 {\tt LOAD AS SQ } & Re-load image with square root transfer function\\
 {\tt LOAD AS LG } & Re-load image with log transfer function\\
 {\tt LOAD AS L2 } & Re-load image with extreme log transfer function\\
 {\tt LOAD AS LN } & Re-load image with linear transfer function\\
 {\tt SET WINDOW}  & Set a sub-image to view\\
 {\tt RESET WINDOW} & Return too viewing the full image\\
 {\tt OFF TRANSF } & Turn off enhancement done with {\tt TVTRANSF}\\
 {\tt OFF COLOR  } & Turn off color enhancements\\
 {\tt TVTRANSF   } & Black \&\ white image enhancement\\
 {\tt TVPSEUDO   } & Color enhancement of numerous sorts\\
 {\tt TVPHLAME   } & Color enhancement of flame type, multiple colors\\
 \Hi{{\tt OFMCOLOR}} & \Hi{Many standard color tables (OFMs)
                     available} \\
 {\tt TVZOOM     } & Interactive zooming and centering of image\\
 {\tt CURVALUE   } & Display value under cursor, mark pixels for list\\
\Hi{{\tt BLOTCH LIST}} & \Hi{Interactively add a blotch area to the
                      pixel list} \\
 {\tt SWAP 1-2   } & Swap solutions for components 1 and 2
                      interactively\\
 {\tt NEXT WINDOW} & Move to next window in large images\\ \hline
\end{tabular}
\end{center}

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{XGAUShic}}
\caption{HI galaxy: center channel image interaction, colored by {\tt
    TVPSEUDO}.}
\label{fig:XGAUS.HIcent}
\end{center}
\end{figure}

Two operations in this menu are different from the usual.  {\tt
  CURVALUE} provides the capability of selecting positions for the
edit ``list.''  During the {\tt CURVALUE} operation, position the
cursor over the desired pixel and press buttons {\tt A}, {\tt B}, or
{\tt C} to add that pixel to the list.  The \Hi{{\tt BLOTCH LIST} and}
{\tt SWAP $n$-$m$} operations use a TV blotch operation like that in the
\AIPS\ verb {\tt TVSTAT} and task {\tt BLSUM}\@.  You are to mark with a
``blotch'' region those pixels in the present image which are to \Hi{be
added to the pixel list or} to have their solutions swapped with those
of the selected component.  Instructions will appear in the message
window as you proceed.  Begin by positioning the cursor at a pixel to
be the first vertex of a connected sequence of vertices and press TV
button {\tt A}\@.  Move to the next vertex and press button {\tt A}
again and repeat until you have marked all vertices for this region.
Then press button {\tt D} if you are done with this region or button
{\tt C} if you need to re-position one of the vertices.  In this case,
move the cursor to the vertex to be re-positioned, press button {\tt A}
and drag the vertex to the corrected position.  Press button {\tt A} or
{\tt B} to fix that vertex and go on to reset another vertex or {\tt D}
to end this region.  You may do this as many times as needed.

When you have finished getting the images just the way you want them,
you may write them out as \AIPS\ image files.  Select the {\tt SET
  DOOUTPUT} option until its value, shown at the top of the screen, is
what you want.  In {\tt XGAUS}, values 1 and 3 cause a residual image
cube to be written, while values 2 and 3 cause images of the parameter
values and their uncertainties to be written.  The baseline and slope
images and their uncertainties are given class {\tt CONST}, {\tt
  SLOPE}, {\tt DCONST}, and {\tt DSLOPE}, while the parameter images
and their uncertainties get class {\tt AMPL$n$}, {\tt CENTR$n$}, {\tt
  WIDTH$n$}, {\tt DAMPL$n$}, {\tt DCENT$n$}, and {\tt DWIDT$n$}, and
the flux and its uncertainty get class {\tt FLUX$n$}, and {\tt
  DFLUX$n$}, where $n$ is the component number.

%\vfill\eject
\section{Zeeman splitting: {\tt ZEMAN}}

When an intrinsically unpolarized spectral line is emitted in the
presence of a magnetic field, the right and left circular
polarizations have their frequencies shifted in opposite directions by
an amount proportional to the magnetic field (at least for modest
magnetic fields).  The traditional analysis of data to measure this
splitting works only for those cases in which the separation of
polarizations is a small fraction of the line width.  In that case,
the function that has been traditional is
\begin{equation}
     V(x,y,z) = A(y,z) I(x,y,z) + 0.5 B(y,z)\,\, \frac{dI(x,y,z)}{dx}
\label{eqn:Zeman}
\end{equation}
where $V(x,y,z)$ is the V Stokes polarization component. $I(x,y,z)$ is
the unpolarized I Stokes component, $x$ is the spectral axis value,
$(y,z)$ is the celestial coordinate value, and $A(y,z)$ and $B(y,z)$
are the parameters to be found by a linear least-squares
method.\footnote{Sault, R.J., Killeen, N. E. B.., Zmuidzinas, J.,
  Loushin, R. 1990, {\it Ap.~J.}, {\bf 74}, 437-461.}  The \AIPS\ task
{\tt ZEMAN} offers this model, with two choices for the method by
which the derivative of $I$ with $x$ is computed.  If the total
intensity spectrum has had {\tt XGAUS} applied, another function may
be used instead:
\begin{equation}
     V(x,y,z) = A(y,z) I(x,y,z) + 0.5 \sum_1^{\tt NGAUSS} B_i(y,z)\,
        \,\frac{dG_i(x,y,z)}{dx}
\label{eqn:ZemanGaus}
\end{equation}
where $G_i(x,y,z)$ is the spectrum of the $i$'th Gaussian component
and one solves for $A(y,z)$ and {\tt NGAUSS} component frequency
separations $B_i(y,z)$.  \Hi{In an early study, VLA data in the C
configuration lacked the spatial resolution to separate two source
components that emit OH lines at different frequencies.  This second
method easily found two different magnetic field values.  This source
has been observed with much better spatial resolution in A configuration
which separates the two components spatially.  The field values found by
{\tt ZEMAN} matched those found in the published, high-resolution data.}
Note that {\tt ZEMAN} returns $B$ in units matching the input $x$ axis
(usually Hz or m/sec).  Any association with magnetic field values is
left to the user.

\subsection{Inputs}

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{ZEMAN.a}}
\caption{\Hi{First model spectrum to fit, initial ``guess'' has all {\tt
    ZEMAN} parameters zero.}}
\label{fig:ZEMAN.init}
\end{center}
\end{figure}

The inputs for {\tt ZEMAN} are very similar to those for {\tt
  XGAUS}\@.  {\tt INNAME} {\it et al.}~specify the V polarization
cube which must be in transposed form with frequency (or velocity) as
the first axis.  {\tt IN2NAME} {\it et al.}~specify the corresponding
I polarization cube, similarly transposed.  The axes of the two images
must match.  {\tt BLC} and {\tt TRC} define the spectral and celestial
coordinate regions of interest for the fit and {\tt FLUX} gives
the lower limit in the I image for the average of three consecutive
channels for the spectra to be fit.  {\tt INVERS} specifies the input
version of the {\tt ZE} table in which the {\tt ZEMAN} results are
stored.  Zero means to make a new one, other values mean to re-visit an
existing solution.  Adverb {\tt OPTYPE} specifies which of the above
formul\ae\ are solved.  {\tt OPTYPE = 'GAUS'} says to use the {\tt XG}
file {\tt IN2VERS} attached to the input I polarization cube to solve
using Equation~\ref{eqn:ZemanGaus}.  {\tt OPTYPE = '2SID'} says to solve
Equation~\ref{eqn:Zeman}, evaluating the derivative by
$$
\frac{dI(x,y,z)}{dx} = 0.5 \left[\, I(x+1,y,z) - I(x-1,y,z) \,\right]
$$
while any other {\tt OPTYPE} value says to evaluate the derivative
with
$$
\frac{dI(x,y,z)}{dx} = I(x,y,z) - I(x-1,y,z)
$$
while solving Equation~\ref{eqn:Zeman}.  {\tt DOOUTPUT} controls what
files are written --- this may be changed interactively \Hi{(except for
  one-dimensional images)} so leave it zero at this point.  Set \Hi{{\tt
  DOTV}$> 0$} to use TV menus to prompt you.  Even though the fitting
operation is linear and so needs little or no guidance, it is best to
watch what is happening so you should never set this adverb false.  If
all seems well, you can turn off the TV after watching a few of the
solutions.  $A(y,z)$ is always fit but should be very close to zero if
you are fitting a cube which has already had the leakage term fit and
removed.  {\tt RMSLIMIT} is an upper limit for the rms of a fit before
the fit is viewed as ``failed'' which causes the TV and interaction to
be turned back on after you have turned it off.  You should get a good
idea of an appropriate value from your initial uses of {\tt ZEMAN} or
from your knowledge of the noise in your data cube.
\vfill\eject

\subsection{Fitting}

The fitting process starts with a plot of the I polarization spectrum
across the top and the V polarization across the bottom.  The data
and axis labels are plotted in graphics channel one (usually yellow)
and the initial guess as lines in graphics channel two (usually
green).  The Gaussians, when {\tt OPTYPE='GAUS'}, are shown as a
smooth line on the I polarization plot, while the initial guess,
plotted in the V polarization, is plotted at the locations of the data
samples.  Then you are offered a menu of options, as shown in the
present figures, on the TV (\Hi{{\tt DOTV} $> 0$}).  The first spectrum
to be fit (\Hi{from the model source used for {\tt XGAUS} in this memo})
is illustrated in Figure~\ref{fig:ZEMAN.init} showing that the first
guess \Hi{is zero for all {\tt ZEMAN} parameters}.  The menu that
appears at this point is\\
\begin{center}
  \hicol
\begin{tabular}{|l|l|}\hline
   {\tt DO FIT}   & To do the fitting \\
   {\tt BAD}      & To mark this voxel bad and go on to the next \\
   {\tt QUIT}     & To exit {\tt ZAMAN} now \\ \hline
\end{tabular}
\end{center}
\hblack

\begin{figure}
\begin{center}
\resizebox{5.5in}{!}{\putfig{ZEMAN.b}}
\caption{First model spectrum to fit, plot after fitting.}
\label{fig:ZEMAN.good}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\resizebox{5.5in}{!}{\putfig{ZEMAN.c}}
\caption{\Hi{Later model spectrum to fit, initial ``guess'' in green
  has Gaussians and values for $B_i(y,z)$, cyan shows actual fit.}}
\label{fig:ZEMAN.init2}
\end{center}
\end{figure}

Note that there is no {\tt RE-GUESS} since there is no need to enter
them with linear least squares.  You select menu options in the same
way as {\tt XGAUS}, with ``button'' {\tt D} giving appropriate
real-time help information and buttons {\tt A}, {\tt B}, and {\tt C}
selecting the highlighted option.  The option to {\tt QUIT} causes the
task to quit at this point.  You may re-start later.  The option {\tt
  BAD} will mark this position as failed and go on to the next position.
The option {\tt DO FIT}, currently highlighted, will cause the task to
attempt the linear fit with the current initial guess.  The result of
that fit is shown in Figure~\ref{fig:ZEMAN.good}.
\vfill\eject

\hicol The menu after each fit contains\\
\begin{center}
\begin{tabular}{|l|l|}\hline
   {\tt GOOD}   & To accept the solution and go on \\
   {\tt TVOFF}  & To accept the solution and go on with no
                  interaction \\
   {\tt HAND}   & To enter a new solution from the terminal \\
   {\tt BAD}    & To mark this solution as bad and go on  \\
   {\tt QUIT}   & To exit {\tt ZAMAN} now \\ \hline
\end{tabular}
\end{center}
What most of these do is obvious.  {\tt HAND} prompts you to enter the
Gain value and then either {\it Enter avg Field in channels} or {\it
  Enter Gaussian xx Field in channels} for {\it xx} $= 1$ to {\tt
  NGAUSS} depending on the {\tt OPTYPE} selected.  When you have entered
all values, {\tt ZEMAN} repeats the plot of Figure~\ref{fig:ZEMAN.init2}
with the new values.  When the task goes on to the next voxel it starts
with the $B_i(y,z)$ of the last fit (Figure~\ref{fig:ZEMAN.init2}) and
so has a much better initial guess.
\hblack

\subsection{Editing and output}

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{ZEMAN.d}}
\caption{\Hi{{\tt ZEMAN} fit (component 1 splitting ``field'') and image
  editing screen.  Note bad pixel due to bad pixel in {\tt XGAUS} fits;
  see Figure~\ref{fig:XGAUS.edit}.}}
\label{fig:ZEMAN.edit}
\end{center}
\end{figure}

Eventually all positions selected by {\tt FLUX} and {\tt BLC} and {\tt
  TRC} will have been fit.  At this point, the task computes images of
the fit parameters $A(y,z)$ and $B_i(y,z)$ and their uncertainties.
\Hi{This ``edit'' stage is omitted if the input image is one
dimensional, but, in that case, channel flagging is available in the
initial fitting of the spectrum.}  It then offers a menu of options
which will allow you to view these images and revisit positions that
seem to have produced incorrect fits.  This ``edit'' menu is illustrated
in Figure~\ref{fig:ZEMAN.edit}.  Note that the size of the region of
interest is quite small as illustrated here.  If {\tt SET WINDOW} is
used to select a small sub-image, then {\tt ZEMAN} will blow it up by
pixel replication to a reasonable size as in the following figure.  The
first column of the menu includes:

\begin{center}
\begin{tabular}{|l|l|}\hline
 {\tt EXIT           } & Exit {\tt XGAUS}, writing output images if
                         {\tt DOOUTPUT} is now $> 0.$\\
 {\tt SET MAX RES    } & Set maximum residual for okay solutions\\
 {\tt SET GAIN RANGE } & Set gain value range(s) for okay solutions\\
 {\tt SET FIELD RANGE } & Set field range(s) for okay solutions\\
 {\tt SET MAX ERR FLD} & Set maximum error(s) in field for okay
                         solutions\\
 {\tt REDO ALL       } & Re-do all solutions which are not okay\\
 {\tt FLAG ALL       } & Mark bad all solutions which are not okay\\
 {\tt OFF ZOOM       } & Turn of TV zoom\\
 {\tt OFF TRANSFER   } & Turn off black \&\ white and color TV
                         enhancements\\
 {\tt RESET WINDOW   } & Display full view of current image \\
 {\tt LABEL WEDGE?   } & Turn on/off labeling of step wedge \\
 {\tt SET DOOUTPUT   } & Increment {\tt DOOUTPUT} in loop 0 -- 3\\
 {\tt ADD TO LIST    } & Type in output pixel coordinates to add to
                         list\\
 {\tt SHOW LIST      } & Display coordinates in list\\
 {\tt REDO LIST      } & Re-do solutions for all pixels in list\\
 {\tt FLAG LIST      } & Flag solutions for all pixels in list\\ \hline
\end{tabular}
\end{center}

There are two kinds of editing implemented here.  In the first, the
user establishes the parameter extrema which should be viewed as
acceptable.  The extrema currently set are shown in the title lines
and include the maximum rms, the range of allowed values for ``gain''
($A(y,z)$), the range of allowed values for ``field'' ($B(y,z)$ in
pixels), and the maximum uncertainty in the field.  Then, {\tt ZEMAN}
may be told to flag all solutions not meeting these criteria, or, more
profitably perhaps, to revisit those positions to see why a poor fit
was obtained.  Note that, unlike the other tasks in this memo, the
linear nature of the fit in {\tt ZEMAN} means that only one solution
is possible at each celestial coordinate \Hi{unless the data are changed
  in some way}.  The other editing method is similar, but acts on a list
of pixel positions.  These may be entered by typing in values, by
clicking on suspect pixel positions in the {\tt CURVALUE} function
described below\Hi{, or by using the {\tt BLOTCH LIST} function also
  described below}.  The contents of the list may be viewed, the
solutions at the positions may be flagged, and they may be revisited to
see why they are suspect.  Note that there are no swapping of solutions
in this task; {\tt XGAUS} establishes which component is which.  After
the flagging or revisiting, the list is cleared.

\begin{figure}
\begin{center}
\resizebox{5.5in}{!}{\putfig{ZEMAN.e}}
\caption{\Hi{{\tt ZEMAN} initial display to re-do a fit.}}
\label{fig:ZEMAN.redo1}
\end{center}
\end{figure}

\hicol
In {\tt 31DEC25} the operations of {\tt REDO ALL} and {\tt REDO LIST}
offer an additional menu option.\footnote{\Hi{This option is available
at the beginning when the input image is one-dimensional.}}  When a
spectrum to be re-done is first displayed it is called {\tt FLAG
  \&\ FIT} as illustrated in Figure~\ref{fig:ZEMAN.redo1}.  If this
option is selected the screen shown in the upper portion of
Figure~\ref{fig:ZEMAN.redo2} appears.  The menu contains\\
\begin{center}
\begin{tabular}{|l|l|}\hline
  {\tt UNDO ALL}   & Undo all channel flagging \\
  {\tt FLAG CHANS} & Interactively flag a range of channels \\
  {\tt UNDO CHANS} & Interactively undo flags in a range of channels \\
  {\tt DO THE FIT} & Leave the flag routine and do the fit \\ \hline
\end{tabular}
\end{center}

\begin{figure}
\begin{center}
\resizebox{6.3in}{!}{\putfig{ZEMAN.f}}
\centerline{\hphantom{MM}}
\resizebox{6.3in}{!}{\putfig{ZEMAN.g}}
\caption{\Hi{Top: Spectrum to be re-done in flag menu.  Some channels on
    the right hand side are flagged and the pink lines show another
    range of channels which will also be flagged.  Bottom: Acceptance
    menu after {\tt FLAG \&\ FIT}.}}
\label{fig:ZEMAN.redo2}
\end{center}
\end{figure}

Flagged channels are plotted in a different color than unflagged
channels (usually reddish and yellow, respectively).  {\tt UNDO ALL}
simply removes all flags in the current spectrum.  {\tt FLAG CHANS} and
{\tt UNDO CHANS} bring up two vertical lines (usually pink in color).
The cursor position controls the position of the lines, button A
switches between the two, button B sets a channel range and continues,
button C sets a range and exits, and button D exits without setting a
new range.  Finally {\tt DO THE FIT} attempts the Zeeman fit and returns
to the menu to accept the result.  That screen is illustrated in the
bottom of Figure~\ref{fig:ZEMAN.redo2}.  Note the {\tt FLAG+REDO} option
added to this menu in the edit phase.  That option takes you back
through the sequence just illustrated.
\hblack

\vfill\eject

The second column of the main edit menu contains
\begin{center}
\begin{tabular}{|l|l|}\hline
 {\tt SHOW IMAGE G } & Enter image interaction with gain\\
 {\tt SHOW IMAGE EG } & Enter image interaction with uncertainty in
                     the gain\\
 {\tt SHOW IMAGE F1 } & Enter image interaction with field of
                     component 1\\
 {\tt SHOW IMAGE EF1} & Enter image interaction with uncertainty in
                     field of component 1\\
 {\tt SHOW IMAGE F2} & Enter image interaction with field of component
                     2\\
 {\tt SHOW IMAGE EF2} & Enter image interaction with uncertainty in
                     field of component 2\\ \hline
\end{tabular}
\end{center}
There is one or, if {\tt OPTYPE='GAUS'} as in the illustrated cases,
the maximum of {\tt NGAUSS} of the {\tt F}$n$ and {\tt EF}$n$ options.
When you select one of the {\tt SHOW} options, the options in yet
another menu appear along with a display of the selected image.  As
illustrated in Figure~\ref{fig:ZEMAN.field1}, these options are

\begin{figure}
\begin{center}
\resizebox{5.1in}{!}{\putfig{ZEMAN.h}}
\caption{{\tt ZEMAN} image of component 1 field $B_1(y,z)$ with display
  options.}
\label{fig:ZEMAN.field1}
\end{center}
\end{figure}

\begin{center}
\begin{tabular}{|l|l|}\hline
 {\tt RETURN     } & Return to the above menus, image stays displayed\\
 {\tt LOAD AS SQ } & Re-load image with square root transfer function\\
 {\tt LOAD AS LG } & Re-load image with log transfer function\\
 {\tt LOAD AS L2 } & Re-load image with extreme log transfer function\\
 {\tt LOAD AS LN } & Re-load image with linear transfer function\\
 {\tt SET WINDOW } & Set a sub-image to view\\
 {\tt RESET WINDOW} & Return to viewing the full image\\
 {\tt OFF TRANSF } & Turn off enhancement done with {\tt TVTRANSF}\\
 {\tt OFF COLOR  } & Turn off color enhancements\\
 {\tt TVTRANSF   } & Black \&\ white image enhancement\\
 {\tt TVPSEUDO   } & Color enhancement of numerous sorts\\
 {\tt TVPHLAME   } & Color enhancement of flame type, multiple colors\\
 \Hi{{\tt OFMCOLOR}} & \Hi{Many standard color tables (OFMs)
                     available} \\
 {\tt TVZOOM     } & Interactive zooming and centering of image\\
 {\tt CURVALUE   } & Display value under cursor, mark pixels for list\\
\Hi{{\tt BLOTCH LIST}} & \Hi{Interactively add a blotch area to the
                      pixel list} \\
 {\tt NEXT WINDOW} & Move to next window in large images\\ \hline
\end{tabular}
\end{center}
Only one of the {\tt LOAD AS} options will appear, namely the next
after the current transfer function from the list of linear, square
root, log, and more extreme log transfer functions.  {\tt CURVALUE}
provides the capability of selecting positions for the edit ``list.''
During the {\tt CURVALUE} operation position the cursor over the
desired pixel and press buttons {\tt A}, {\tt B}, or {\tt C} to add
that pixel to the list.  \Hi{{\tt BLOTCH LIST} is described in detail
in the {\tt XGAUS} and {\tt RMFIT} sections of this Memo.}  The other
options are familiar as {\tt AIPS} verbs.  Instructions for interaction
will appear on the terminal and button {\tt D} in the menu may be used
to obtain a helpful display on the terminal.

When you have finished getting the images just the way you want them,
you may write them out as \AIPS\ image files.  Select the {\tt SET
  DOOUTPUT} option until its value, shown at the top of the screen, is
what you want.  In {\tt ZEMAN}, values 4, 5, 6, and 7 cause a residual
(data-model) image cube to be written, while values 2, 3, 6, and 7
cause images of the gain and field values and their uncertainties to
be written and values 1, 3, 5, and 7 cause an image of the V cube to
be written with the gain times the I-polarization cube subtracted.
The corrected V-polarization image gets the class specified by {\tt
  OUTCLASS}, the residual image gets class {\tt VRESID}, the gain and
its uncertainty get classes {\tt GAIN} and {\tt DGAIN}, and the
field(s) and their uncertainties get classes {\tt FIELD}$n$ and {\tt
  DFELD}$n$, with $n = 1$ to the maximum {\tt NGAUSS}\@.
%\vfill\eject

\section{Polarization fitting: {\tt RMFIT}}

Faraday rotation occurs when a polarized source of radiation
$F_0(\Phi_1)$ has its radiation pass through a region of electrons (a
``Faraday screen'').  The complex polarization ($Q + iU$) is rotated
in position angle by an amount proportional to the wavelength squared,
\ie\ it is multiplied by $\exp (i 2 \phi_1 \lambda^2)$.  If this
radiation is then joined by another source of radiation $F_1(\Phi_2)$
and the two then pass through additional screens $\phi_2$, $\phi_3$,
and $\phi_4$ after which there is a third source of radiation and one
final screen, the net polarization observed  is then
\begin{eqnarray}
P(\lambda^2) & = & F_0(\Phi_1)\, e^{i2\Phi_1\lambda^2} + F_1(\Phi_2)\,
                   e^{i2\Phi_2\lambda^2} + F_2(\Phi_3)\,
                   e^{i2\Phi_3\lambda^2},\\
\noalign{\mbox{where}}
\Phi_1 & = & \phi_1 + \phi_2 + \phi_3 + \phi_4 + \phi_5 \nonumber\\
\Phi_2 & = & \phi_2 + \phi_3 + \phi_4 + \phi_5 \nonumber\\
\Phi_3 & = & \phi_5\nonumber\\
\noalign{\mbox{or, in the limit:}}
P(\lambda^2) & = & \int F(\Phi)\, e^{i2\Phi\lambda^2} d\Phi \label{eqn:FARS}
\end{eqnarray}

``Faraday Rotation Measure Synthesis'' is then the Fourier analysis of
Equation~\ref{eqn:FARS} as initially described by Brentjens and de
Bruyn.\footnote{Brentjens, M. A. and de Bruyn, A. G. 2005, {\it
    Astron.\ \& Ap.}, 1217.}  \AIPS\ task {\tt FARS} takes as input
two matching spectral cubes, one in Q and one in U polarization.  It
does a suitable Fourier transform in wavelength squared and can do a
one-dimensional complex Clean and restoration at each spatial
coordinate.  A sample of a two-component model is shown in
Figure~\ref{fig:FARS model}.  Although the model contains data from
frequencies from 1 GHz to 3 GHz, the figure shows that it is not
possible with this analysis to separate the two components even though
they are at a physically realistic and interesting separation.

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{FARS}}
\caption{{\tt FARS} output with model data.  Source has components at
100 and 140 radians per $m^2$.  Blue is without Clean, green is Clean
restored with very narrow components, and pink is Clean with a normal
restoration.  Note that Clean finds components at 120 plus about 68
and 172 radians per $m^2$.}
\label{fig:FARS model}
\end{center}
\end{figure}

Therefore a new fitting task, called {\tt RMFIT}, was written.  Its
job is to do a non-linear fit of
\begin{eqnarray}
  Q(i) & = & \sum  P_j \cos (2\theta_j + 2 R_j \lambda(i)^2)
              F(\beta_j,, \lambda(i)^2) \\
  U(i) & = & \sum  P_j \sin (2\theta_j + 2 R_j \lambda(i)^2)
              F(\beta_j,, \lambda(i)^2)
\end{eqnarray}
where $P_j$ is the polarized flux of component $j$, $\theta_j$ is the
polarization angle of component j, $R_j$ is the rotation measure of
component j, $\lambda(i)$ is the wavelength of spectral channel $i$,
$\beta_j$ is a ``thickness'' parameter and $F$ is one of
\begin{eqnarray*}
F(\beta_j,\lambda(i)^2) & = & \left(
    \frac{\lambda_1^2}{\lambda(i)^2}\right) ^ {\frac{\beta_j}{2}} \\
  & = & \frac{\sin (\beta_j \lambda(i)^2)}{\beta_j \lambda(i)^2} \\
  & = & e ^ {-\ln(2)\, \left(\,\frac{\beta_j \lambda(i)^2}{1.8954}\,\right)^2} \\
  & = & e ^ {-\ln(2)\, \frac{\beta_j \lambda(i)^2}{1.8954}}
\end{eqnarray*}
where $\lambda_1$ is the wavelength at the fiducial frequency of 1 GHz
and the thickness models will be referred to by the names spectral
index, slab, Gaussian, and exponential, respectively.  The peculiar
number $1.8954$ was chosen so that the last three models reach the 0.5
point at the same value of $\beta_j$.  Rotation measure thickness, in
general, makes the polarized signal decrease at longer wavelengths.
It is one of the reasons that radio sources at low frequency tend to
have little or no polarized flux density.

\subsection{Inputs}

{\tt RMFIT} expects five input cubes which contain, in order, the Q
polarization, the U polarization, the I polarization, the {\tt FARS}
output real or amplitude, and the {\tt FARS} output imaginary or phase
images.  All five must align appropriately, the first three must be
transposed so that frequency (or frequency ID number) is the first
axis (as they go into {\tt FARS}), and the last two have the rotation
measure axis first (as they come out of {\tt FARS})\@.  The {\tt FARS}
images are required to allow the user to make initial guesses for the
model.  Do some Cleaning in {\tt FARS} and restore with a ``Clean
beam'' reasonably close to the default.  The I polarization image is
optional and may be a continuum image rather than a spectral cube.

The rest of the inputs to {\tt RMFIT} are much like those of {\tt
  XGAUS}\@.  {\tt INVERS} defines the version of the {\tt RM} table to
be used to hold the model fit, with zero causing a new table to be
used.  {\tt BLC} and {\tt TRC} define the pixel ranges to be used in
the current execution, where {\tt BLC(1)} and {\tt TRC(1)} control the
spectral channels that will be fit and {\tt BLC(2)}, {\tt TRC(2)},
{\tt BLC(3)}, and {\tt TRC(3)} control the area in celestial
coordinates over which the fitting is done.  {\tt BLC(1)} and {\tt
TRC(1)} do not apply to the {\tt FARS} images, but they do control
which channels are used to find the peak polarization and unpolarized
brightnesses for new {\tt RM} tables.  The {\tt RM} table now contains
all pixels on axes 2 and 3 of the input images, so there should never
be a problem resuming fitting with a new region of interest.  Adverbs
{\tt YINC} and {\tt ZINC} control the stride taken in the first pass
through the cube; a second pass will then fit all pixels not fit in
the first pass.  Adverbs {\tt PCUT} and {\tt ICUT} control which
pixels are fit.  Any spectra having average total polarization
($\sqrt{Q^2 + U^2}$) greater than {\tt PCUT} and, if an I image is
given, total intensity greater than {\tt ICUT} will be fit.

{\tt DOOUTPUT} controls what files are written --- this may be
changed interactively so leave it zero at this point.  \Hi{{\tt DOTV}
is {\it required} to be true} and will use TV menus to prompt you.  The
badly named adverb {\tt NGAUSS} controls the maximum number of
components to fit and is limited to 4.  Adverb {\tt DOSPIX} controls the
choice of thickness model to be fit, with 1 for spectral index, 2 for
slab, 3 for Gaussian, 4 for  exponential, and anything else for no
thickness.  Polarization cubes may have quite different noise levels in
different spectral channels.   Adverb {\tt DOWEIGHT} controls whether a
weighted fit is done.  The weights may be found by robust rms methods
from the Q and U images or read from a text file specified by adverb
{\tt INFILE}\@.  {\tt RMSLIMIT} is an upper limit for the rms of a fit
before the fit is viewed as ``failed'' which causes the TV and
interaction to be turned back on after you have turned it off.  You
should get a good idea of an appropriate value from your initial uses of
{\tt RMFIT} or from your knowledge of the noise in your data cube.

A recommended strategy, when creating an {\tt RM} table, would be to
set the the spectral region {\tt BLC(1)} through {\tt TRC(1)} to
encompass all reliable spectral channels so that the full range of
polarized and unpolarized brightnesses will be found.  The task begins
by creating this table, including the entire input image, and populating
each row with the average unpolarized and polarized brightnesses in each
spectrum.  {\tt RMFIT} reads the table every {\tt YINC} rows and {\tt
  ZINC} planes and, for those with brightnesses exceeding {\tt PCUT} and
{\tt ICUT}, attempts a fit.  Your interaction with this fit will be
described below.  After the first pass, the task loops over every row
and plane fitting those positions which have enough brightness and which
have not already been fit.  Finally, after all pixels above the cutoffs
have been fit, the task goes into an ``edit'' mode.  It constructs
images of each fit parameter and Q and U at zero wavelength and of the
uncertainties in these parameters.  You may view these images, select
positions explicitly or by their parameter values or rms and revisit the
fits of those positions.  This stage will be described in detail below.

At any time you may exit the task and then re-start it using the
same {\tt RM} table.  Good reasons for doing this include fitting
smaller regions with each pass using the appropriate number of
components for that region.  Doing small regions which will have
similar parameter values helps a great deal with the initial guessing
done by the task (mostly using the previous solution).  You might also
fit the cube initially with a high value of {\tt PCUT} and then
re-start with a lower value to extend the areas fit.

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{RMFIT.a}}
\caption{Model field first pixel {\tt FARS} rotation measure spectrum
  to fit.  With {\tt NGAUSS=1} the initial guess is pretty good, but it
  is clear that this position has 2 components.}
\label{fig:RMFIT.init}
\end{center}
\end{figure}

\subsection{Fitting}

The fitting process starts with a plot of the {\tt FARS} amplitude and
phase as a function of rotation measure at the first coordinate with
adequate signal.  The data and axis labels are plotted in graphics
channel one (usually yellow) and the initial guess as X's in graphics
channel two (usually green).  Then you are offered a menu of options,
as shown in the present figures, on the TV.  The first spectrum to be
fit is illustrated in Figure~\ref{fig:RMFIT.init} showing that the first
guess for one component is sensible.  The menu that appears at this
point is\\

\begin{center}
  \hicol
\begin{tabular}{|l|l|}\hline
   {\tt DO FIT}   & {\tt RMFIT} will attempt the fit \\
   {\tt RE-GUESS} & {\tt RMFIT} will prompt you to point out the
                    RM values of the components\\
   {\tt BAD}      & {\tt RMFIT} will mark this fit as failed and
                    go on\\
   {\tt QUIT}     & {\tt RMFIT} will exit now \\ \hline
\end{tabular}
\end{center}
\hblack

You select a menu option by moving the cursor to the desired option
with the mouse and registering that move with the TV by clicking the
left mouse button.  The selected menu item will change color as shown
in the figure.  If you press TV ``button'' {\tt D} at this point
(actually keyboard character D), helpful information about the
selected item will appear on your terminal window.  If you press one
of TV ``buttons'' {\tt A}, {\tt B}, or {\tt C} (actually keyboard
characters A, B, C), the selected function will be performed.  The
option to {\tt QUIT} causes the task to quit at this point.  You may
re-start later.  The option {\tt BAD} will mark this position as failed
and go on to the next position.  The option {\tt DO FIT} will cause the
task to attempt the non-linear fit with the current initial guess.  The
selected option in Figure~\ref{fig:RMFIT.init} \Hi{is {\tt DO FIT}
  since, for one component, the guess is as good as it can be.}

After a fit has been attempted, the U spectrum is plotted in the top
half of the screen and the Q spectrum is plotted in the lower half.
If you have selected {\tt SHOW A\&P} previously, then the polarization
angle spectrum is plotted in the top half and the total polarization
spectrum is plotted in the lower half.  The data and labels are in
graphics channel 1 (usually yellow), the initial guess is plotted in
channel 2 (usually green), and the fit model is plotted in channel 4
(usually cyan).  The fit values and rmses are displayed on your
terminal.  The user chose to fit only \Hi{one} component initially,
getting the unhappy result shown in Figure~\ref{fig:RMFIT.1bad}.

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{RMFIT.b}}
\caption{\Hi{Model field first pixel Q and U spectra, fit with only one
  component.  Initial guess (green) and fit (cyan) are close but miss
  the data (yellow) badly.  User selects {\tt QUIT} to restart with 2
  components.}}
\label{fig:RMFIT.1bad}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\resizebox{6.3in}{!}{\putfig{RMFIT.c}}
\centerline{\hphantom{MM}}
\resizebox{6.3in}{!}{\putfig{RMFIT.d}}
\caption{\Hi{Top: RM spectrum with {\tt NGAUSS = 2}. Initial guess is not
    good so {\tt RE-GUESS} is selected.  Bottom: Display when entering a
    new initial guess is in a different color.}}
\label{fig:RMFIT.redo2}
\end{center}
\end{figure}

\hicol
The user, having found that there are actually two components, quit and
then restarted {\tt RMFIT} setting {\tt NGAUSS=2} and {\tt INVERS} to
point at the new {\tt RM} file.  Figure~\ref{fig:RMFIT.redo2} shows the
initial RM spectrum over again, but with two components the automatic
initial guess is not useful.  Selecting {\tt RE-GUESS}, the display in
the bottom of Figure~\ref{fig:RMFIT.redo2} appears in a different color
to signify that an initial guess is requested.  Position the cursor at
the peak of each component and press any TV button.  The resulting fit
is much better as seen in Figure~\ref{fig:RMFIT.2good}.

\begin{figure}
\begin{center}
\resizebox{6.3in}{!}{\putfig{RMFIT.e}}
\centerline{\hphantom{MM}}
\resizebox{6.3in}{!}{\putfig{RMFIT.f}}
\caption{\Hi{Top: Model field first pixel Q and U spectra, fit with two
  components.  Initial guess not as close to fit, but fit now matches
  the data.  Bottom: Spectrum re-plotted as amplitude and phase.}}
\label{fig:RMFIT.2good}
\end{center}
\end{figure}

The menu that appears at this point offers the option to {\tt QUIT}
which causes the task to quit at this point.  You may re-start later.
The option {\tt GOOD} will save this solution and go on to the next
position.  The option {\tt BAD} will mark this position as failed and go
on to the next position.   The option {\tt RE-GUESS} will loop back to
prompt you for a new guess and repeat the fit.  Options {\tt 1}, {\tt
  2}, $\ldots$, {\tt NGAUSS} will loop back to request an initial guess
with the selected number of Gaussians.  Option {\tt HAND} will prompt
you to enter using the terminal the parameters for each component.
Enter on one line for each component, the polarization brightness (in
image units), the polarization angle (in degrees), the rotation measure
and, if you are fitting a thickness, the thickness width (in radians per
meter squared).  {\tt RMFIT} will then repeat the display in
Figure~\ref{fig:RMFIT.2good} to see if you made a good guess.
Immediately after a {\tt HAND} operation only, the option {\tt DO FIT}
is offered to go back with the hand-entered values as the initial guess
for a new fit.

\begin{center}
\begin{tabular}{|l|l|}\hline
   {\tt GOOD}     & {\tt RMFIT} saves the result and goes on \\
   {\tt DO FIT}   & {\tt RMFIT} re-does the fit\\
   {\tt RE-GUESS} & {\tt RMFIT} requests a new initial guess \\
   {\tt TVOFF}    & {\tt RMFIT} saves the result and turns off
                    interactivity \\
   {\tt HAND}     & {\tt RMFIT} requests a new guess with the terminal
                    \\
   {\tt BAD}      & {\tt RMFIT} marks result bad and goes on \\
   {\tt SHOW Q\&U}& {\tt RMFIT} switches plot to Q and U \\
   {\tt SHOW P\&A} & {\tt RMFIT} switches plot to amplitude and phase\\
   {\tt 1}        & {\tt RMFIT} requests a new initial guess with one
                    component \\
   {\tt 2}        & {\tt RMFIT} requests a new initial guess with two
                    components \\
   {\tt QUIT}     & {\tt RMFIT} exits immediately \\ \hline
\end{tabular}
\end{center}

The user may select one of the {\tt SHOW} options to switch between
displays of Q and U and displays of total polarization and angle as in
Figure~\ref{fig:RMFIT.2good}.  Option {\tt TVOFF} allows you to turn off
interactivity, allowing the task to run using its own initial guesses
until it finds a completely unreasonable solution or one with an rms
greater than {\tt RMSLIMIT}\@.  When that happens, you are shown the
offending fit parameters and the task resumes with the plot at the
bottom of Figure~\ref{fig:RMFIT.redo2}, already prompting you to enter a
new guess, to allow you to try to fix things.
\hblack

\subsection{Editing and output}

Eventually all positions selected by {\tt PCUT}, {\tt ICUT}, {\tt
  BLC}, and {\tt TRC} will have been fit.  At this point, the task
computes images of the fit parameters plus the Q and U polarization at
zero wavelength and their uncertainties.  It then offers a lengthy
menu of options which will allow you to view these images and revisit
positions that seem to have produced incorrect fits.  If {\tt NGAUSS}
$> 1$, options to swap portions of image $n$ with corresponding
portions of image $m$ are also offered.  This ``edit'' menu is
illustrated in Figure~\ref{fig:RMFIT.edit2} for the {\tt NGAUSS = 2}
and {\tt DOSPIX = 3} case.  {\tt RMFIT} will replicate pixels in both
directions if needed to make small images large enough to see, but it
can also handle images larger than the TV display area.  The {\tt SET
  SPINX RANGE} appears only when {\tt DOSPIX = 1} and the {\tt SET
  THICKNESS MAX} appears only when {\tt DOSPIX =} 2, 3, or 4, the
first column of the menu shows the editing functions.

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{RMFIT.i}}
\caption{\Hi{Model field showing rotation measure for component one.
    There are pixels needing correction.
    It has been colored by {\tt OFMCOLOR}, the menu is labeled, and
    the acceptable range for component two has been set.}}
\label{fig:RMFIT.edit2}
\end{center}
\end{figure}

There are three kinds of editing implemented here.  In the first, the
user establishes the parameter extrema which should be viewed as
acceptable.  The extrema currently set are shown in the title lines.
Then, {\tt RMFIT} may be told to flag all solutions not meeting these
criteria, or, more profitably, to revisit those positions to see if a
better fit can be obtained.  The other editing methods are similar,
but act on a list of pixel positions.  These may be entered by typing
in values, by clicking on suspect pixel positions in the {\tt CURVALUE}
function described below \Hi{or by using the {\tt BLOTCH LIST} function
  also described below}.  The contents of the list may be viewed, the
solutions at the positions may be flagged, or they may be revisited to
attempt for a better solution, or the solutions at the listed positions
may be swapped between components $n$ and $m$.  The menu will offer only
appropriate swaps between components, thus 1 and 2 in our {\tt NGAUSS=2}
example.  Swapping may be required if {\tt RMFIT} gets confused as to
which component you want to call number 1 and which number 2, as
illustrated in Figure~\ref{fig:RMFIT.edit2}.  After the flagging,
revisiting, or swapping, the list is cleared.

\begin{center}
\begin{tabular}{|l|l|}\hline
 {\tt EXIT           } & Exit {\tt RMFIT}, writing output images if
                         {\tt DOOUTPUT} is now $> 0.$\\
 {\tt SET MIN S/N    } & Set minimum amplitude S/N(s) for okay
                         solutions\\
 {\tt SET MIN P1     } & Set minimum polarization at 1 GHz for okay
                         solutions\\
 {\tt SET MAX RESID  } & Set maximum residual for okay solutions\\
 {\tt SET RM RANGE   } & Set rotation measure value range(s) for okay
                         solutions\\
 {\tt SET MAX THETA ER} & Set maximum error(s) in polarization angle
                         for okay solutions\\
 {\tt SET SPINX RANGE} & Set spectral index range(s) for okay
                         solutions\\
 {\tt SET THICKNESS MAX} & Set thickness maximum value(s) for okay
                         solutions\\
 {\tt REDO ALL       } & Re-do all solutions which are not okay\\
 {\tt FLAG ALL       } & Mark bad all solutions which are not okay\\
 {\tt OFF ZOOM       } & Turn off TV zoom\\
 {\tt OFF TRANSFER   } & Turn off black \&\ white and color TV
                         enhancements\\
 {\tt RESET WINDOW   } & Display full view of current image \\
 {\tt LABEL WEDGE?   } & Turn on/off labeling of step wedge \\
 {\tt SET DOOUTPUT   } & Increment {\tt DOOUTPUT} in loop 0-3\\
\Hi{{\tt PIXRANGE}} & \Hi{To enter new plot ranges with the terminal}
                          \\
 {\tt ADD TO LIST    } & Type in output pixel coordinates to add to
                         list\\
                       & \Hi{for Q and U and for total polarization and
                           angle}\\
 {\tt SHOW LIST      } & Display pixel coordinates in list\\
 {\tt REDO LIST      } & Re-do solutions for all pixels in list\\
 {\tt FLAG LIST      } & Flag solutions for all pixels in list\\
 {\tt SWAP LIST 1-2  } & Swap solutions for components 1 and 2 for all
                         pixels in list\\ \hline
\end{tabular}
\end{center}

The second (and potentially third) menu columns contain {\tt NGAUSS}
sets of functions

\begin{center}
\begin{tabular}{|l|l|}\hline
 {\tt SHOW IMAGE P1\_1} & Enter image interaction with total
                     polarization at 1 GHz of component 1\\
 {\tt SHOW IMAGE TH\_1} & Enter image interaction with polarization
                     position angle of component 1\\
 {\tt SHOW IMAGE RM\_1} & Enter image interaction with rotation
                     measure of component 1\\
 {\tt SHOW IMAGE BE\_1} & Enter image interaction with RM thickness
                     ($\beta$) of component 1\\
 {\tt SHOW IMAGE SP\_1} & Enter image interaction with spectral index
                     of component 1\\
 {\tt SHOW IMAGE Q0\_1} & Enter image interaction with Q at 0
                     wavelength of component 1\\
 {\tt SHOW IMAGE U0\_1} & Enter image interaction with U at 0
                     wavelength of component 1\\
 {\tt SHOW IMAGE EP1\_1} & Enter image interaction with uncertainty in
                     total polarization at 1 GHz\\
%                         & of component 1\\
 {\tt SHOW IMAGE ETH\_1} & Enter image interaction with uncertainty in
                     polarization position angle\\
%                         & of component 1\\
 {\tt SHOW IMAGE ERM\_1} & Enter image interaction with uncertainty in
                     rotation measure of component 1\\
 {\tt SHOW IMAGE UBE\_1} & Enter image interaction with uncertainty in
                     RM thickness ($\beta$) of component 1\\
 {\tt SHOW IMAGE ESP\_1} & Enter image interaction with uncertainty in
                     spectral index of component 1\\
 {\tt SHOW IMAGE EQ0\_1} & Enter image interaction with uncertainty in
                     Q at 0 wavelength of component 1\\
 {\tt SHOW IMAGE EU0\_1} & Enter image interaction with uncertainty in
                     U at 0 wavelength of component 1\\ \hline
\end{tabular}
\end{center}
The {\tt SP} options will appear only if {\tt DOSPIX = 1} was specified
and the {\tt BE} options will appear only if {\tt DOSPIX =} 2, 3, or 4
was specified.

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{RMFIT.h}}
\caption{\Hi{Image screen showing model field rotation measure of
    component 1.  Pixels in the central ``source'' which were solved
    incorrectly.}}
\label{fig:RMFIT.badpix}
\end{center}
\end{figure}

When you select one of these functions most of the following operations
will appear in yet another menu.  This menu is illustrated in
Figure~\ref{fig:RMFIT.badpix}.  Only one of the {\tt LOAD AS} options
will appear, with the next one in the sequence offered when the current
one has been invoked.  The {\tt SWAP $n$-$m$} options will appear as
needed when {\tt NGAUSS} $ > 1$.  The {\tt NEXT WINDOW} option appears
when needed to display an image too large to fit on the display screen.
These options mostly invoke familiar functions from \AIPS\ to control
the {\tt FUNCTYPE} used in loading the image to the display, to enhance
the image intensities, to color the enhanced image intensities, and to
zoom the display.

\begin{center}
\begin{tabular}{|l|l|}\hline
 {\tt RETURN     } & Return to the above menus, image stays displayed\\
 {\tt LOAD AS SQ } & Re-load image with square root transfer function\\
 {\tt LOAD AS LG } & Re-load image with log transfer function\\
 {\tt LOAD AS L2 } & Re-load image with extreme log transfer function\\
 {\tt LOAD AS LN } & Re-load image with linear transfer function\\
 {\tt SET WINDOW } & Set a sub-image to view\\
 {\tt RESET WINDOW} & Return to viewing the full image\\
 {\tt OFF TRANSF } & Turn off enhancement done with {\tt TVTRANSF}\\
 {\tt OFF COLOR  } & Turn off color enhancements\\
 {\tt TVTRANSF   } & Black \&\ white image enhancement\\
 {\tt TVPSEUDO   } & Color enhancement of numerous sorts\\
 {\tt TVPHLAME   } & Color enhancement of flame type, multiple colors\\
 \Hi{{\tt OFMCOLOR}} & \Hi{Many standard color tables (OFMs)
                     available} \\
 {\tt TVZOOM     } & Interactive zooming and centering of image\\
 {\tt CURVALUE   } & Display value under cursor, mark pixels for list\\
\Hi{{\tt BLOTCH LIST}} & \Hi{Interactively add a blotch area to the
                      pixel list} \\
 {\tt SWAP 1-2   } & Swap solutions for components 1 and 2
                     interactively\\
 {\tt NEXT WINDOW} & Move to next window in large images\\ \hline
\end{tabular}
\end{center}

Two operations in this menu are different from the usual.  {\tt
  CURVALUE} provides the capability of selecting positions for the
edit ``list.''  During the {\tt CURVALUE} operation position the
cursor over the desired pixel and press buttons {\tt A}, {\tt B}, or
{\tt C} to add that pixel to the list.  The \Hi{{\tt BLOTCH LIST} and}
{\tt SWAP $n$-$m$} operations use a TV blotch operation like that in the
\AIPS\ verb {\tt TVSTAT} and task {\tt BLSUM}\@.  You are to mark with a
``blotch'' region those pixels in the present image which are to \Hi{be
  added to the pixel list or} to have their solutions swapped with those
of the selected component.  Instructions will appear in the message
window as you proceed.  Begin by positioning the cursor at a pixel to be
the first vertex of a connected sequence of vertices and press TV button
{\tt A}\@.  Move to the next vertex and press button {\tt A} again and
repeat until you have marked all vertices for this region.  Then press
button {\tt D} if you are done with this region or button {\tt C} if you
need to re-position one of the vertices.  In this case, move the cursor
to the vertex to be re-positioned, press button {\tt A} and drag the
vertex to the corrected position.  Press button {\tt A} or {\tt B} to
fix that vertex and go on to reset another vertex or {\tt D} to end this
region and \Hi{list or} swap the solutions.  You may do this as many
times as needed.

When you have finished getting the images just the way you want them,
you may write them out as \AIPS\ image files.  Select the {\tt SET
  DOOUTPUT} option until its value, shown at the top of the screen, is
what you want.  In {\tt RMFIT}, values 1 and 3 cause residual image
cubes to be written, while values 2 and 3 cause images of the
parameter values and their uncertainties to be written.  The Q and U
residual images get classes {\tt QRESID} and {\tt URESID},
respectively.  The polarization at 1 GHz, polarization angle, rotation
measure, and spectral index or thickness (if any) images and their
uncertainties get classes {\tt PPOL$n$}, {\tt THETA$n$}, {\tt
  ROTME$n$}, {\tt SPIX$n$}, {\tt THICK$n$}, {\tt DPPOL$n$}, {\tt
  DTHET$n$}, {\tt DROTM$n$}, {\tt DSPIX$n$}, and {\tt DTHIC$n$},
respectively, and Q and U at zero wavelength get classes {\tt
  Q0\_$n$}, {\tt U0\_$n$}, {\tt DQ0\_$n$}, and  {\tt DU0\_$n$}, where
$n$ is the component number.
%\vfill\eject

\section{Post-fit plotting}

The images produced by {\tt XGAUS}, {\tt ZEMAN} and {\tt RMFIT} may be
displayed using all the usual tools such as {\tt KNTR}, {\tt PROFL},
and numerous other tasks.  However, the display of the spectral data
and the various fits to them required new tasks.

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{XG2PL}}
\caption{{\tt XG2PL} example spectrum from a circular area for a
  Gaussian and Zeeman fit.}
\label{fig:XG2PL}
\end{center}
\end{figure}

{\tt XG2PL} plots a spectrum for a single pixel or for a rectangular
or circular region about a single pixel.  For each pixel included in
the average, the task reads the I polarization image to obtain the
data and the appropriate line(s) in the {\tt XG} table to obtain the
{\tt XGAUS} solution for that pixel.  It then computes the spectrum of
each component in the model, plus the sum of the components, and the
residual (data-model).  Each of these then enter into the average of
that parameter.  Finally, the task plots a user-selected number of the
parameters and, optionally, prints all of them to a text file.  The
plot may appear on the TV or be placed in a standard plot file
attached to the I polarization image.

Optionally, {\tt XG2PL} will also add the spectrum of the V
polarization data and the results of the fitting done by {\tt ZEMAN}
at the same pixels as the I polarization Gaussians.  The task reads
the V polarization image for the data spectrum and the {\tt ZE} table
for the Zeeman-splitting model (either using the Gaussians or the
simpler ones using the I polarization data).  It computes the spectrum
of each component (including the gain term in each), the net model
(sum of the components but including the gain term only once), and the
residual.  Finally, the task plots a user-selected number of the
parameters and, optionally, prints all of them to a text file.  In
general, the I polarization spectrum appears in the upper part of the
plot and the V polarization spectrum appears in the lower portion.
Either portion may be omitted under control of the adverbs.  The
output of {\tt XG2PL} is illustrated in Figure~\ref{fig:XG2PL}.

The inputs for {\tt XG2PL} begin with the I polarization image which
is required and then the V polarization image which may be omitted.
{\tt INVERS} and {\tt IN2VERS} give the version numbers of the {\tt
  XG} and {\tt ZE} tables, respectively.  {\tt APARM} provides the
central pixel coordinates, the plot intensity ranges for I and for V,
the size of the rectangle (or circle) over which to average, and a
flag limiting the average to those pixels having a model fit,  {\tt
  BPARM} is a set of flags selecting which parameters are plotted.
{\tt CPARM} selects the channel range to be plotted, the relative size
of the V and I plots, the type of the horizontal axis (channels,
frequency, velocity), whether the channels are plotted in reverse
order, and whether the data are plotted with stepped or directly
connected lines.  {\tt OUTTEXT} specifies the output text file, if
any.  The usual {\tt XYRATIO}, {\tt LTYPE}, {\tt DOTV}, and {\tt
  GRCHAN} adverbs control the scale, labeling, choice of TV versus
plot file, and, if TV, which graphic channel(s) are used.  If {\tt
  GRCHAN} is zero, graphics channel 1 is used for data and labeling
(usually yellow), channel 2 is used for the full model (usually
green), channel 3 is used for the residual (usually pink), and channel 4
is used for the model component(s) (usually cyan).
\vfill\eject

\begin{figure}
\begin{center}
\resizebox{6.0in}{!}{\putfig{RM2PL}}
\caption{{\tt RM2PL} example spectrum from a circular area for a
  2-component rotation measure fit.}
\label{fig:RM2PL}
\end{center}
\end{figure}

{\tt RM2PL} plots a spectrum for a single pixel or for a rectangular
or circular region about a single pixel.  For each pixel included in
the average, the task reads the Q and U polarization images to obtain
the data and the appropriate line(s) of the {\tt RM} table to obtain the
{\tt RMFIT} solution for that pixel.  It then computes the spectrum of
each component in the model, plus the sum of the components, and the
residual (data-model).  Each of these then enter into the average of
that parameter.  Finally, the task plots a user-selected number of the
parameters and, optionally, prints all of them to a text file.  The
plot may appear on the TV or be placed in a standard plot file
attached to the Q polarization image.  The user may choose to omit all
of the Q or all of the U polarization plot.  The output of {\tt RM2PL}
is illustrated in Figure~\ref{fig:RM2PL}.

The inputs for {\tt RM2PL} begin with the Q polarization image which
is required and then the U polarization image which may be omitted.
{\tt INVERS} gives the version number of the {\tt RM} table.  {\tt
  APARM} provides the central pixel coordinates, the plot intensity
ranges for Q and for U, the size of the rectangle (or circle) over
which to average, and a flag limiting the average to those pixels
having a model fit.  {\tt BPARM} is a set of flags selecting which
parameters are plotted.  {\tt CPARM} selects the channel range to be
plotted, the relative size of the U and Q plots, the type of the
horizontal axis (channels, frequency, wavelength, wavelength squared),
whether the channels are plotted in reverse order, and whether the
data are plotted with stepped or directly connected lines.  {\tt
  OUTTEXT} specifies the output text file, if any.  The usual {\tt
  XYRATIO}, {\tt LTYPE}, {\tt DOTV}, and {\tt   GRCHAN} adverbs
control the scale, labeling, choice of TV versus plot file, and, if
TV, which graphic channel(s) are used.  If {\tt GRCHAN} is zero,
graphics channel 1 is used for data and labeling (usually yellow),
channel 2 is used for the full model (usually green), channel 3 is
used for the residual (usually pink), and channel 4 is used for the
model component(s) (usually cyan).
\vfill\eject

\hblack
\section{Model creation}

\subsection{Visibility data models}

There are two tasks which will modify an existing visibility (``$uv$'')
data set, adding models and noise to the existing data or replacing
those data.  {\tt UVMOD} can apply a full set of calibration and data
selection adverbs to the input $uv$ data set.  Then the data are
scaled by {\tt FACTOR}, which can be zero to eliminate the input data
entirely.  Finally, a model is added.  The model can be a simple set
of up to 4 {\tt NGAUSS} components specified by adverbs {\tt CTYPE},
{\tt FMAX}, {\tt FPOS}, and {\tt FWIDTH}\@.  Alternatively, adverb
{\tt INLIST} may specify a text file containing up to 9999 components.
Each non-comment line in this file specifies
\vfill\eject
\begin{center}
\begin{tabular}{|r|l|}\hline
 1. & I-polarization flux (Jy at header frequency)\\
 2. & X-shift from reference (arc seconds)\\
 3. & Y-shift from reference (arc seconds)\\
 4. & Major axis (arc seconds)\\
 5. & Minor axis (arc seconds)\\
 6. & Position angle (degrees CCW from North)\\
 7. & Type code: 1 Gaussian, 2 solid disk, 3 solid rectangle,\\
    & 4 optically thin sphere, 5 exponential, else point\\
 8. & Spectral index (see note below)\\
 9. & Spectral index curvature\\
10. & Q-polarization flux (Jy at header frequency)\\
11. & U-polarization flux (Jy at header frequency)\\
12. & V-polarization flux (Jy at header frequency) \\ \hline
\end{tabular}
\end{center}
where trailing zeros may be omitted.  If spectral index is used, note
the definition of spectral index and flux.  If $x$ is the logarithm,
base 10, of the frequency divided by the header frequency, then
$$ F = F_0 10^{(\alpha x + \beta x^2)} $$
where $\alpha$ is the spectral index and $\beta$ is its curvature.
Thus the flux you must specify in the model is the flux at the
frequency of the data.

{\tt SPMOD} is a very similar task designed to do spectral-line
modeling rather than spectral index and polarization.  In this case
{\tt INLIST} is required and each of the up to 9999 non-comment lines
contains
\begin{center}
\begin{tabular}{|r|l|}\hline
 1. & RR-polarization flux (Jy)\\
 2. & X-shift from reference (arc seconds)\\
 3. & Y-shift from reference (arc seconds)\\
 4. & Major axis (arc seconds)\\
 5. & Minor axis (arc seconds)\\
 6. & Position angle (degrees CCW from North)\\
 7. & Type code: 1 Gaussian, 2 solid disk, 3 solid rectangle,\\
    & 4 optically thin sphere, 5 exponential, else point\\
 8. & Spectral line center (channels, but see note below)\\
 9. & Spectral line width (FWHM channels)\\
10. & LL-polarization flux (Jy)\\ \hline
\end{tabular}
\end{center}
where trailing zeros may be omitted.  The LL flux will be set equal to
the RR flux if it is omitted.  But note that both the RR and LL flux
may be specified to have any value positive, zero, or negative.  This
will allow Zeeman-splitting to be modeled with any degree of
splitting.  Spectral channels are counted from {\tt BCHAN} through
{\tt ECHAN} as $1$ to $N$ in the first spectral window ({\tt BIF}),
then $N+1$ through $2 N$ in the next spectral window, and so forth.
This allows you to put lines in any spectral window you desire.

\subsection{Image data models}

There are three tasks designed to modify an image adding a specified
model and noise.  {\tt IMMOD} is intended for continuum images.  The
existing image may be scaled (or eliminated) and noise added.  It adds
{\tt NGAUSS} components up to four using adverbs {\tt OPCODE} to
specify type, and {\tt FMAX}, {\tt FPOS}, and {\tt FWIDTH} to specify
peak brightness and pixel position and size.  Alternatively, adverb
{\tt INLIST} may be used to specify up to 9999 components, one per
line in the text file. Each non-comment line specifies
\vfill\eject

\begin{center}
\begin{tabular}{|r|l|}\hline
 1. & Peak brightness (Jy/beam)\\
 2. & Component $X$ center (pixels)\\
 3. & Component $Y$ center (pixels)\\
 4. & Major axis (pixels)\\
 5. & Minor axis (pixels)\\
 6. & Position angle (degrees CCW from North)\\ \hline
\end{tabular}
\end{center}

The more recent task {\tt MODIM} is specifically designed for
polarization modeling and has been used extensively in testing {\tt
  RMFIT}\@.  It takes in three image cubes, one for I polarization,
one for Q, and one for U.  Alternatively, it can create new images
with no input images used.  Up to 9999 components are taken from the
text file specified by {\tt INLIST}\@.  Each non-comment line in this
file specifies {\it every one of}
%\vfill\eject
\begin{center}
\begin{tabular}{|r|l|}\hline
 1. & I-polarization brightness (Jy/beam at 1 GHz)\\
 2. & Q-polarization brightness (Jy/beam at 1 GHz)\\
 3. & U-polarization brightness (Jy/beam at 1 GHz)\\
 4. & Spectral index (see note below)\\
 5. & Rotation measure (radians $m^{-2}$)\\
 6. & Rotation measure thickness (radians $m^{-2}$)\\
 7. & Component $X$ center (pixels)\\
 8. & Component $Y$ center (pixels)\\
 9. & Major axis (pixels)\\
10. & Minor axis (pixels)\\
11. & Position angle (degrees CCW from North)\\
12. & Type code: 1 Point, 2 Gaussian, 3 solid disk, 4 solid rectangle,\\
    & 5 optically thin sphere, 6 exponential, else point\\ \hline
\end{tabular}
\end{center}
Note that all parameters must be specified for every component.
Adverb {\tt OPTYPE} specifies the type of rotation measure thickness
model, with values {\tt 'SLAB'}, {\tt 'GAUS'}, and {\tt 'EXP '}\@.
The brightness of a component at frequency $\nu$ in GHz is $F = F_1
\nu^\alpha$\@.  Input (if any) and output images are in the
transposition that places celestial coordinates first and frequency on
the third axis.

A third modeling task, {\tt MODSP}, models I and V image cubes from
models similar to those of {\tt MODIM}, but suitable for spectral lines.
It takes in two image cubes, one for I or RR polarization and one for V
or LL polarization.  Alternatively, it can create new images with no
input images used.  Up to 9999 components are taken from the text file
specified by {\tt INLIST}\@.  Each non-comment line in this file
specifies

\begin{center}
\begin{tabular}{|r|l|}\hline
 1. & R-polarization brightness at line center (Jy/beam)\\
 2. & L-polarization brightness at line center (Jy/beam)\\
 3. & Component $X$ center (pixels)\\
 4. & Component $Y$ center (pixels)\\
 5. & Major axis (pixels)\\
 6. & Minor axis (pixels)\\
 7. & Position angle (degrees CCW from North)\\
 8. & Type code: 1 Point, 2 Gaussian, 3 solid disk, 4 solid rectangle,\\
    & 5 optically thin sphere, 6 exponential, else Gaussian\\
 9. & Line center (channels)\\
10. & Line width (FWHM of Gaussian in channels)\\
11. & Change of line center per $X$ pixel (channels/pixel)\\
12. & Change of line center per $Y$ pixel (channels/pixel)\\
13. & Change of line width per $X$ pixel (channels/pixel)\\
14. & Change of line width per $Y$ pixel (channels/pixel)\\ \hline
\end{tabular}
\end{center}

Note that parameters 1--10 must be specified for every component.  The
remaining position and width derivatives will be taken as zero if
omitted.  Input (if any) and output images are in the transposition
that places celestial coordinates first and frequency on the third
axis.  Input images may be in either I and V polarizations or RR and
LL, output images will be in I and V.  Input images are scaled by {\tt
  FACTOR} before the model is added.  {\tt MODSP} allows the spectral
channel parameters to be linear functions of position around the
center of the component.  Thus, the line center channel and width are
given by
\begin{eqnarray*}
  C(x,y) & = & C_i + \frac{dC_i}{dx} (x - X_i) +
                     \frac{dC_i}{dy} (y - Y_i) \\
  W(x,y) & = & W_i + \frac{dW_i}{dx} (x - X_i) +
                     \frac{dW_i}{dy} (y - Y_i) \\
\end{eqnarray*}
where, for component $i$, $C_i$ is parameter 9 above,
$\frac{dC_i}{dx}$ is parameter 11, $\frac{dC_i}{dy}$ is parameter 12,
$X_i$ is parameter 2, $Y_i$ is parameter 3. $W_i$ is parameter 10,
$\frac{dW_i}{dx}$ is parameter 13, and $\frac{dW_i}{dy}$ is parameter
14.

\end{document}
