Adding STMAG and AB mag to flux conversion

.travis.yml

@@ -27,7 +27,8 @@ before_install:
 - ls
 - conda update --yes conda
 install:
-- conda install --yes python=$TRAVIS_PYTHON_VERSION numpy scipy nose matplotlib
+- conda install --yes python=$TRAVIS_PYTHON_VERSION numpy scipy nose matplotlib astropy
+- conda install --yes -c astropy astroquery
 - pip install codecov pypandoc
 - python setup.py install
 script:

aotools/astronomy/_astronomy.py

@@ -1,5 +1,11 @@
 import numpy
+import astropy.constants as c
+from ftplib import FTP
+from io import BytesIO
+from astroquery.vizier import Vizier
+from .intrinsic_colors import colors_dict, colors_columns
 
+# Source : https://www.cfa.harvard.edu/~dfabricant/huchra/ay145/mags.html
 # Dictionary of flux values at the top of the atmosphere
 #                 band, lamda, dLamda, m=0 flux (Jy)
 FLUX_DICTIONARY = {'U': [0.36, 0.15, 1810],
@@ -45,6 +51,105 @@ def photons_per_mag(mag, mask, pixel_scale, wvlBand, exposure_time):
     return photons
 
 
+def flux_AB_mag(mag, wvl, wvlBand):
+    """
+    Calculates the photon flux for a given magnitude in a given wave band using AB definition
+
+    Parameters:
+        mag (float): magnitude in AB system
+        wvl (float): central wavelength in any length unit
+        wvlBand (float): length of wavelength band in the same unit as wvl
+
+    Returns:
+        float: flux of photons (ph/m2/s)
+    """
+    return 10 ** (-mag / 2.5) * 10 ** (-48.60/2.5 - 7 + 4) * 1/c.h.value * wvlBand / wvl
+
+
+def flux_STMAG_mag(mag, wvl, wvlBand):
+    """
+    Calculates the photon flux for a given magnitude in a given wave band using STMAG definition
+
+    Parameters:
+        mag (float): magnitude in AB system
+        wvl (float): central wavelength in nanometer
+        wvlBand (float): length of wavelength band in nanometres
+    Returns:
+        float: flux of photons (ph/m2/s)
+    """
+    return 10 ** (-mag / 2.5) * 10 ** (-21.10/2.5 - 7 + 4 + 8) * c.h.value * c.c.value * wvlBand * wvl
+
+
+def get_magV(star_type, mag, band):
+    """
+    Translates magnitude in the specified band in V band
+
+    Parameters:
+        star_type (str): only the first letter and number
+        mag (float): magnitude in the band defined by "band"
+        band (str): band in which the magnitude is expressed, can be B, J, H, K, L, M, N, R, I
+
+    Returns:
+        float: magnitude in V band
+    """
+    if band == 'B':
+        magV = mag - colors_dict[star_type[:2] + '.0'][1]
+    else:
+        ind_col = numpy.where(numpy.asarray(colors_columns) == '(V-' + band + ')')[0][0]
+        magV = colors_dict[star_type[:2] + '.0'][ind_col] + mag
+    return magV
+
+
+def flux_star_type(mag, band, star_type, wvl, wvlBand):
+    """
+    Calculates the flux over a custom wave band, for a given star of a given magnitude in one of the standard bands (I, H, J, etc.)
+
+    Parameters:
+        mag (float): magnitude in the band defined by "band"
+        band (str): band in which the magnitude is expressed, can be B, J, H, K, L, M, N, R, I
+        star_type (str): must be the full star type, i.e. "A0V"
+        wvl (float): central wavelength in microns
+        wvlBand (float): length of wavelength band in microns
+
+    Returns:
+        float: flux of photons (ph/m2/s)
+    """
+    # get star spectrum, will correspond to a V=0
+    Vizier.ROW_LIMIT = -1
+    catalog = Vizier.get_catalogs('J/PASP/110/863/synphot')[0]
+    pickles_index = numpy.where(catalog["SpType"] == star_type)[0][0] + 1
+
+    ftp = FTP('ftp.stsci.edu')
+    ftp.login()
+
+    ss = BytesIO()
+    message = ftp.retrbinary('RETR cdbs/grid/pickles/dat_uvk/pickles_uk_' + str(pickles_index) +
+                             '.ascii', ss.write)
+    ftp.quit()
+
+    # spectrum is in erg/s/cm2/A, first column are wavelengths in A, second column is the spectrum
+    erg_spec = numpy.asarray([numpy.asarray(row.split('  '), dtype='float')
+                              for row in str(ss.getvalue()).split("\\n")[39:-1]], dtype='float')
+    wavelengths = erg_spec[:, 0] * 1e-4
+
+    # translate spectrum in photon/s/m2/micron
+    phot_spec = erg_spec[:, 1] * wavelengths / (c.h.value * c.c.value) * 1e-5
+
+    # translate magnitude of band in V
+    magV = get_magV(star_type, mag, band)
+
+    # normalise spectrum to correspond to that magnitude
+    fac = 10 ** (-magV / 2.5)
+    norm_spec = phot_spec * fac
+
+    # integrate over band
+    ind = numpy.where((wavelengths >= wvl - wvlBand/2) & (wavelengths <= wvl + wvlBand/2))
+    selection = norm_spec[ind]
+    flux = (numpy.sum(selection) - (selection[0] + selection[-1]) / 2) * \
+           numpy.mean(numpy.diff(wavelengths[ind]))
+    return flux
+
+
 def photons_per_band(mag, mask, pxlScale, expTime, waveband='V'):
         '''
         Calculates the photon flux for a given aperture, star magnitude and wavelength band
@@ -88,7 +193,7 @@ def magnitude_to_flux(magnitude, waveband='V'):
     """
 
     flux_Jy = FLUX_DICTIONARY[waveband][2] * 10 ** (-0.4 * magnitude)
-    flux_photons = flux_Jy * 1.51E7 * FLUX_DICTIONARY[waveband][1]  # photons sec^-1 m^-2
+    flux_photons = flux_Jy * 1.51E7 * FLUX_DICTIONARY[waveband][1] / FLUX_DICTIONARY[waveband][0]  # photons sec^-1 m^-2
     return flux_photons
 
 
@@ -104,6 +209,6 @@ def flux_to_magnitude(flux, waveband='V'):
         float: Apparent magnitude
     """
 
-    flux_Jy = flux / (1.51E7 * FLUX_DICTIONARY[waveband][1])
+    flux_Jy = flux / (1.51E7 * FLUX_DICTIONARY[waveband][1] / FLUX_DICTIONARY[waveband][0])
     magnitude = float(-2.5 * numpy.log10(flux_Jy / FLUX_DICTIONARY[waveband][2]))
     return magnitude

aotools/astronomy/intrinsic_colors.py

@@ -0,0 +1,55 @@
+
+
+colors_columns = ['U-B','(B-V)','(V-J)','(V-H)','(V-K)','(V-L)','(V-M)','(V-N)','(V-R)','(V-I)']
+colors_dict = {
+'B0.0': [-1.08,-0.3,-0.8,-0.92,-0.97,-1.13,-1,-9.99,-0.22,-0.44],
+'B0.5': [-1,-0.28,-0.77,-0.89,-0.95,-1.11,-0.99,-9.99,-0.2,-0.43],
+'B1.0': [-0.95,-0.26,-0.73,-0.85,-0.93,-1.08,-0.96,-9.99,-0.18,-0.42],
+'B1.5': [-0.88,-0.25,-0.7,-0.82,-0.91,-1.05,-0.94,-9.99,-0.17,-0.41],
+'B2.0': [-0.81,-0.24,-0.67,-0.79,-0.89,-1.02,-0.92,-9.99,-0.15,-0.4],
+'B2.5': [-0.72,-0.22,-0.64,-0.76,-0.86,-0.97,-0.88,-0.96,-0.14,-0.39],
+'B3.0': [-0.68,-0.2,-0.6,-0.72,-0.82,-0.92,-0.84,-0.91,-0.13,-0.38],
+'B3.5': [-0.65,-0.19,-0.58,-0.7,-0.8,-0.9,-0.82,-0.87,-0.12,-0.37],
+'B4.0': [-0.63,-0.18,-0.56,-0.68,-0.77,-0.86,-0.79,-0.84,-0.11,-0.35],
+'B4.5': [-0.61,-0.17,-0.54,-0.65,-0.74,-0.83,-0.76,-0.8,-0.11,-0.34],
+'B5.0': [-0.58,-0.16,-0.51,-0.62,-0.71,-0.78,-0.73,-0.75,-0.1,-0.33],
+'B6.0': [-0.49,-0.14,-0.46,-0.57,-0.64,-0.7,-0.65,-0.66,-0.09,-0.29],
+'B7.0': [-0.43,-0.13,-0.41,-0.51,-0.57,-0.61,-0.58,-0.58,-0.08,-0.26],
+'B7.5': [-0.4,-0.12,-0.39,-0.48,-0.54,-0.57,-0.54,-0.53,-0.08,-0.24],
+'B8.0': [-0.36,-0.11,-0.36,-0.45,-0.49,-0.52,-0.49,-0.48,-0.07,-0.22],
+'B8.5': [-0.27,-0.09,-0.31,-0.4,-0.43,-0.43,-0.42,-0.39,-0.07,-0.18],
+'B9.0': [-0.18,-0.07,-0.26,-0.34,-0.33,-0.34,-0.34,-0.3,-0.06,-0.14],
+'B9.5': [-0.1,-0.04,-0.22,-0.29,-0.26,-0.27,-0.26,-0.22,-0.03,-0.11],
+'A0.0': [-0.02,-0.01,-0.16,-0.19,-0.17,-0.18,-0.18,-0.14,-0.01,-0.05],
+'A1.0': [0.01,0.02,-0.11,-0.12,-0.11,-0.12,-0.13,-0.08,0.01,-0.03],
+'A2.0': [0.05,0.05,-0.07,-0.04,-0.05,-0.07,-0.08,-0.02,0.03,0],
+'A3.0': [0.08,0.08,-0.02,0.03,0.01,-0.01,-0.02,0.03,0.05,0.02],
+'A4.0': [0.09,0.12,0.03,0.11,0.08,0.05,-0.04,0.09,0.07,0.07],
+'A5.0': [0.09,0.15,0.09,0.19,0.15,0.12,0.1,0.16,0.1,0.12],
+'A6.0': [0.1,0.17,0.13,0.3,0.21,0.17,0.15,0.21,0.11,0.16],
+'A7.0': [0.1,0.2,0.18,0.32,0.27,0.23,0.2,0.26,0.13,0.19],
+'A8.0': [0.09,0.27,0.25,0.42,0.36,0.33,0.29,0.34,0.16,0.26],
+'A9.0': [0.08,0.3,0.31,0.49,0.44,0.41,0.36,0.41,0.18,0.31],
+'F0.0': [0.03,0.32,0.37,0.57,0.52,0.49,0.43,0.48,0.21,0.36],
+'F1.0': [0,0.34,0.43,0.64,0.58,0.57,0.49,0.54,0.23,0.4],
+'F2.0': [0,0.35,0.48,0.71,0.66,0.66,0.56,0.6,0.25,0.45],
+'F5.0': [-0.02,0.45,0.67,0.93,0.89,0.9,0.77,0.8,0.33,0.57],
+'F8.0': [0.02,0.53,0.79,1.06,1.03,1.06,0.91,0.91,0.37,0.64],
+'G0.0': [0.06,0.6,0.87,1.15,1.14,1.18,1.01,1.01,0.41,0.7],
+'G2.0': [0.09,0.63,0.97,1.25,1.26,1.31,1.12,1.11,0.45,0.75],
+'G3.0': [0.12,0.65,0.98,1.27,1.28,1.33,1.14,1.13,0.45,0.76],
+'G5.0': [0.2,0.68,1.02,1.31,1.32,1.38,1.18,1.17,0.47,0.78],
+'G8.0': [0.3,0.74,1.14,1.44,1.47,1.55,1.34,1.3,0.52,0.85],
+'K0.0': [0.44,0.81,1.34,1.67,1.74,1.85,1.61,1.54,0.61,0.97],
+'K1.0': [0.48,0.86,1.46,1.8,1.89,2.02,1.78,1.68,0.67,1.05],
+'K2.0': [0.67,0.92,1.6,1.94,2.06,2.21,1.97,1.84,0.73,1.14],
+'K3.0': [0.73,0.95,1.73,2.09,2.23,2.4,2.17,2.01,0.8,1.25],
+'K4.0': [1,1,1.84,2.22,2.38,2.57,2.36,2.15,0.86,1.34],
+'K5.0': [1.06,1.15,2.04,2.46,2.66,2.87,2.71,2.44,0.97,1.54],
+'K7.0': [1.21,1.33,2.3,2.78,3.01,3.25,3.21,2.83,1.13,1.86],
+'M0.0': [1.23,1.37,2.49,3.04,3.29,3.54,3.65,3.16,1.26,2.15],
+'M1.0': [1.18,1.47,2.61,3.22,3.47,3.72,3.95,3.39,1.36,2.36],
+'M2.0': [1.15,1.47,2.74,3.42,3.67,3.92,4.31,3.66,1.46,2.62],
+'M3.0': [1.17,1.5,2.84,3.58,3.83,4.08,4.62,3.89,1.56,2.84],
+'M4.0': [1.07,1.52,2.93,3.74,3.98,4.22,4.93,4.11,1.65,3.07]
+}

appveyor.yml

@@ -21,7 +21,7 @@ install:
   - conda config --set always_yes yes --set changeps1 no
   - conda update -q conda
   - conda info -a
-  - "conda create -q -n test-environment python=%PYTHON_VERSION% numpy scipy nose matplotlib"
+  - "conda create -q -n test-environment python=%PYTHON_VERSION% numpy scipy nose matplotlib astropy"
   - activate test-environment
 
 test_script:

setup.py

@@ -23,6 +23,8 @@
     install_requires=[
         'numpy',
         'scipy',
+        'astropy',
+        'astroquery',
     ],
     classifiers=[
         "Programming Language :: Python",

test/test_astronomy.py

@@ -26,3 +26,25 @@ def test_magnitude_to_flux_and_flux_to_magnitude():
 def test_photons_per_band():
     photons = astronomy.photons_per_band(5.56, circle(2, 5), 0.5, 0.001)
     assert type(photons) == float
+
+
+def test_flux_AB_mag():
+    flux = astronomy.flux_AB_mag(10, 0.91, 0.13)
+    assert type(flux) == float
+
+
+def test_flux_STMAG_mag():
+    flux = astronomy.flux_STMAG_mag(10, 0.91, 0.13)
+    assert type(flux) == float
+
+
+def test_get_magV():
+    magV = astronomy.get_magV('A0', 0, 'H')
+    assert type(magV) == float
+    assert magV == -0.19
+
+
+def test_flux_star_type():
+    flux = astronomy.flux_star_type(10, 'H', 'A0V', 1.1, 0.2)
+    assert type(flux) == astronomy.numpy.float64
+