; ZEMAN
;---------------------------------------------------------------
;! Fits 1-dimensional Zeeman model to data
;# TASK ANALYSIS PLOT SPECTRAL ONED MODELING TV-APPL POLARIZATION
;-----------------------------------------------------------------------
;;  Copyright (C) 2013
;;  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
;-----------------------------------------------------------------------
ZEMAN     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
ZEMAN     Fits 1-dimensional Zeeman model to images of I and V
INNAME                             Input V image name (name)
INCLASS                            Input V image name (class)
INSEQ             0.0     9999.0   Input V image name (seq. #)
INDISK            0.0        9.0   Input V image disk unit #
IN2NAME                            Input I image name (name)
IN2CLASS                           Input I image name (class)
IN2SEQ            0.0     9999.0   Input I image name (seq. #)
IN2DISK           0.0        9.0   Input I image disk unit #
OUTNAME                            Output image name (name)
OUTCLASS                           Output image name (class)
OUTSEQ           -1.0     9999.0   Output image name (seq. #)
OUTDISK           0.0        9.0   Output image disk unit #.
BLC                                Bottom left corner of input
TRC                                Top right corner of input
FLUX              0.0              Flux cutoff (see HELP)
INVERS            0.0              ZE table versions
OPTYPE                             'GAUS', '2SID', other
IN2VERS           0.0              Mode GAUS: input XG table
DOOUTPUT         -1.0        7.0   Write corrected image: bit 1
                                   write parameter img: bit 2
                                   Write residual image: bit 3
DOTV             -1.0        1.0   Plot data on TV
DOMAX           -1.0         1.0   Solve for constant scaling I
BADDISK                            Disk to avoid for scratch
----------------------------------------------------------------
ZEMAN
Task:  Fits one-dimensional Zeeman model to spectral cubes.  The
       simplest model is
           V(c) = A * I(c) + 0.5 * B * dI(c)/dc
       where c is spectral channel or frequency and A and B are
       determined for each celestial coordinate by a least-squares
       fitting routine.  If the I image has been through XGAUS, then a
       more complicated model may be fit:
           V(c) = A * I(c) + 0.5 * sum [ B_i * dI_i/dc ]
       where the i loops from 1 to the number of Gaussians fit at
       each celestial pixel.  I is the actual observed unpolarized
       brightness while I_i is the Gaussian brightness of the i'th
       component.

       ZEMAN can be run in a batch mode in which it simply finds the
       solution to the above models and writes out images of the
       parameters, of a corrected V cube, and of the residual.  If
       DOTV is true, it will use an interactive mode to display each
       fit and ask whether the solution should be kept or marked bad.
       When you get tired of the interaction, you can turn off the TV
       and fit the rest of the cube in the batch mode.  Interactivity
       will then resume, displaying images of the parameters and
       allowing you to flag more solutions interactively.

       Like XGAUS, ZEMAN works with a table of solutions called a ZE
       table for ZEMAN.  It begins by filling in the table with null
       solutions and then calls a routine to determine solutions pixel
       by pixel for those pixels which are strong enough (FLUX).  This
       can be run interactively or in batch and one can turn off
       interactivity after establishing that all is working well.
       After all pixels are fit, the task will go into an image
       display and editing phase using the TV, but only if the task is
       started with DOTV true.  You can refit or flag selected
       solutions at this stage.  Finally, selected output images are
       written from the fit.  Note, during interaction you may simply
       quit for now and resulme at some later time using the table
       which has already been partially (or fully) filled in.

Adverbs:
  INNAME.....Transposed V image name (name).       Standard defaults.
  INCLASS....Transposed V image name (class).      Standard defaults.
  INSEQ......Transposed V image name (seq. #).     0 => highest.
  INDISK.....Transposed V image disk #.   0 => any.
  IN2NAME....Transposed I image name (name).       Standard defaults.
  IN2CLASS...Transposed I image name (class).      Standard defaults.
  IN2SEQ.....Transposed I image name (seq. #).     0 => highest.
  IN2DISK....Transposed I image disk #.   0 => any.
  OUTNAME....Output image name (name).      Standard defaults.
  OUTCLASS...Output image name (class).     Standard defaults.
             Used only for the corrected V cube.  The OUTCLASSes for
             the fit parameters are GAIN and FIELDn for A and B_n.
  OUTSEQ.....Output image name (seq. #).   0 => highest unique.
  OUTDISK....Disk drive # of output image. 0 => highest number with
             sufficient space.
  BLC........Bottom right corner in input image of desired subimage.
             Default is entire image.  NOTE: in GAUS mode, the BLC and
             TRC of XGAUS are used instead.
  TRC........Top right corner in input image of desired subimage.
             Default is entire image.
  FLUX.......A flux cutoff in the same units as the input I image
             (i.e. Jy/beam).  If a row does not have three consecutive
             points above this level, no fit is done to the row.
  INVERS.....ZE table version: 0 => make a new one.
  OPTYPE.....'GAUS' use a pre-computed XG table (from XGAUSS) to do
                    the second model above.
             '2SID' do the first model using a "2-sided" derivative:
                    dI(c)/dc = 0.5 * (I(c+1) - I(c-1))
             else:  do the first model using a "1-sided" derivative
                    dI(c)/dc = I(c) - I(c-1)
  IN2VERS....For OPTYPe='GAUS' only: Input XG table version 0 -> high
  DOOUTPUT...= 1,3,5, or 7 requests that the V image be corrected by
                      the solution for A and written out using class
                      OUTCLASS
             = 2,3,6, or 7 requests that the model images (A, B_i) be
                      written out as cataloged images with classes
                      GAIN and FIELDi.
             = 4,5,6, or 7 requests that the residual (data-model)
                      image be written out as a cataloged image with
                      class 'VRESID'
             Note that DOOUTPUT can be changed interactively during
             the imaging portion of the task (if DOTV > 0) and the
             value of the adverb at exit is all that matters.
  DOTV.......False (<= 0) causes the task to run in a batch mode.
             True (> 0) implies plot the data on the TV
  DOMAX......<= 0. means hold this parm fixed (NGAUSS > 1 only).
  BADDISK....Disk drives to avoid for scratch files.
----------------------------------------------------------------
ZEMAN:  One-dimensional Zeeman fitting of I and V cubes
Documenter:  E. W. Greisen NRAO
Related Programs: CUBIT, GAL, IMFIT, JMFIT, SAD, SLFIT

Beginning with the 31DEC13 version of AIPS, ZEMAN has been overhauled
into the task it should have been.  The interactive Tektronix option
was removed since Tek emulation in xterm windows has proved unreliable
and since fill TV menus were added to the later stages of the task.
At present, ZEMAN assumes that there are at most 3 axes with more than
one pixel in the input image.

ZEMAN begins by making an XG table containing one row for every pixel
in the output plane from BLC(2,3) to TRC(2,3).  Each row of the table
is initialized with the peak value in each row of the input image,
where the "peak" is the largest value for 3 adjacent pixels.  Then
ZEMAN begins fitting those rows having a peak value above FLUX, doing
every YINC pixel in the Y axis and ZINC pixel in the Z axis.  Initial
guesses for each fit are found from the row data when NGAUSS = 1 or
from a previous fit for NGAUSS > 1.  This is one of the reasons the
interactive mode is recommended strongly for NGAUSS > 1.  After each
input image row is fit, the results are added to the XG table
immediately.  This allows the interactive user to quit at any time and
then restart the process later (set INVERS to point at the
pre-existing XG table which you intend to modify).  After the YINC by
ZINC pass through the data, a second pass through the data is done on
every pixel.  Those pixels which have already been fit are skipped,
but the results from them contribute to the next initial guess.

Finally, when all pixels have been fit, ZEMAN enters an interactive
routine designed to improve the results before any output images are
written.  However, before describing that function, we need to provide
details of the interactive fitting of each row of the input image.
This process is followed whenever a row is fit in the first YINC/ZINC
function, in the second every pixel function, and in the third result
editing function.

For each row, the first step in the fitting process is to determine an
initial guess non-interactively.  If ZEMAN is currently in interactive
mode, the next steps are:
    1. If DOSLICE > 0 (STRONGLY recommended), the input data are
       plotted on the TV.  If DOSLICE > 1 (also recommended), the
       initial guess is added to the plot of the data.
    2. If DOSLICE > 1 and either NGAUSS > 1 or a fitting has been
       forced, ZEMAN asks in the AIPS window for instructions.  If
       your answer begins with E or e, ZEMAN will ask you to enter a
       new initial guess (step 3).  If the answer begins with B or b,
       ZEMAN will mark this pixel as bad and, if the answer begins
       with Q or q, ZEMAN will simply exit, allowing you to restart it
       at a later time.  If the answer is anything else, e.g. a simple
       carriage return ("Enter"), ZEMAN goes on to step 4 below.
    3. To enter a new initial guess, you will be prompted to position
       the TV cursor at the center and height of each Gaussian (note
       that the X and Y positions matter) and at the half-width of each
       Gaussian (note that only the X position matters).  To mark each
       point, hit any button (A, B, C, or D).  ZEMAN then returns to
       step 1 above to plot the new guess and ask again.
    4. Once an acceptable initial guess has been found, ZEMAN proceeds
       to call a non-linear least squares fitting routine to determine
       the Gaussian parameters that appear to fit the data best.  the
       answers are then checked to see if they are "reasonable" -
       negative components, components centered outside the input
       data, and components < 1.5 cells wide (FWHM) are thought to be
       unreasonable.
    5. If the result is unreasonable and the current mode is not
       interactive, but it started as interactive, the interactive
       mode is turned back on (with a message and a display of the
       answers in the message terminal window) and ZEMAN returns to
       step 1.  If the task did not start as interactive and the
       result is unreasonable, the solution for that input row is
       marked bad and the task goes on to the next row.
    6. If the current mode is interactive, the residuals are added to
       the plot (if DORESID > 0) and the final model is also added to
       the plot (if DOMODL > 0).  In the input terminal, you are then
       told if the answers are "unreasonable" and are shown the
       answers, reasonable or not.  You are then offered a variety of
       choices.  If your response begins with:
          a) B or b - the solution is marked as bad and ZEMAN goes on
                      to the next input row.
          b) Q or q - ZEMAN exits cleanly, leaving the XG table to be
                      worked on again at a later time.
          c) R or r - ZEMAN returns to step 1 to try again with the
                      current number of Gaussians (<= NGAUSS)
          c) E or e - ZEMAN returns to step 1 to try again with the
                      current number of Gaussians (<= NGAUSS)
          e) 1, 2, 3, or 4 - ZEMAN returns to step 1 to try again with
                      the specified number of Gaussians (note that
                      only numbers <= NGAUSS are respected)
          f) H or h - ZEMAN prompts you to enter the Gaussian
                      parameters for each component.  You enter the
                      peak in input image units, the center in pixels
                      wrt the reference pixel, and width in pixels.  All
                      3 numbers must appear in one line.  It will
                      prompt for all NGAUSS components, but will
                      change the prompt if the current number of
                      Gaussians is < NGAUSS.  If you enter something
                      besides 0 in that case, the the current number
                      of Gaussians is increased appropriately.  ZEMAN
                      then returns to the start of this step 6 to
                      allow you to see if you made a good guess.
          g) T or t - ZEMAN turns off the interactive mode and does
                      solutions in a batch-like fashion until an
                      unreasonable solution is found or it finishes
                      the current function.
          h) Anything else - ZEMAN saves the answer (with
                      uncertainties) and goes on to the next input
                      row.
       Note, on retries (R, E, 1, 2, 3, and 4), step 2 will skip the
       questions and go on to step 3 directly (as if you had answered
       E).
Recommended adverb values are DOSLICE = 2; DOMODEL = 1; GRCHAN = 0;
and probably DORESID = 1.  This allows full information and
interactivity and uses different graphics planes for the different
plots.  When GRCHAN is not zero, old plots of guesses, residuals, and
fits may still be visible when new solutions are plotted.  (Otherwise
the plot of the input data will be erased.)

After all pixels that are strong enough have been given a solution,
ZEMAN enters a menu-driven function.  The menu has in the left column:
-------------------
| EXIT            |   Exit ZEMAN, writing output images if DOOUTPUT is
                      now > 0.
| SET MIN S/N     |   Set minimum amplitude S/N(s) for okay solutions
| SET MAX RES     |   Set maximum residual for okay solutions
| SET OFFX RANGE  |   Set offset range(s) for okay solutions
| SET WIDTH RANGE |   Set width range(s) for okay solutions
| SET MAX ERR WID |   Set maximum error(s) in width for okay solutions
| REDO ALL        |   Re-do all solutions which are not okay
| FLAG ALL        |   Mark bad all solutions which are not okay
| OFF ZOOM        |   Turn of TV zoom
| OFF TRANSFER    |   Turn off black & white and color TV enhancements
| SET DOOUTPUT    |   Increment DOOUTPUT in loop 0-3 - with 1 and 3
                      causing residual images and 2 and 3 causing
                      parameter images to be written on EXIT
| ADD TO LIST     |   Type in output pixel coordinates to add to list
| SHOW LIST       |   Display coordinates in list
| REDO LIST       |   Re-do solutions for all pixels in list
| FLAG LIST       |   Flag solutions for all pixels in list
| SWAP LIST 1-2   |   Swap solutions for components 1 and 2 for all
                      pixels in list
| SWAP LIST 1-3   |   Swap solutions for components 1 and 3 for all
                      pixels in list
| SWAP LIST 2-3   |   Swap solutions for components 2 and 3 for all
                      pixels in list
| SWAP LIST 1-4   |   Swap solutions for components 1 and 4 for all
                      pixels in list
| SWAP LIST 2-4   |   Swap solutions for components 2 and 4 for all
                      pixels in list
| SWAP LIST 3-4   |   Swap solutions for components 3 and 4 for all
                      pixels in list
-------------------
Only appropriate SWAP LIST options are actually displayed.

Three editing concepts are present in this menu.  The first is to
establish what parameter values (S/N in peak, peak residual, center
values, width values, and error in widths) for each component
constitute an "okay" solution.  The selected limits are displayed
above the menu.  Then you can choose to try to REDO ALL solutions
which are not okay or simply FLAG ALL, marking them as bad.  The
second concept is to create a list of up to 1000 pixels which are
thought to have some problem.  You can enter these pixels by hand (or
remove pixels from the list by hand) with ADD TO LIST.  A faster way
to add to the list is to show the image of some parameter and then
select interesting pixels while doing CURVALUE (see below).  The list
of pixels can be re-done (REDO LIST) or marked as bad (FLAG LIST).
The list is cleared by these operations so you can make a new list as
needed.  The third editing method available here is to take the pixel
list and swap the solutions at those pixels between component N with
component M.

The right hand column of options include:
------------------
| SHOW IMAGE A1  |   Enter image interaction with peak value of
                     component 1
| SHOW IMAGE C1  |   Enter image interaction with center pixel of
                     component 1
| SHOW IMAGE W1  |   Enter image interaction with width of
                     component 1
| SHOW IMAGE F1  |   Enter image interaction with "flux" of
                     component 1
| SHOW IMAGE EA1 |   Enter image interaction with uncertainty in
                     peak value of component 1
| SHOW IMAGE EC1 |   Enter image interaction with uncertainty in
                     center pixel of component 1
| SHOW IMAGE EW1 |   Enter image interaction with uncertainty in
                     width of component 1
| SHOW IMAGE EF1 |   Enter image interaction with uncertainty in
                     "flux" of component 1
------------------
For NGAUSS > 1, appropriate additional choices are offered.

When you select one of the above options, the above menus are turned
off, the image plane maintained in memory is displayed on the TV
screen, and a new menu is displayed.  If offers:
---------------
| RETURN      |   Return to the above menus, image stays displayed
| LOAD AS SQ  |   Re-load image with square root transfer function
| LOAD AS LG  |   Re-load image with log transfer function
| LOAD AS L2  |   Re-load image with extreme log transfer function
| LOAD AS LN  |   Re-load image with linear transfer function
| OFF TRANSF  |   Turn off enhancement done with TVTRANSF
| OFF COLOR   |   Turn off color enhancements
| TVTRANSF    |   Black & white image enhancement
| TVPSEUDO    |   Color enhancement of numerous sorts
| TVPHLAME    |   Color enhancement of flame type, multiple colors
| TVZOOM      |   Interactive zooming and centering of image
| CURVALUE    |   Display value under cursor, mark pixels for list
| SWAP 1-2    |   Swap solutions for components 1 and 2 interactively
| SWAP 1-3    |   Swap solutions for components 1 and 3 interactively
| SWAP 2-3    |   Swap solutions for components 2 and 3 interactively
| SWAP 1-4    |   Swap solutions for components 1 and 4 interactively
| SWAP 2-4    |   Swap solutions for components 2 and 4 interactively
| SWAP 3-4    |   Swap solutions for components 3 and 4 interactively
| NEXT WINDOW |   Move to next window into large images
---------------
For NGAUSS < 4, appropriate SWAP options are suppressed.

Only one of the LOAD AS xx options is offered - SQ when the current
function type is LN, LG when the current function type is SQ, L2 when
the current function type is LG, and LN when the current type is L2.
Selecting this option, reloads the image with the newly selection
function type and changes the menu option accordingly.  The options
OFF TRANSF through TVZOOM are essentially the same as the
corresponding verbs in AIPS.  CURVALUE is like the verb of that name,
but, when you press buttons A or B, the pixel under the cursor is
added to the editing list.

The SWAP n-m options invoke the interactive polygon setting operation
used by the AIPS verb TVSTAT.  This starts in the "set polygon n" mode
where button A selects a vertex in the polygon, B sets the last vertex
in the polygon and gets ready to set polygon n+1, C sets the last
vertex and enters a vertex editing mode, and D sets the last vertex
and exits the polygon setting.  In the vertex editing mode, move the
cursor to a vertex to be moved and press button A or B to reset its
position.  After moving it to the desired point, press button A or B
to fix that point and restart the vertex editing mode.  Press button C
to fix that point and then start creating polygon n+1.  Press button D
to fix that point and exit the polygon setting.  Once button D has
been set, the selected areas of pixels have their solutions swapped.
The images are updated and reloaded automatically.

Very large images may not be able to fit on your TV.  When the image
exceeds the size of the TV, the top line will show the component
number followed by a subimage number in parentheses and the NEXT
WINDOW option will appear.  The first sub-images displayed is called
number 0 and shows the full image every n'th pixel in X and Y.
Sub-image 1 begins at the lower left, moves right, then back to the
left and up, and so forth until the top right is reached.  Every pixel
is displayed in these sub-images.  Use the NEXT WINDOW option to step
through the sub-images in a circular fashion.
