; BDEPO
;---------------------------------------------------------------
;! computes depolarization due to rotation measure gradients
;# Task Polarization Analysis
;-----------------------------------------------------------------------
;;  Copyright (C) 1995
;;  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
;-----------------------------------------------------------------------
;---------------------------------------------------------------
BDEPO     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
BDEPO: Task to find beam depolarization from RM gradient
USERID     -32000.0     32000.0    User ID. 0=current user
                                     32000=>all users
INNAME                             First RM grad image name
INCLASS                            First RM grad image class
INSEQ           0.0      9999.0    First RM grad image seq. no.
INDISK          0.0         9.0    First RM grad image disk no.
IN2NAME                            Second RM grad image name
IN2CLASS                           Second RM grad image class
IN2SEQ          0.0      9999.0    Second RM grad image seq. no.
IN2DISK         0.0         9.0    Second RM grad image disk no.
BLC                                BLC of input image(s)
TRC                                TRC of input image(s)
OUTNAME                            Output DPR image name
OUTCLASS                           Output DPR image class
OUTSEQ          0.0      9999.0    Output DPR image seq. #
OUTDISK         0.0         9.0    Output DPR image disk drive
CPARM                              User supplied parameters
---------------------------------------------------------------
BDEPO
Type: Task
Use:  Depolarization can be produced if the Rotation Measure
      gradient is sufficiently large across the synthesized
      beam.  BDEPO tries to predict the depolarization ratio
      (DPR) between       two frequencies owing to a linear
      Rotation Measure (RM) gradient across the beam.  The
      gradient image(s) may be calculated from your RM image
      (with NINER), or may be models of a foreground screen
      that you have generated in some other way.  The output
      DPR image can be plotted (with IMVIM) against an observed
      DPR image in order to asses the importance of beam
      depolarization.  The observed DPR should be m(F1)/m(F2)
      where F1 < F2 are the frequencies, and m is the degree of
      polarization.

      It assumed that the intrinsic degree of polarization is
      constant across the beam.  A beam of arbitrary
      orientation must be present in the header and you must
      supply the two frequencies between which to calculate the
      DPR.  If the RM structure is still possibly unresolved in
      your highest resolution RM image, you will run the risk
      of underpredicting the beam depolarization.

      This task runs in 2 modes.   First, the more general case
      is treated if you input 2 RM gradient images.   The first
      must be the x (RA) direction gradient, and the second the
      y (DEC) direction gradient.   The second mode is invoked
      if you input only one RM gradient image.  In this case it
      is assumed that the x and y components of the RM gradient
      are equal and given by the input image.

      The units of the gradient images should be reduceable to
      RAD/M/M/PIXEL by scale factors which you may supply via
      the CPARM(3) and CPARM(4) parameters. For example,
      consider BDEPO working in the second mode mentioned
      above.  A generalised gradient image can be produced by
      NINER with opcode 'SOBL'.  If DPARM(1) = 0, NINER makes a
      gradient of the form :

                G = 8 * (g_x**2 + g_y**2)**0.5

      where g_x and g_y are x and y directional derivatives.
      Now, BDEPO assumes g = g_x = g_y so that
      g = G / (8 * sqrt(2)).  Thus, a total scale factor of
      CPARM(3) = 1 / (8 * sqrt(2)) would be required.  This
      factor is built in as the default when CPARM(3) = 0 and
      only one input image is specified.

      If two input images are specified, CPARM(3,4) = 0 =>
      CPARM(3,4) = 1.0

Adverbs:
USERID........The ID of the owner of the images. 0=>current
              user.  32000=>all users.
INNAME........First (generalised or x-axis) RM gradient image
              name
INCLASS.......First RM gradient image class
INSEQ.........First RM gradient image seq. #
INDISK........First RM gradient image disk drive #
IN2NAME.......Second (y-axis) RM gradient image name
IN2CLASS......Second RM gradient image class
IN2SEQ........Second RM gradient image seq. #
IN2DISK.......Second RM gradient image disk drive #
BLC...........BLC of input image(s)
TRC...........TRC of input image(s)
OUTNAME.......Output predicted DPR image name
OUTCLASS......Output predicted DPR image class
OUTSEQ........Output predicted DPR image seq. #
OUTDISK.......Output predicted DPR image disk drive #
CPARM.........Frequencies (Hz) at which to calculate the DPR
              are stored in CPARM(1) and (2).  The order is
              unimportant.  CPARM(3) and CPARM(4) are scale
              factors by which the gradient images are
              multiplied.  If there is no second image,
              CPARM(3) = 0 => invoke scale factor appropriate to
              a NINER 'SOBL' gradient image.  If there is a
              second image, then both CPARM(3) and CPARM(4) = 0
              give scale factors of unity.
----------------------------------------------------------------
TASK :             BDEPO
DOCUMENTOR :       Neil Killeen, NRAO
RELATED PROGRAMS : RM, COMB, NINER, IMVIM

                                PURPOSE
                              -----------
Depolarization can be produced by differential rotation across
the beam owing to a Rotation Measure gradient.  This rotation
can occur anywhere along the line of sight.  Thus, "beam"
depolarization is distinguished from "internal" depolarization
in which the rotating material is mixed in with the emitting
material.   BDEPO tries to predict the beam depolarization
ratio (DPR) between two frequencies owing to a linear Rotation
Measure (RM) gradient across the beam. The depolarization ratio
is defined as DPR = m(F1) / m(F2) where m is the degree of
polarization at the freqeuncy F.  The convention used in this
program is that F1 < F2 so that, in the absence of noise, DPR <
1.0.  The output DPR image can be plotted (with IMVIM) against
an observed DPR image in order to asses the importance of beam
depolarization.  The gradient image(s) may be calculated from
your RM images (with NINER), or perhaps may be models of a
foreground screen that you have generated in some other way.
Note that because we generally make images with about 3 pixels
per beam, you should set the XINC and YINC adverbs in IMVIM to
2 or 3 so that the plotted points are more or less independent
of each other.  Also. if you use NINER, remember that it
operates on your image with a 3 x 3 matrix, so that there is a
loss of independence between points here as well.

BDEPO assumes that the intrinsic degree of polarization is
constant across the beam.   If there is unresolved and
significant structure in your RM gradient image, then BDEPO
will underpredict the beam depolarization.


MODE 1:
-------
Let us call the x (RA) and y (DEC) components of the RM gradient
g_x and g_y.  If we also assume that m is constant across the
beam,  then the DPR between wavelengths L_1 > L_2  is :

       --                             --                    ----
       |-g_x**2 * (L_1**4 - L_2**4)   |  [A * g_y/g_x + C]**2| |
DPR=exp|--------------------------- * |1+--------------------| |
       |        4 * ln2 * A           |     [A*B - C**2]     | |
       --                             --                    ----

where

            a**2 B_min**2 + b**2 B_maj**2
       A =  -----------------------------
                  B_maj**2 B_min**2

            a**2 B_maj**2 + b**2 B_min**2
       B =  -----------------------------
                  B_maj**2 B_min**2

            ab (B_maj**2 - b**2 B_min**2)
       C =  -----------------------------
                  B_maj**2 B_min**2

       B_maj = FWHM of the beam major axis
       B_min = FWHM of the beam minor axis
       a = cos(theta)
       b = sin(theta)
       theta = angle between beam major axis and positive X-axis
               measured positive counter clockwise, and with
               0-180 degree limits.

This is BDEPO's most general mode.   The x and y gradient images
should be input as the first and second images, respectively.
Before the DPR is calculated, g_x and g_y (as read directly from
the input images) can be multipled by factors stored in CPARM(3)
and CPARM(4), respectively.  That is,

     g_x = CPARM(3) * First_RM_gradient_image_value
     g_y = CPARM(4) * Second_RM_gradient_image_value

If these CPARMS are 0.0 on starting up BDEPO, they default to
1.0.  They might be used to account for weighting in the input
images, or conversion to radians from degrees.  After
application of the scale factors, the units of the images must
be RAD/M/M/PIXEL.

In the case that the beam is round, this expression simplifies
to :
          --                                                  --
          | - B_maj**2 * (g_x**2 + g_y**2) * (L_1**4 - L_2**4) |
DPR = exp | -------------------------------------------------- |
          |                    4 * ln2                         |
          --                                                  --


MODE 2:
-------
In the case that g = g_x = g_y and the beam is not round, we
have :

        --                               --                -- --
        | - g**2 * (L_1**4 - L_2**4)     |      [A + C]**2  |  |
DPR= exp|   -----------------------   *  |1 +  ------------ |  |
        |         4 * ln2 * A            |     [A*B - C**2] |  |
        --                               --               --  --

To invoke this mode in BDEPO, supply only the first gradient
image.  The scale factor CPARM(3) must then reduce the units to
RAD/M/M/PIXEL

     g = CPARM(3) * First_RM_gradient_image_value

If the beam is round, then :
          --                                         --
          |    - (g * B_maj)**2 * (L_1**4 - L_2**4)   |
DPR = exp |    ------------------------------------   |
          |                   2 * ln2                 |
          --                                         --

As an example, consider the output produced by NINER with opcode
'SOBL' and DPARM(1) = 0.0.   This is a generalised gradient
formed from a 3 x 3 matrix as follows :

                 | A1      A2      A3 |
                 | A8      A0      A4 |
                 | A7      A6      A5 |

NINER computes G = sqrt (X**2 + Y**2)
where     X = (A3 + 2A4 + A5) - (A1 + 2A8 + A7)
          Y = (A1 + 2A2 + A3) - (A7 + 2A6 + A5)

Now, we can associate X with g_x and Y with g_y such that

          g_x ~ X / (4*2)    and  g_y ~ Y / (4*2)

The factor of 4 for weighting, the factor of two because a
gradient = DELTA Y / DELTA X and DELTA X = 2 pixels so that

          G ~ 8 sqrt (g_x**2 + g_y**2)

Now in MODE 2, BDEPO assumes g = g_x = g_y so that

          G ~ 8 * sqrt (2) * g

Thus, the scale factor in BDEPO input via CPARM(3) that is
necessary for a NINER 'SOBL' (DPARM(1)=0)) gradient image is

          1 / (8 sqrt(2))

This is the MODE 2 default scale factor applied when CPARM(3)=0


Other parameters:

BLC, TRC     :  The usual windows for the input images.  These
                default to the entire image

CPARM(1 & 2) :  Specify the frequencies (Hz) between which to
                calculate the DPR. The order is unimportant and
                there are no defaults.
----------------------------------------------------------------

