.. C-Munipack - User's manual

   Copyright 2012 David Motl
   
   Permission is granted to copy, distribute and/or modify this document under the 
   terms of the GNU Free Documentation License, Version 1.2 or any later version published 
   by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and 
   no Back-Cover Texts.

   $Id: photometry.rst,v 1.2 2015/07/12 07:44:57 dmotl Exp $

.. _photometry: 
   
.. index:: Aperture photometry
   pair: aperture; photometry
   pair: point-spread; function

Aperture photometry
===================

Aperture photometry determines brightness of a stellar object by integrating 
the signal in a small area on a frame. This area, usually called the aperture, is 
chosen in such a way that it includes all pixels that were exposed by the light 
from the measured object, but not the light from another source. 
             
The point spread function that represents the spatial distribution of the signal 
from an ideal stellar object is rotationally symmetric, so the aperture has circular 
shape. The conditions described above define requirements of the ideal aperture size.
If it were too small, a non-negligible part of the signal would be omitted; on the
other hand, the bigger the aperture the bigger is the chance that another object appears
in the area.

Supposing the frames were properly calibrated, the frame contains the superposition 
of two sources - the sky background and the signal from objects. When we sum
the pixel values in an aperture, the sum contains these both of these. To get a 
signal from the object only, we need to subtract the background level.
 
.. index:: sky annulus

It is not possible to measure the background at the positions of the objects directly. 
We need to estimated local background level from neighboring pixels. For this
estimation we use another aperture that has an annular shape. Its center is 
aligned to the same point as the integrating aperture. The central blank part of 
the annulus masks the signal from the measured object. The mean value of pixels
in the annulus is used to estimate background level at its center. In an ideal 
case, the annulus would contain only pixels that represent the background, without 
objects. In practice this is not generally the case, therefore a robust algorithm 
for estimation of the mean level must be applied - such an algorithm must effectively
reject pixels exposed by stellar objects in the annulus. See the chapter
:ref:`robust-mean` for its implementation in the C-Munipack software.

.. figure:: images/aphot.png
   :align: center
   :alt: Aperture photometry
   
   Aperture photometry. (a) Aperture *i*: *F* -- image data samples, *I* -- net
   intensity of a stellar object, *B\,.\,A* -- background signal, *r_i* -- aperture radius;
   (b) Sky annulus: *B* -- background mean level, :math:`r_{IS}` -- inner radius of sky annulus,
   :math:`r_{OS}` -- outer radius of sky annulus.

The C-Munipack software uses the same algorithm for aperture photometry as
Stetson's DAOPHOT software. You can see [stetson11]_ for the documentation 
that comes with the original code.  

.. rubric:: Deriving brightness of an object

We have stated that the sum *S* of pixels in a small area *A* around an object
is a sum of the object's net intensity *I* plus background intensity :math:`B\,.\,A`:

.. math::
   :label: s
   
	S = I + B\,.\,A

The values of *S* and *B* are derived from the source frame, the area *A* is determined
as the area of circle of radius *r*, where *r* is the size of the aperture in pixels.
It is then easy to compute the net intensity *I* of an object in ADU:

.. math::
   :label: object_intensity
   
	I = S - B\,.\,A

.. index::
   pair: Pogson's; law	
   pair: instrumental; magnitude
	
Supposing that the net intensity *I* is proportional to the observed flux *F*, we can derive
the apparent magnitude *m* of the object, utilizing the Pogson's law:

.. math::
   :label: Pogson_law
    
	m = -2.5 \log_{10}\left(\frac{I}{I_0}\right)

.. index::
   pair: differential; magnitude
	
where :math:`I_0` is a reference intensity from an object of some chosen reference flux :math:`F_0`. The 
ratio between flux and intensity is not known, however it is legitimate to choose any reference 
intensity :math:`I_0` value, providing only the magnitude difference between two objects is required - 
the difference is independent of choice of the reference intensity :math:`I_0`. In the C-Munipack 
software, the reference intensity was set to :math:`10^{10}` ADU.

.. rubric:: When the brightness cannot be measured?

In some situations, brightness of an object cannot be measured:

* The distance between the object's center and the image border is smaller than the aperture radius.  
* There is an overexposed pixel or a bad pixel inside the aperture.
* The net intensity *I* is zero or negative.

When the output data are saved to an DAOPHOT-compatible photometry file, unmeasured
objects are indicated by brightness value 99.9999. The reader should recognize any
brightness greater than 99.0 magnitudes as an invalid measurement.     
                                                                         
.. rubric:: Estimating the measurement error

Once we have derived the raw instrumental brightness of an object, we will try
to estimate its standard error. First of all, we will recall a few general rules that apply
to the standard error and its propagation. This is a general rule for error propagation
through a function *f* of uncertain value *X*:

.. math::
   :label: var_x
   
	\operatorname{Var}(f(X)) = (\frac{df}{dx})^2 \operatorname{Var}(X)

Using this general rule, we derive two laws of error propagation. In the first case,
the uncertain value *X* is multiplied by a constant *a* and shifted by a constant
offset *b*. This law can also be used in the case where only a multiplication or only 
an offset occurs.    
                                                                                   
.. math::
   :label: var_axb
   
	\operatorname{Var}(aX + b) = a^2 \operatorname{Var}(X)

The second law defines the error of a logarithm of uncertain value *X*:   

.. math::
   :label: var_log_bx

	\operatorname{Var}(\log(\pm bX)) = \frac{\operatorname{Var}(X)}{\bar{X}^2} 

Please note, that the *log* function here is the natural logarithm, while the Pogson's 
formula (see above) incorporates the base-10 logarithm. The following equation helps
us to deal with this difference:

.. math::
   :label: log_b_x
   
	\log_b(x) = \frac{\log_k(x)}{\log_k(b)}

Putting these two equations together we get:

.. math::
   :label: var_log_dx
   
	\operatorname{Var}(\log_{10}(\pm bX)) = \frac{\operatorname{Var}(X)}{\bar{X}^2\,\log(10)^2} 
   
If we have two uncorrelated uncertain variables *X* and *Y*, the variance of their 
sum is the sum of their variances, this equation is known as Bienaym\'{e} formula. 

.. math::
   :label: sum_uncertain_variables
   
	\operatorname{Var}(X + Y) = \operatorname{Var}(X) + \operatorname{Var}(Y)

From this formula, we can also derive the standard error of a sample mean. If
we have *N* observations of random variable *X* with sample-based estimate of
the standard error of the population *s*, then the standard error of a sample 
mean estimate of the population mean is 

.. math::
   :label: se_x
   
	SE_{\bar{X}} = \frac{s}{\sqrt{N}}

Armed with this knowledge, we can start thinking about the estimation of standard
error of object brightness. We will consider the following three sources of uncertainty:
(1) random noise inside the star aperture that includes the thermal noise of the detector,
read-out noise of the signal amplifier and the analog-to-digital converter, (2) Poisson
statistics of counting of discreet events (photons incident on a detector) that
occur during a fixed period of time and (3) the error of estimation of mean sky level.

For the estimation of mean sky level, we have used the robust mean algorithm. It
allows to estimate its sample variance :math:`\sigma_{pxl}^2`. This is a pixel-based variance
and because we have summed together *A* pixels in the star aperture, the Bienaym\'{e} 
formula applies, the sum *S* is a sum of *A* uncorrelated random variables, each of
which has variance :math:`\sigma_{pxl}^2`. For the variance of the first source of error we
get:

.. math::
   :label: sigma_2_1
   
	\sigma_1^2 = A\,\sigma_{pxl}^2

where *A* is a number of pixels in the star aperture. 

From Poisson statistics we can derive a variance that occur due to counting of discreet
events, photons incident on a detector, that occur during a fixed period of time, the
exposure. We will again need to use the gain *p* of the detector to convert a intensity
in ADU to a number of photons. If the measured net intensity of an object is *I* we compute
the mean number of photons :math:`lambda` as

.. math::
   :label: lambda
	
   \lambda = I\,p

Then, the variance of intensity due to Poisson statistics is equal to its mean value.

.. math::
   :label: sigma_ph
   
	\sigma_{ph}^2 = \operatorname{Var}(\operatorname{Pois}(\lambda)) = \lambda = I\,p

The variance is in photons, we have to convert it back to ADU to get the variance
in units :math:`ADU^2`.

.. math::
   :label: sigma_22

	\sigma_2^2 = \frac{\sigma_{ph}^2}{p^2} = \frac{I\,p}{p^2} = \frac{I}{p}

We have derived the sky level as a sample mean of pixel population in the sky annulus.
Because each pixel in the annulus has variance :math:`\sigma_{pxl}^2`, the variance of sample 
mean is

.. math::
   :label: s_sky_2
   
	s_{sky}^2 = \frac{\sigma_{pxl}^2}{n_{sky}}

where :math:`n_{sky}` is the number of pixels in sky annulus.

From equation :eq:`sum_uncertain_variables` we compute the variance of object's 
intensity as

.. math::
   :label: sum_noises
   
	\sigma_{ADU}^2 = \sigma_1^2 + \sigma_2^2 + A^2\,s_{sky}^2

Note, that in equation :eq:`object_intensity` the sky level is multiplied by *A*, 
so we have to multiply its variance by :math:`A^2` - see the equation :eq:`sum_noises`.
Now, we use the law of error propagation for the logarithm adopted to match the formula
of the Pogson's law.

.. math::
   :label: Pogson_law_sigma
    
	\sigma_{mag}^2 = (\frac{-2.5}{I\,\log(10)})^2\,\sigma_{ADU}^2

Putting equations :eq:`Pogson_law_sigma` and :eq:`sum_noises` together, 
we can derive the standard error of the object's brightness in magnitudes as

.. math::
   :label: sigma_mag
   
	\sigma_{mag} = \frac{1.08574}{I}\,\sqrt{\sigma_{ADU}^2}

.. rubric:: Intensities

In some cases when further processing of the photometric data is involved, ensemble
photometry for example, it is necessary to get intensity (brightness) values of 
individual stars instead of the standard differential magnitudes. The intensity :math:`I_{ADU}` 
is defined as integral of the light in an aperture minus contribution of a local 
background. The intensity can be transformed into magnitudes using Pogson's law 
and a reference intensity. In the C-Munipack software, this value is called the 
raw instrumental magnitude.

The net intensity *I* in ADU can be computed from the raw instrumental magnitude *m* using an 
inversion of the Pogson's law (see :eq:`Pogson_law`). As was stated above,
the reference flux was set to :math:`10^{10}` ADU. The resulting relation is:

.. math:
   :label: Pogson_law_inverted
   
	I = I_0.10^{\frac{m}{-2.5}} = 10^{10}.10^{\frac{m}{-2.5}} = 10^{10 - 0.4\,m}

The error estimation :math:`\sigma_I` for the net intensity *I* can be computed from the error 
estimation :math:`\sigma_{mag}` using an inversion of the equation :eq:`Pogson_law_sigma`:

.. math::
   :label: sigma_adu
   
	\sigma_{ADU} = \frac{I}{1.08574}\,\sigma_{mag}
