=============================
The arts\_scat Python module
=============================

:author: Cory Davis

.. contents::


The SingleScatteringData class
==============================

The main component of the arts_scat module is the SingleScatteringData
class. The data members of this object are identical to the class of
the same name in ARTS; it includes all the single scattering
properties required for polarized radiative transfer calculations: the extinction
matrix, the phase matrix, and the absorption coefficient vector.  The
angular, frequency, and temperature grids for which these are defined
are also included.  Another data member - *ptype*, describes the
orientational symmetry of the particle ensemble, which determines the format of
the single scattering properties.  The data structure of the ARTS
SingleScatteringData class is described on pages 80-83 of the ARTS User
Guide.

The methods in the arts\_scat SingleScatteringData class enable the
calculation of the single scattering properties, and the output of the
SingleScatteringData structure in the ARTS XML format (see example file).  The main
methods are:-

 * **calc** (*precision* =0.001): calculates the phase matrix,
   extinction matrix, and absorption coefficient vector. The optional argument
   *precision* determines the precision of internal *T*-matrix and
   numerical integration (only for ptype = 30) computations.  The
   precision parameter is recorded  in the output XML file.

 * **file\_gen** (*filename*): outputs the single scattering data to
   *filename* in ARTS XML format

 * **generate** (*precision* = 0.001): performs *calc()* and *file\_gen*
   with a filename generated from particle parameters

A SingleScatteringData object is initialised with a dictionary of
\{*keyword*:value\} parameters.  If this dictionary is omitted, or if
any required parameters are missing from the dictionary, default
parameters are used. *This is a little dangerous - I might disable
the defaults*.  After initialisation these parameters (eg:
*equiv\_radius*, *aspect\_ratio*,...) can be changed by modifying
**SingleScatteringData.params**, although I prefer to create a new
SingleScatteringData object for each set of parameters.

The SingleScatteringData class also has a **load(filename)**
method, which loads a SingleScatteringData object from an existing
file. Note that this can only import data members that are actually in
the file - so the scattering properties won't be consistent with the
*params* data member. This method is useful when combined with
the **combine** and **compress** functions described below.

Examples of Use
---------------

1. Simple generation of single scattering data from a python
   interpreter session.


First import the required modules

::

  bash-2.05b$ python
  Python 2.2.2 (#1, Feb 24 2003, 19:13:11)
  [GCC 3.2.2 20030222 (Red Hat Linux 3.2.2-4)] on linux2
  Type "help", "copyright", "credits" or "license" for more information.
  >>> from Numeric import * #import the numeric python module
  >>> import arts_scat      #import the arts_scat module
  \end{verbatim}%$

and define the necessary parameters

::

 >>> scat_params={
 ... "NP":-1,    #particle type as in Mischenko's T-matrix code;
 ...                       #-1=spheroids,-2=cylinders
 ... 'phase':'ice'         #The phase of the scattering particle 'ice'/'liquid'
 ... 'ptype': 20,          #random orientation (30=horizontal)
 ... 'equiv_radius': 200,  #equivalent volume radius in microns
 ... 'aspect_ratio': 3,    #oblate (<1=prolate)
 ... "T_grid":[250],       #temperature is needed to calculate the
 ...                 #complex refractive index of ice/liquid water. 
 ...                 #Currently only 
 ...                 #one value is supported in ARTS.
 ... 'za_grid':arange(0,181,10),    #zenith angle grid
 ... 'aa_grid':arange(0,181,10),    #azimuth angle grid
 ... 'f_grid': [240e9, 242e9]       #frequency grid
 ... }


Now create a SingleScatteringData object, calculate scattering
properties, and write an ARTS XML file.

::

 >>> my_scat_data=arts_scat.SingleScatteringData(scat_params)
 >>> my_scat_data.calc() #use the default precision value
 >>> my_scat_data.file_gen("an_example_file.xml")

Now try the *generate()* method and see where the data ended up

>>> my_scat_data.generate()
>>> print my_scat_data.filename
/home/cory/data/p20f240T250r200NP-1ar3ice.xml

Extras
======

The arts\_scat module also includes some useful functions for
manipulating and testing SingleScatteringData objects.

combine and compress
--------------------

If for a given 1D profile you have more particle types than non-zero
ice water content (IWC) levels, it makes sense to create one
artificial SingleScatteringData object for each of these levels.  For example,
in the JPL code there are 40 size-bins, and  profile 2101 in
the 1996 simulated data set, which is a very big cloud, there are only
7 non-zero IWC levels.  In ARTS radiative transfer calculations, the
reduction from 40 to 7 particle types gives a considerable reduction
in memory use *and* CPU time.  The *combine* and *compress*
functions in the arts\_scat module enable this. 

Tests
-----

The arts\_scat module includes several tests, which, along with Pythons
unittest module, can demonstrate that everything is working properly
and also test the quality of the single scattering properties.  The
entire suite of tests can either be run from the python interpreter by
typing \`arts\_scat.run\_tests()', or from the shell with the command
\`python [python library path]/unittest.py arts\_scat.test\_suite'.

Algorithm Description and Theoretical Basis
===========================================

Single Scattering Properties
----------------------------

The single scattering data is calculated using the *T*-matrix method.
For details of the *T*-matrix method, please consult the references
mentioned below.

ptype=20
~~~~~~~~

In the ptype=20 case, where we have completely random orientation,
we use a slightly modified version of Mishchenko's random orientation
T-matrix code [mishtrav98]_.  Mischenko's source code has been
altered to provide a python extension module,
**tmd.so**, which contains the function *tmd*.  This function
returns the scattering cross-section, `C_{sca}`, the extinction
cross-section, `C_{ext}`, and the scattering matrix elements,
`F_{11}`, `F_{22}`, `F_{33}`, `F_{44}`, `F_{12}`, `F_{34}`. This
function is called by the SingleScatteringData.calc() method, and it
is not intended for
{\tt tmd.tmd } to be called directly by the user.
SingleScatteringData.calc() calls {\tt tmd.tmd} with arguments suitable for a monodispersion of particle sizes, as
this allows a consistent interface for both ptype=20 and
ptype=30. Size distributions of both these particle types can be
handled by the {\tt clouds}
module.  However, if you want to use the size distribution capabilities
of Mishchenko's code the {\tt tmd.tmd} function can be called directly.
See the pydoc documentation for the argument list and
[mishtrav98]_ for the argument definitions.

ptype=30
~~~~~~~~~~~~~~~~~~~~~~~~

For ptype=30, where there is horizontal alignment, but random
azimuthal orientation, we use a modified version the fixed orientation
code of Mishchenko [mishchenko00]_.  Mischenko's source code has been
altered to provide a python extension module,
**tmatrix.so**, which contains the functions *tmatrix* and
*ampld*.  *tmatrix.tmatrix* calculates the *T*-matrix for
given particle and radiation parameters, and
stores this in the data structure
*tmatrix.tmat*. The function *tmatrix.ampl*, uses the
`T`-matrix data to calculate the scattering amplitude function
`\mathbf{S(n,n',\alpha,\beta)}`, for given incident (`\mathbf{n}`) and
scattered (`\mathbf{n'}`) directions, and orientation angles `\alpha`,
and `\beta`.  This arrangement makes use of the fact that the
`T`-matrix need only be calculated once for a given particle and
frequency, to get single scattering properties for a range of
directions and particle orientations.  From
`\mathbf{S(n,n',\alpha,\beta)}`, it is straightforward to get the
extinction matrix `\mathbf{K(n)}`, and phase matrix
`\mathbf{Z(n,n')}`; for details see \cite{mishchenkoetal00}.  Again,
the *tmatrix*
module is not intended for the user; the functions mentioned above are
called within SingleScatteringData.calc().

The purpose of the ptype=30 case in the SingleScatteringData python
class is to represent scattering by *horizontally aligned*
particles.  This means that oblate particles (aspect ratio > 1) have their rotation axis
parallel to the local zenith.  In the notation of [mishchenko00]_,
this corresponds to an orientation angle `\beta=0`, which makes the
orientation angle `\alpha` irrelevant due to the rotational symmetry
of the particle. Conversely, prolate particles have the
axis of rotation perpendicular to the local zenith.  This means that in
the case of horizontally aligned prolate particles, scattering
properties must be averaged over all possible azimuth orientations,
`\alpha`, with `\beta=\pi/2`. 

Orientation Averaging
---------------------

This section only applies to horizontally aligned prolate
particles. The orientationally averaged extinction matrix is obtained
from the averaged *T*-matrix,
which can be calculated \`exactly' from a single *T*-matrix calculation
according to the analytic method described in [mishchenko91]_.
This is implemented in the function *tmatrix.avgTmatrix*.
Unfortunately the orientationally average *T*-matrix is not useful for
calculating the orientationally averaged phase matrix.  In short this
is because unlike the extinction matrix, phase matrix elements can not
be expressed as linear expansions of *T*-matrix elements.  Therefore
`\left<\mathbf{Z(n,n')}\right>` must be obtained by numerical
integration.


.. texmath:: \left<\mathbf{Z(n,n')}\right>=\frac{1}{\pi}\int^\pi_0\mathbf{Z(n,n'},\beta=\pi/2,\alpha)d\alpha\label{orient_av}


Several quadrature routines have been trialed for this integration. To
date, by far the best in terms of accuracy and speed has been Gauss
Legendre quadrature [pressetal92]_.  In this case we use a 10
point quadrature.  This is implemented by the
*gauss\_leg* function in the *arts\_math* module.

Calculation of the absorbtion coefficient vector
------------------------------------------------

FIXME

Refractive Index of Ice and Liquid Water
----------------------------------------

FIXME

References
==========

.. [mishtrav98] need to do something about references
.. [mishchenko00] need to do something about references
.. [mishchenko91] need to do something about references
.. [pressetal92] need to do something about references