.. role:: m(raw) :format: html latex .. role:: l(raw) :format: latex ============================================================== PyARTS User Guide, Algorithm Description and Theoretical Basis ============================================================== :author: Cory Davis .. include:: version.rst :email: cdavis@staffmail.ed.ac.uk .. contents:: :depth: 2 :l:`\newpage` Introduction ============ About this document ------------------- This document has two purposes; it aims to serve as a user guide, with descriptions of important functions and classes and examples of their use, and also where necessary, there are algorithm descriptions, and theoretical arguments justifying these algorithms. Much of the user guide components of this document have been automatically extracted from the *docstrings* in the PyARTS source code ( see http://epydoc.sourceforge.net/docstrings.html ). This ensures consistency between the user guide and on-line help, keeps the user guide up-to-date, and also improves the quality of the on-line help. Algorithm description and theoretical basis (ATBD) content has only been included in cases where the algorithm in question is novel and complex. The most mathematically arduous component of this package is the *T* -matrix code of Mishchenko. Since this code has been included in only a very slightly modified form, there is no *T* -matrix ATBD information included in this document. Instead the user is directed to the appropriate papers by Mishchenko *et al* (see `The arts_scat module`_ for exact references), which are all available from http://www.giss.nasa.gov/~crmim/publications. Why Python? ----------- The following are some of the reasons that make Python an attractive language for scientific computing - straightforward incorporation of FORTRAN code - elegant syntax - Fully object oriented - Interactive - comprehensive standard library - powerful, and freely available third-party scientific packages (SciPy, matplotlib) - platform independent - straightforward package distribution - **FREE** In the case of PyARTS, the choice of Python was primarily based on the first point, which was very important given the reliance upon pre-existing fortran code, and also the last point, which removes an important obstruction in the sharing of code between scientists. Interactivity, and the availability of high level, well documented libraries, contributed to the rapid development of the PyARTS package. .. include:: ../README.txt :l:`\newpage` .. include:: PyARTS.rst :l:`\newpage` The ``arts`` module =================== .. include:: arts.rst ``ArtsRun`` objects ------------------- .. include:: arts.ArtsRun.rst Selected methods ++++++++++++++++ .. include:: arts.ArtsRun.run.rst .. include:: arts.ArtsRun.start.rst .. include:: arts.ArtsRun.process_out_stream.rst .. include:: arts.ArtsRun.run_parallel.rst Selected Functions ------------------ .. include:: arts.pnd_fieldCalc.rst .. include:: arts.ppathCalc.rst .. include:: arts.scat_data_monoCalc.rst .. include:: arts.xml_ascii_to_binary.rst .. include:: arts.xml_binary_to_ascii.rst .. _arts.create_incoming_lookup: .. include:: arts.create_incoming_lookup.rst .. include:: arts.IWP_cloud_opt_pathCalc.rst :l:`\newpage` The ``clouds`` module ===================== .. include:: clouds.rst ``Cloud`` objects ----------------- .. include:: clouds.Cloud.rst Selected methods +++++++++++++++++ .. include:: clouds.Cloud.addHydrometeor.rst .. _Cloud.pnd_field_gen: .. include:: clouds.Cloud.pnd_field_gen.rst .. _Cloud.scat_file_gen: .. include:: clouds.Cloud.scat_file_gen.rst Hydrometeor objects ------------------- ``Crystal`` +++++++++++ .. include:: clouds.Crystal.rst ``Droplet`` +++++++++++ .. include:: clouds.Droplet.rst ``Gamma`` +++++++++ .. include:: clouds.Gamma.rst ``MonoCrystal`` +++++++++++++++ .. include:: clouds.MonoCrystal.rst Selected functions ------------------ .. include:: clouds.boxcloud.rst .. include:: clouds.gamma_dist.rst .. _clouds.mh97: .. include:: clouds.mh97.rst .. include:: clouds.nioku.rst Algorithm Description and Theoretical Basis ------------------------------------------- How Cloud and Hydrometeor objects interact ++++++++++++++++++++++++++++++++++++++++++ The Cloud class, and the various hydrometeor classes (see `Hydrometeor objects`_), aim to make it easy to construct arts scenarios from IWC/LWC fields, such as those obtained from NWP/GCM model output, and to be flexible and simple in the application of different microphysical regimes. Examination of the Cloud source code reveals a very simple class. A Cloud object contains a data member *hydrometeors* which is a list of hydrometeor objects. The `Cloud.scat_file_gen`_ and `Cloud.pnd_field_gen`_ methods simply iterate through this list, calling the *scat_calc*, and *pnd_calc* methods of each hydrometeor object. This results in a list of scattering data files, and particle number density field for the scenario. The only requirement for the hydrometeor objects is that they have the methods *scat_calc*, and *pnd_calc* methods, which have the same arguments as e.g. Crystal objects... .. include:: clouds.Crystal.scat_calc.rst .. include:: clouds.Crystal.pnd_calc.rst This system should allow for the straightforward implementation of new user defined microphysical regimes. .. figure:: clouds.png :scale: 25 The Cloud class. Using Laguerre - Gauss Quadrature to represent scattering properties of particle polydispersions ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Previously the method for calculating particle number densities (PND) has been sub-optimal. We arbitrarily chose a set of particle sizes, and took bin boundaries between them to give our size bins. The particle size distribution function was then integrated over the size bin to give the particle number density for each size. These PNDs were then all scaled so that IWC was conserved. This method is inelegant : there is no satisfactory way of determining the sizebins / bin points, which led to the choice of a large number (40) of size bins for safety, and is unnecessarily costly. The theory of Gaussian quadrature states that for an *N* point method, the approximation, .. raw:: latex html \begin{equation} \int_0^\infty x^a\exp(-x)f(x)dx\cong\sum_{i=1}^N w_if(x_i) \end{equation} , is *exact* if :m:`$f(x)$` is a polynomial of order up to :m:`$2N-1$`. The weighting function on the left is closely related to Gaussian distribution and modified Gaussian distributons often found in cloud particle size distributions. The :m:`$x^a$` term can accomodate some of the radial dependency (eg :m:`$r^2$`, :m:`$r^3$`, :m:`$r^6$`) of single scattering properties. Given that our particle number densities are used to calculate some single scattering property :m:`$\Phi$`, for a polydisperion with some size distribution function :m:`$n(r)$`, then in ARTS we will be calculating .. _Eq. 2: .. raw:: latex html \begin{eqnarray} \int_0^\infty n(r) \Phi(r) dr&=& \int_0^\infty \frac{n(r)}{x^a\exp(-x)} x^a\exp(-x)\Phi(r)\frac{dr}{dx}dx\\\nonumber &\cong&\sum_{i=1}^N \frac{w_in(r_i)}{x_i^a\exp(-x_i)}\left(\frac{dr}{dx}\right)_i\Phi(r_i) \label{eq2} \end{eqnarray} , where :m:`$r$` and :m:`$x$` are related by a simple transformation, the exact form of which is determined by the size distribution. The method will be most successful if :m:`$\frac{n(r)}{x^a\exp(-x)}\frac{dr}{dx}\Phi(r)$` can be well approximated by a polynomial. `Eq. 2`_ suggests that we represent the polydispersion using a set of *N* particles with sizes given by the Gauss-Laguerre abscissa, :m:`$x_i$`, and for each particle, :m:`$i$`, the particle number density is given by .. raw:: latex html \begin{equation} PND_i=\frac{w_in(r_i)}{x_i^a\exp(-x_i)}\left(\frac{dr}{dx}\right)_i \end{equation} Calculation of abscissas and weights for Gauss-Laguerre quadrature is done using the ``scipy`` function ``special.laguerre`` A demonstration ~~~~~~~~~~~~~~~ .. figure:: phasLaguerre.pdf The scattering matrix for liquid spheres, with a liquid water content of 1 gm\ :m:`$^{-3}$`, using a modified gamma size distribution. The :m:`$x$`- axes represent the number of quadrature points, and the different lines on each plot are for different scattering angles. Figure 1 indicates that 3 quadrature points (and hence particle types in ARTS) is sufficient for calculating the single scattering properties of liquid water clouds obeying a modified gamma distribution. Reducing the number of particles needed in ARTS simulations improves performance of both ARTS-MC and ARTS-DOIT scattering modules. Implementation in ``Droplet`` and ``Gamma`` objects ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ``Gamma`` hydrometeors, and ``Droplet`` hydrometeors, which use a modified gamma size distribution, are economical because the non-linear factors in the size distribution function are considered independent of the atmospheric field variables (IWC, LWC, and *T*). This means clouds have the same normalised size distribution at all positions, where size distribution is then scaled to give the correct LWC or IWC. This allows us to use a single *particle type*, with the scattering properties corresponding to an IWC/LWC of 1 gm\ :sup:`-3`. The PND field is then identical to the IWC/LWC field. Implementation in ``Crystal`` objects ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ``Crystal`` hydrometeor objects use the MacFarquhar and Heymsfield (1997) size distribution which was obtained from aircraft measurements in tropical cirrus. This correlation (see `clouds.mh97`_) is clearly more complicated than the exponential form best suited for Laguerre Gauss quadrature. However Laguerre Gauss quadrature still seems a good choice given MH97's use of the gamma distribution for small particles. For ARTS we require a finite, and as small as possible, number of particle types. In `Eq. 2`_ we use the transformation :m:`$x=2\alpha r$`. Since the exponential term in the MH97 gamma component depends on IWC we have to choose a suitable value for :m:`$\alpha$`, that will give accurate quadrature for the range of IWC encountered, using a minimum number of quadrature points (particle types). The likely range of values for :m:`$\alpha_{<100}$` in MH97, led to the choice of :m:`$\alpha=0.02$`. A simple test, involving calculating IWC using Laguerre Gauss quadrature for a range of input IWC, showed that :m:`$\alpha=0.02$` resulted in errors of 1\% for IWC=1 gm\ :m:`$^{-3}$` and N=4, and 1\% for IWC=0.1 gm\ :m:`$^{-3}$` and N=7. By default the Crystal class uses *N*\ =10 quadrature points (particle types). :l:`\newpage` The arts_scat module ==================== .. include:: arts_scat.rst SingleScatteringData objects ---------------------------- .. _arts_scat.SingleScatteringData: .. include:: arts_scat.SingleScatteringData.rst .. figure:: SingleScatteringData.png :scale: 25 The SingleScatteringData class. Selected methods +++++++++++++++++ .. include:: arts_scat.SingleScatteringData.__init__.rst .. include:: arts_scat.SingleScatteringData.calc.rst .. include:: arts_scat.SingleScatteringData.generate.rst .. include:: arts_scat.SingleScatteringData.load.rst .. include:: arts_scat.SingleScatteringData.save.rst Selected functions ------------------- .. include:: arts_scat.batch_generate.rst .. include:: arts_scat.combine.rst .. include:: arts_scat.refice.rst .. _arts_scat.refliquid: .. include:: arts_scat.refliquid.rst .. _tmat_fxd: .. include:: arts_scat.tmat_fxd.rst .. include:: arts_scat.tmat_rnd.rst .. include:: arts_scat.phasmat.rst .. include:: arts_scat.extmat.rst Algorithm 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 in the following sections. 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, :m:`$C_{sca}$`, the extinction cross-section, :m:`$C_{ext}$`, and the scattering matrix elements, :m:`$F_{11}$`, :m:`$F_{22}$`, :m:`$F_{33}$`, :m:`$F_{44}$`, :m:`$F_{12}$`, :m:`$F_{34}$`. This function is called by the SingleScatteringData.calc() method, and it is not intended for ``tmd.tmd`` to be called directly by the user. ``SingleScatteringData.calc()`` calls ``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 clouds module`_. However, if you want to use the size distribution capabilities of Mishchenko's code, the ``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 :m:`$\mathbf{S(n,n',\alpha,\beta)}$`, for given incident, :m:`$\mathbf{n}$`, and scattered, :m:`$\mathbf{n'}$`, directions, and orientation angles :m:`$\alpha$`, and :m:`$\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 :m:`$\mathbf{S(n,n',\alpha,\beta)}$`, it is straightforward to get the extinction matrix :m:`$\mathbf{K(n)}$`, and phase matrix :m:`$\mathbf{Z(n,n')}$`; for details see [mishchenko00]_. 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 :m:`$\beta=0$`, which makes the orientation angle :m:`$\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, :m:`$\alpha$`, with :m:`$\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 :m:`$\left<\mathbf{Z(n,n')}\right>$` must be obtained by numerical integration. .. raw:: latex html \[\left<\mathbf{Z(n,n')}\right>=\frac{1}{\pi}\int^\pi_0\mathbf{Z(n,n'},\beta=\pi/2,\alpha)d\alpha\] 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 absorption coefficient vector ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Calculation of the absorption coefficient vector is the most taxing part of the arts_scat.SingleScatteringData calculations, particularly for oblate p30 particles. For p20 particles we have simply the absorption cross-section, :m:`$K_{a1}=C_{ext}-C_{sca}$`, where the values on the RHS are obtained directly from ``tmd.tmd`` However, for p30 particles, the absorption coefficient vector is given by .. _eq. 4: .. raw:: html latex \begin{eqnarray} K_{ai}&=&\left-\int_{4\pi}d\mathbf{(n')}\left\\\nonumber &=&\left-2\int_{\pi}d\Delta\phi\int_{\pi}d\theta'\leftsin(\theta') \label{Kabs} \end{eqnarray} The integration is performed using multi-dimensional Gauss-Legendre quadrature, which is implemented as the `multi_gauss_leg`_ function in `the arts_math module`_. In the case of prolate particles, the evaluation of :m:`$\left$` requires integration over azimuthal orientation. For this reason, the integration in `Eq. 4`_ is done using 6 point Gaussian quadrature, whereas for oblate particles we use the 10 point method. Refractive Index of Ice and Liquid Water ++++++++++++++++++++++++++++++++++++++++ The calculation of single scattering properties for ice and liquid hydrometeors requires knowledge of the complex refractive index of the material in question. For both ice and liquid water the complex refractive index is a function of temperature and frequency. PyArts incorporates the fortran code, REFICE.f of Stephen Warren, Warren Wiscombe, and Bo-Cai Gao to calculate the refractive index of ice at a given frequency and temperature. This is most easily accessed by the function ``refice`` (see above) in the arts_scat module. This function looks up tables based mainly on the tabulated data of [Warrenetal84]_. REFICE.f has incorporated data published since Warren's paper, but this is not in the mm-submm range. Stephen Warren has suggested that the data of [MaetzlerWegmueller87]_, be consulted for the microwave region. This has **not** been implemented in the REFICE extension model. Figure 2 was generated by the script ``plot_refr_ind.py``, which is in the "examples" directory of the PyARTS distribution. .. figure:: ../examples/eosmlsrefice :scale: 30 Output of examples/plot_refr_ind.py For liquid cloud droplets, the complex refractive index is calculated according to the model described in the EOS-MLS Cloudy-Sky ATBD, which is based on the empirical model of [Liebeetal89]_ and [Hufford91]_. This is implemented in the `arts_scat.refliquid`_ function described above. The range of convergent size parameters and aspect ratios for ice crystal optical property generation ----------------------------------------------------------------------------------------------------- The original *T*-matrix fortran codes have been modified to call subroutines in the `LAPACK `_ library. This gives the same extended range of convergent size parameters and aspect ratios as the optional NAG enhancements described in Mishchenko's papers. Figure 3 shows the minimum integer size parameters, for a range of aspect ratios, that cause convergence failures in the PyARTS implementation of the T-matrix codes. Here the complex refractive index is given by ... >>> from PyARTS.arts_scat import * >>> T=240 #temperature >>> f=300e9 #frequency >>> m=refice(f,T) #refractive index of ice using Warren(1984) >>> print m (1.78117084503+0.00504761422053j) The convergence failure parameters shown in Figure 3 were obtained by the script ``tmat_limits.py``, which resides in the test folder of the PyARTS dsitribution. For size parameters and aspect ratios on or above the curves in Fig. 3, the *T*-matrix code will fail to converge. .. figure:: tmatlim :scale: 50 *T*-matrix convergence failure parameters for :m:`m=(1.78117084503+0.00504761422053j)` :l:`\newpage` The ``arts_types`` module ========================= .. include:: arts_types.rst ``GriddedField3`` objects ------------------------- .. include:: arts_types.GriddedField3.rst Methods +++++++ .. include:: arts_types.GriddedField3.__init__.rst .. include:: arts_types.GriddedField3.__call__.rst .. include:: arts_types.GriddedField3.expandTo3D.rst .. include:: arts_types.GriddedField3.pad.rst .. include:: arts_types.GriddedField3.load.rst .. include:: arts_types.GriddedField3.save.rst ``GasAbsLookup`` objects ------------------------ .. include:: arts_types.GasAbsLookup.rst Methods +++++++ .. include:: arts_types.GasAbsLookup.__init__.rst .. include:: arts_types.GasAbsLookup.__call__.rst .. include:: arts_types.GasAbsLookup.calc.rst .. include:: arts_types.GasAbsLookup.load.rst .. include:: arts_types.GasAbsLookup.save.rst .. include:: arts_types.GasAbsLookup.set_slidata.rst :l:`\newpage` The ``plotting`` module ======================== .. include:: plotting.rst ``SentinelMap`` objects ----------------------- .. include:: plotting.SentinelMap.rst Selected functions ------------------ .. include:: plotting.hotcoldmap.rst .. include:: plotting.myPcolor.rst .. include:: plotting.mySubplot.rst .. include:: plotting.drawCloudBox.rst .. include:: plotting.drawPpath.rst .. include:: plotting.drawSurface.rst .. include:: plotting.setDataAspectRatioByAxisPos.rst .. include:: plotting.setDataAspectRatioByFigSize.rst .. include:: plotting.shiftaxes.rst Demonstration +++++++++++++ Figure 4 shows the output of examples/geometry.py, which uses . ``plotting.drawCloudBox``, ``plotting.drawPpath``, and ``plotting.drawSurface``. .. figure:: geometry.pdf :scale: 50 the output of examples/mc_incoming_gengeometry.py, which uses . ``plotting.drawCloudBox``, ``plotting.drawPpath``, and ``plotting.drawSurface``. :l:`\newpage` The ``artsXML`` module ========================= .. include:: artsXML.rst ``XMLfile`` objects ------------------- .. include:: artsXML.XMLfile.rst Selected Methods ++++++++++++++++ .. include:: artsXML.XMLfile.addArray.rst .. include:: artsXML.XMLfile.addArrayOfArray.rst .. include:: artsXML.XMLfile.addArrayOfString.rst .. include:: artsXML.XMLfile.addString.rst .. include:: artsXML.XMLfile.addTensor.rst .. include:: artsXML.XMLfile.close.rst Selected functions ------------------ .. include:: artsXML.saveTensor.rst .. include:: artsXML.saveArray.rst .. include:: artsXML.saveString.rst .. include:: artsXML.load.rst :l:`\newpage` The ``arts_math`` module ========================= .. include:: arts_math.rst Selected functions ------------------ .. _gauss_leg: .. include:: arts_math.gauss_leg.rst .. _multi_gauss_leg: .. include:: arts_math.multi_gauss_leg.rst .. include:: arts_math.gridmerge.rst .. include:: arts_math.locate.rst .. include:: arts_math.nlogspace.rst :l:`\newpage` The ``general`` module ========================= .. include:: general.rst Selected functions ------------------ .. include:: general.multi_thread.rst .. include:: general.multi_thread2.rst .. include:: general.quickpickle.rst .. include:: general.quickunpickle.rst :l:`\newpage` The ``sli`` module ========================= .. include:: sli.rst ``SLIData2`` objects -------------------- .. include:: sli.SLIData2.rst Selected methods ++++++++++++++++ .. include:: sli.SLIData2.__init__.rst .. include:: sli.SLIData2.refine.rst .. include:: sli.SLIData2.interp.rst .. include:: sli.SLIData2.plot.rst .. include:: sli.SLIData2.load.rst .. include:: sli.SLIData2.save.rst Demonstration +++++++++++++ Figure 5 shows the output of examples/mc_incoming_gen.py, which makes use of the ``SLIData2 class`` to create an optimised 2D (altitude, zenith angle) lookup table of incoming clear-sky radiance, for MonteCarlo scattering simulations with a pseudo-3D atmosphere (3D cloud, but a 1D atmosphere outside the cloudbox). .. figure:: sli.jpg :scale: 25 the output of examples/mc_incoming_gen.py, which makes use of the ``SLIData2 class``. :l:`\newpage` The ``arts1`` module ==================== .. include:: arts1.rst Selected functions ------------------ .. include:: arts1.abs_file.rst .. include:: arts1.arts1_file_from_command_list.rst .. include:: arts1.arts1_load.rst .. include:: arts1.arts1_save.rst :l:`\newpage` .. [Changetal97] Chang, J. Z., J. P. Allebach, C. A. Bouman, Sequential Linear Interpolation of Multidimensional Functions, *IEEE Trans. Image Proc*, 6(9),1997. (`pdf `__) .. [Hufford91] Hufford, G., A model for the complex permittivity of ice at frequencies below 1 THz. *Int. J. Infrared Millimeter Waves*, 12, 677-682, 1991. .. [Liebeetal89] Liebe, H. J., T. Manabe, and G. A. Hufford, Millimeter-wave attenuation and delay rates due to fog/cloud conditions. *IEEE Trans. Ant. Prop.*, 37, 1617-1623, 1989. .. [MaetzlerWegmueller87] Maetzler and Wegmueller, *J. Phys. D*. 20, 1623-1630, 1987 .. [McFarquharHeymsfield97] G.M. McFarquhar and A.J. Heymsfield, Parametrization of tropical ice crystal size distributions and implications for radiative transfer: Results from CEPEX, *J. Atmos. Sci.*, 54, 2187-2200, 1997. .. [mishchenko91] M. I. Mishchenko, Extinction and polarization of transmitted light by partially aligned nonspherical grains, *Astrophysical Journal*, 367, 561-574, 1991. (`pdf `__) .. [mishtrav98] Mishchenko, M.I. and Travis, L.D, Capabilities and limitations of a current FORTRAN implementation of the *T*-matrix method for randomly oriented, rotationally symmetric scatterers., *J. Quant. Spectrosc. Radiat. Transfer*, 60(3), 309-324,1998. (`pdf `__) .. [mishchenko00] M.I. Mishchenko, Calculation of the amplitude matrix for a nonspherical particle in a fixed orientation, *Applied Optics*, 39(6), 1026-1031,2000. (`pdf `__) .. [pressetal92] H.P Press, S.A Teukolksky, W.T. vetterling, and B.P Flannery, *Numerical Recipes in C: The Art of Scientific Computing*, Cambridge University Press, 1992. .. [Warrenetal84] S.G. Warren, Optical constants of ice from the ultraviolet to the microwave, *Applied Optics*, 23(8), 218-229, 1984.