Merge pull request #132 from JuanCab/fix_gain_dammit

Fixed documentation for calculate_noise and fixed bad terminology

stellarphot/differential_photometry/aij_rel_fluxes.py

@@ -16,7 +16,7 @@ def add_in_quadrature(array):
 def calc_aij_relative_flux(star_data, comp_stars,
                            in_place=True, coord_column=None,
                            star_id_column='star_id',
-                           flux_column_name='aperture_net_flux'):
+                           counts_column_name='aperture_net_counts'):
     """
     Calculate AstroImageJ-style flux ratios.
 
@@ -45,9 +45,9 @@ def calc_aij_relative_flux(star_data, comp_stars,
         If provided, use this column to match comparison stars to coordinates.
         If not provided, the coordinates are generated with SkyCoord.
 
-    flux_column_name : str,  optional
-        If provided, use this column to find flux.
-        If not provided, the column 'aperture_net_flux' is used.
+    counts_column_name : str,  optional
+        If provided, use this column to find counts.
+        If not provided, the column 'aperture_net_counts' is used.
 
     star_id_column : str,  optional
         Name of the column that provides a unique identifier for each
@@ -103,7 +103,7 @@ def calc_aij_relative_flux(star_data, comp_stars,
     # Check whether any of the comp stars have NaN values and,
     # if they do, exclude them from the comp set.
     check_for_nan = Table(data=[star_data[star_id_column].data,
-                                star_data[flux_column_name].data],
+                                star_data[counts_column_name].data],
                           names=['star_id', 'net_counts'])
     check_for_nan = check_for_nan.group_by('star_id')
     check_for_nan['good'] = ~np.isnan(check_for_nan['net_counts'])
@@ -121,7 +121,7 @@ def calc_aij_relative_flux(star_data, comp_stars,
     # Make a small table with just counts, errors and time for all of the comparison
     # stars.
 
-    comp_fluxes = star_data['date-obs', flux_column_name, error_column_name][good]
+    comp_fluxes = star_data['date-obs', counts_column_name, error_column_name][good]
     # print(np.isnan(comp_fluxes[flux_column_name]).sum(),
     #       np.isnan(comp_fluxes[error_column_name]).sum())
     # print(star_data[good][flux_column_name][np.isnan(comp_fluxes[flux_column_name])])
@@ -130,12 +130,12 @@ def calc_aij_relative_flux(star_data, comp_stars,
     # and convert to regular column...eventually
 
     comp_fluxes = comp_fluxes.group_by('date-obs')
-    comp_totals = comp_fluxes.groups.aggregate(np.sum)[flux_column_name]
-    comp_num_stars = comp_fluxes.groups.aggregate(np.count_nonzero)[flux_column_name]
+    comp_totals = comp_fluxes.groups.aggregate(np.sum)[counts_column_name]
+    comp_num_stars = comp_fluxes.groups.aggregate(np.count_nonzero)[counts_column_name]
     comp_errors = comp_fluxes.groups.aggregate(add_in_quadrature)[error_column_name]
 
-    comp_total_vector = np.ones_like(star_data[flux_column_name])
-    comp_error_vector = np.ones_like(star_data[flux_column_name])
+    comp_total_vector = np.ones_like(star_data[counts_column_name])
+    comp_error_vector = np.ones_like(star_data[counts_column_name])
 
     if len(set(comp_num_stars)) > 1:
         raise RuntimeError('Different number of stars in comparison sets')
@@ -144,9 +144,9 @@ def calc_aij_relative_flux(star_data, comp_stars,
 
     # Have to remove the flux of the star if the star is a comparison
     # star.
-    is_comp = np.zeros_like(star_data[flux_column_name])
+    is_comp = np.zeros_like(star_data[counts_column_name])
     is_comp[good] = 1
-    flux_offset = -star_data[flux_column_name] * is_comp
+    flux_offset = -star_data[counts_column_name] * is_comp
 
     # This seems a little hacky; there must be a better way
     for date_obs, comp_total, comp_error in zip(comp_fluxes.groups.keys,
@@ -155,11 +155,11 @@ def calc_aij_relative_flux(star_data, comp_stars,
         comp_total_vector[this_time] *= comp_total
         comp_error_vector[this_time] = comp_error
 
-    relative_flux = star_data[flux_column_name] / (comp_total_vector + flux_offset)
+    relative_flux = star_data[counts_column_name] / (comp_total_vector + flux_offset)
     relative_flux = relative_flux.flatten()
 
-    rel_flux_error = (star_data[flux_column_name] / comp_total_vector *
-                      np.sqrt((star_data[error_column_name] / star_data[flux_column_name])**2 +
+    rel_flux_error = (star_data[counts_column_name] / comp_total_vector *
+                      np.sqrt((star_data[error_column_name] / star_data[counts_column_name])**2 +
                               (comp_error_vector / comp_total_vector)**2
                               )
                      )

stellarphot/differential_photometry/tests/test_aij_rel_fluxes.py

@@ -51,7 +51,7 @@ def _raw_photometry_table():
                             _repeat(star_dec, n_times), _repeat(fluxes, n_times),
                             _repeat(errors, n_times),
                             _repeat(star_ids, n_times)],
-                      names=['date-obs', 'RA', 'Dec', 'aperture_net_flux',
+                      names=['date-obs', 'RA', 'Dec', 'aperture_net_counts',
                              'noise-aij', 'star_id'])
 
     return expected_flux_ratios, expected_flux_error, raw_table, raw_table[1:4]
@@ -111,22 +111,22 @@ def test_bad_comp_star(bad_thing):
         coord_bad_ra = coord_inp.ra + 3 * u.arcsecond
         input_table['RA'][-1] = coord_bad_ra.degree
     elif bad_thing == 'NaN':
-        input_table['aperture_net_flux'][-1] = np.nan
+        input_table['aperture_net_counts'][-1] = np.nan
 
     output_table = calc_aij_relative_flux(input_table, comp_star,
                                           in_place=False)
 
-    old_total_flux = comp_star['aperture_net_flux'].sum()
-    new_flux = old_total_flux - last_one['aperture_net_flux']
+    old_total_flux = comp_star['aperture_net_counts'].sum()
+    new_flux = old_total_flux - last_one['aperture_net_counts']
     # This works for target stars, i.e. those never in comparison set
     new_expected_flux = old_total_flux / new_flux * expected_flux
 
     # Oh wow, this is terrible....
     # Need to manually calculate for the only two that are still in comparison
-    new_expected_flux[1] = (comp_star['aperture_net_flux'][0] /
-                            comp_star['aperture_net_flux'][1])
-    new_expected_flux[2] = (comp_star['aperture_net_flux'][1] /
-                            comp_star['aperture_net_flux'][0])
+    new_expected_flux[1] = (comp_star['aperture_net_counts'][0] /
+                            comp_star['aperture_net_counts'][1])
+    new_expected_flux[2] = (comp_star['aperture_net_counts'][1] /
+                            comp_star['aperture_net_counts'][0])
 
     new_expected_flux[3] = expected_flux[3]
     if bad_thing == 'NaN':

stellarphot/io/aij.py

@@ -392,7 +392,7 @@ def generate_aij_table(table_name, comparison_table):
     by_source_columns = {
         'xcenter': 'X(IJ)',
         'ycenter': 'Y(IJ)',
-        'aperture_net_flux': 'Source-Sky',
+        'aperture_net_counts': 'Source-Sky',
         'aperture_area': 'N_Src_Pixels',
         'noise-aij': 'Source_Error',
         'snr': 'Source_SNR',

stellarphot/photometry.py

@@ -367,24 +367,24 @@ def add_to_photometry_table(phot, ccd, annulus, apertures, fname='',
     night = Time(ccd.header['DATE-OBS'], scale='utc')
     night.format = 'mjd'
     phot['night'] = int(np.floor(night.value - 0.5))
-    phot['aperture_net_flux'] = (phot['aperture_sum'] -
+    phot['aperture_net_counts'] = (phot['aperture_sum'] -
                                  (phot['aperture_area'] *
                                   phot['sky_per_pix_avg']))
 
     # This can happen, for example, when the object is faint
     # and centroiding is bad.
-    bad_flux = phot['aperture_net_flux'] < 0
+    bad_flux = phot['aperture_net_counts'] < 0
     all_bads = bad_flux | bad_fwhm
 
-    phot['aperture_net_flux'][all_bads] = np.nan
+    phot['aperture_net_counts'][all_bads] = np.nan
 
     if observatory_location.lower() == 'feder':
         phot['BJD'] = find_times(phot['date-obs'][0], phot['exposure'][0],
                                  ra=bjd_coords.ra, dec=bjd_coords.dec)
 
     if camera is not None:
         phot['mag_inst'] = \
-            (-2.5 * np.log10(camera.gain * phot['aperture_net_flux'].value /
+            (-2.5 * np.log10(camera.gain * phot['aperture_net_counts'].value /
                              phot['exposure'].value))
 
     metadata_to_add = ['AIRMASS', 'FILTER']
@@ -580,45 +580,42 @@ def photometry_on_directory(directory_with_images, object_of_interest,
 
     gain = camera.gain
 
+    # Compute and save noise
     noise = calculate_noise(gain=camera.gain, read_noise=camera.read_noise,
                             dark_current_per_sec=camera.dark_current,
-                            flux=all_phot['aperture_net_flux'],
+                            counts=all_phot['aperture_net_counts'],
                             sky_per_pix=all_phot['sky_per_pix_avg'].value,
                             aperture_area=all_phot['aperture_area'],
                             annulus_area=all_phot['annulus_area'],
                             exposure=all_phot['exposure'].value,
                             include_digitization=False)
+    all_phot['noise'] = noise  # Noise in electrons
+    all_phot['noise-aij'] = noise / gain  # Noise in counts
 
-    snr = gain * all_phot['aperture_net_flux'] / noise
-
-    all_phot['mag_error'] = 1.085736205 / snr
-    all_phot['noise'] = noise
-    # AstroImageJ includes a factor of gain in the noise. IMHO it is part of the
-    # flux but, for convenience, here it is
-    all_phot['noise-aij'] = noise / gain
+    # Compute and save SNR
+    snr = gain * all_phot['aperture_net_counts'] / noise
     all_phot['snr'] = snr
+    all_phot['mag_error'] = 1.085736205 / snr
 
     return all_phot
 
 
 def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0,
-                    flux=0.0, sky_per_pix=0.0,
-                    aperture_area=0, annulus_area=0,
-                    exposure=0,
-                    include_digitization=False):
+                    counts=0.0, sky_per_pix=0.0, aperture_area=0, annulus_area=0,
+                    exposure=0, include_digitization=False):
     """
     Computes the noise in a photometric measurement.
 
-    This function computes the noise (in units of gain) in a photometric measurement using the
+    This function computes the noise (in units of electrons) in a photometric measurement using the
     revised CCD equation from Collins et al (2017) AJ, 153, 77.  The equation is:
 
     .. math::
 
-        \\sigma = \\sqrt{G \\cdot F + A \\cdot \\left(1 + \\frac{A}{B}\\right)\\cdot \\left[ G\\cdot S + D \\cdot t + R^2 + (0.289 G)^2\\right]}
+        \\sigma = \\sqrt{G \\cdot C + A \\cdot \\left(1 + \\frac{A}{B}\\right)\\cdot \\left[ G\\cdot S + D \\cdot t + R^2 + (0.289 G)^2\\right]}
 
-    where :math:`\sigma` is the noise, :math:`G` is the gain, :math:`F` is the flux,
+    where :math:`\sigma` is the noise, :math:`G` is the gain, :math:`C` is the source counts,
     :math:`A` is the aperture area in pixels, :math:`B` is the annulus area in pixels,
-    :math:`S` is the sky per pixel, :math:`D` is the dark current per second,
+    :math:`S` is the sky counts per pixel, :math:`D` is the dark current per second,
     :math:`R` is the read noise, and :math:`t` is exposure time.
 
     Note: The :math:`(0.289 G)^2` term is "digitization noise" and is optional.
@@ -627,19 +624,19 @@ def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0,
     ----------
 
     gain : float, optional
-        Gain of the CCD. In units of electrons per DN.
+        Gain of the CCD. In units of electrons per count.
 
     read_noise : float, optional
         Read noise of the CCD in electrons.
 
     dark_current_per_sec : float, optional
         Dark current of the CCD in electrons per second.
 
-    flux : float, optional
-        Flux of the source in electrons.
+    counts : float, optional
+        Counts of the source.
 
     sky_per_pix : float, optional
-        Sky per pixel in electrons.
+        Sky per pixel in counts per pixel.
 
     aperture_area : int, optional
         Area of the aperture in pixels.
@@ -657,7 +654,7 @@ def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0,
     -------
 
     noise : float
-        The noise in the photometric measurement.
+        The noise in the photometric measurement in electrons.
     """
 
     try:
@@ -670,8 +667,8 @@ def calculate_noise(gain=1.0, read_noise=0.0, dark_current_per_sec=0.0,
     else:
         area_ratio = aperture_area * (1 + aperture_area / annulus_area)
 
-    # Units in here are an absolute mess
-    poisson_source = gain * flux
+    # Convert counts to electrons
+    poisson_source = gain * counts
 
     try:
         poisson_source = poisson_source.value

stellarphot/tests/test_photometry.py

@@ -21,7 +21,7 @@ def test_calc_noise_source_only(gain, aperture_area):
     expected = np.sqrt(gain * counts)
 
     np.testing.assert_allclose(calculate_noise(gain=gain,
-                                               flux=counts,
+                                               counts=counts,
                                                aperture_area=aperture_area),
                                expected)
 
@@ -103,7 +103,7 @@ def test_calc_noise_messy_case(digit, expected):
     read_noise = 12
 
     np.testing.assert_allclose(
-        calculate_noise(flux=counts,
+        calculate_noise(counts=counts,
                         gain=gain,
                         dark_current_per_sec=dark_current,
                         read_noise=read_noise,