sersic 1.0.2

Exact deprojection of Sérsic surface brightness profiles.

# -*- coding: utf-8 -*-
========
Sérsic
========

This is an implementation of the exact deprojection of Sérsic surface
brightness profiles as described in:

"Analytical expressions for the deprojected Sérsic model"
Baes and Gentile
Astronomy and Astrophysics, Volume 525:A136 (2011)
http://arxiv.org/abs/1009.4713

This code depends on the mpmath python library
(http://mpmath.googlecode.com) for an implementation of the Meijer G
function required by the Baes and Gentile (hereafter B+G) formulas for
rational values of the Sérsic index. B+G also give formulas for
irrational Sérsic indices, but I could not find an implementation of
the required Fox H function. Therefore this code requires rational
Sérsic indices, but any irrational number can be approximated
arbitrarily well by some rational number, so... no problem.

The code also depends on scipy (http://scipy.org), but the dependence
is mostly for testing: doing numerical integrals and such. If you
trust that the code passes the tests and don't want to run them
yourself (ha!), the only dependence on scipy is finding the root
(scipy.optimize.newton) in the function bg_constants(), which would be
pretty easy to replace.

This was implemented for use use in the paper:
"On Galaxies and Homology"
Novak et. al.
Monthly Notices of the Royal Astronomical Society, 424:635 (2012)
http://arxiv.org/abs/1205.2533

The B+G formulas are not simple but not terribly complex. However, I
wanted to be quite sure that I had correctly implemented the formulas
(and that the formulas themselves were correct), so I put significant
effort into comprehensive tests. That testing is the primary virtue
of this particular implementation of the B+G formulas, and is my
reason for releasing this code.

To test the implementation of the B+G formulas, I also implemented the
approximate formulas for deprojected Sérsic profiles described in:

"Dark matter in elliptical galaxies - I. Is the total mass density
profile of the NFW form or even steeper?"
Mamon and Lokas
Monthly Notices of the Royal Astronomical Society, Volume 362:95 (2005)
http://arxiv.org/abs/astro-ph/0405466

Additional tests use the parameterization of the Sérsic profile given
in:

Galactic Astronomy
Binney and Merrifield
Princeton University Press, 1998

For more information:
http://pypi.python.org/pypi/sersic/

For the source code:
http://launchpad.net/sersic

Installation
============

First install the mpmath python library (http://mpmath.googlecode.com).

To run the tests and avoid implementing your own root finding routine,
install numpy (http://www.numpy.org) and scipy (http://www.scipy.org).
But you've got those installed already, right?

The sersic module can be installed by any one of the following
standard incantations:

pip sersic
or
easy_install sersic
or
python setup.py install

If you want to install in your home directory, you can add the --user
flag to any of the above.

To run the tests:
python test.py

Usage
=====

Very likely the only function any user will care about is:

def luminosity(pp, qq, reff=1.0, lum=1.0) where pp and qq are the
numerator and denominator of the Sersic index (both integers) so that
nn=pp/qq, reff is the projected half-light radius, and luminosity is
the total luminosity. This returns a function that takes a radius and
returns a luminosity density.

>>> lum = luminosity(5,3)
>>> lum(1.1)
>>> lum([1.1, 2.2, 3.3])

All of the other functions in the module are just to help make sure
that the Baes + Gentile expressions are implemented correctly.

If you find that the calculations are taking a long time, you can try
reducing the precision of the calculation of the Meijer G function.
The mpmath library is designed to be able to work to arbitrary
precision, and is by default set to ~double precision. For
information on changing the precision, see:

http://mpmath.googlecode.com/svn/trunk/doc/build/basics.html#setting-the-precision

Keep in mind that small values of the numerator and denominator are
faster than large ones (e.g. luminosity(2,1) is fast,
luminosity(200,99) is slow.

The function returns arbitrary precision floats defined by the mpmath
library. They should behave as normal floats, but you may want to
cast them to Python floats to avoid invoking mpmath's more careful
arithmetic.

License
=======

The code is released under the MIT license, so you should be able to
do whatever you want with it.

If you use this software to produce a scientific publication or if you
incorporate this code into a larger project, I would appreciate it if
you send me a note at [e-mail]