Source code for slugpy.photometry_convert

"""
Function to convert photometric data between the photometric systems
that SLUG knows about.
"""

import os
import os.path as osp
import numpy as np
import scipy.constants as physcons
from .int_tabulated import int_tabulated
from .int_tabulated import int_tabulated2
from .read_filter import read_filter

# Constants and units conversions to cm
try:
    c = physcons.c*100.0
except:
    # This exception is here to deal with readthedocs not having scipy
    c = 3.0e10
Angstrom = 1e-8
pc = 3.0856775814671918e18

[docs]def photometry_convert(photsystem, phot, units, wl_cen=None, filter_last=False, filter_names=None, filter_dir=None): """ Function to convert photometric data between photometric systems. Parameters photsystem : string The photometric system to which to convert. Allowable values are 'L_nu', 'L_lambda', 'AB', 'STMAG', and 'Vega', corresponding to the options defined in the SLUG code. If this is set and the conversion requested involves a conversion from a wavelength-based system to a frequency-based one, wl_cen must not be None. phot : array array of photometric data; if the array has more than one dimension, the first dimension is assumed to represent the different photometric filters (unless filter_last is True, in which case the last dimension is represents the array of filters) units : iterable of strings iterable listing the units of the input photometric data. On return, strings will be changed to the units of the new system. wl_cen : array central wavelengths of the filters, in Angstrom; can be left as None if the requested conversion doesn't require going between wavelength- and frequency-based systems. filter_last : bool If the input data have more than one dimension, by default it is assumed that the first dimension contains values for the different photometric filters. If this keyword is set to True, it will instead be assumed that the last dimension contains the values for the different filters. filter_names : iterable of strings Names of all filters, used to read the filter response functions from disk; only needed for conversions to and from Vega magnitudes, and ignored otherwise filter_dir : string Directory where the filter data files can be found. If left as None, filters will be looked for in the $SLUG_DIR/lib/filters directory. This parameter is used only for conversions to and from Vega magnitudes. Returns Nothing Raises ValueError, if wl_cen is None but the requested conversion requires going between wavelength- and frequency-based systems """ # Set target units based on requested conversion if photsystem == 'L_lambda': target_units = 'erg/s/A' elif photsystem == 'L_nu': target_units = 'erg/s/Hz' elif photsystem == 'AB': target_units = 'AB mag' elif photsystem == 'STMAG': target_units = 'ST mag' elif photsystem == 'Vega': target_units = 'Vega mag' else: raise ValueError('Unknown photometric system ' + repr(photsystem)) # If Vega magnitudes is either the current unit or the desination # unit, compute the magnitude of Vega in each of our filters if photsystem == 'Vega' or 'Vega mag' in units: # Make sure we have filter names if filter_names is None: raise ValueError("must set filter_names keyword for " "conversions to or from Vega mag") # Look in SLUG_DIR for Vega spectrum from astropy.io import ascii if 'SLUG_DIR' in os.environ: slugdir = os.environ['SLUG_DIR'] fname = osp.join(slugdir, 'lib', 'atmospheres', 'A0V_KURUCZ_92.SED') try: vegadata = ascii.read(fname) except IOError: vegadata = None # Look in cwd if that fails if vegadata is None: try: vegadata = ascii.read('A0V_KURUCZ_92.SED', 'r') except IOError: raise IOError(1, "could not find Vega spectrum "+ "file A0V_KURUCZ_92.SED in "+ "SLUG_DIR/lib/atmospheres or "+ "in current directory") # Record Vega data vega_logwl = np.log(np.array(vegadata['col1'])) vega_spec = np.array(vegadata['col2']) vega_F_nu = Angstrom * vega_spec * np.exp(2.0*vega_logwl) / c # Get the convolution of Vega's spectrum with the Johnson V # filter, and convert that to magnitudes fdat = read_filter('Johnson_V', filter_dir=filter_dir) fnorm = int_tabulated(np.log(fdat.wl), fdat.response) vega_Fbar_nu \ = int_tabulated2(vega_logwl, vega_F_nu, np.log(fdat.wl), fdat.response) / fnorm vega_mag_V = -2.5*np.log10(vega_Fbar_nu) # Now loop over other filters, and compute the magnitude of # Vega in each one, normalizing to give a magnitude of 0.02 in # V vega_mag = np.zeros(len(filter_names)) for i in range(len(filter_names)): # Skip special values if (units[i] == 'Lsun') or (units[i] == 'phot/s'): continue fdat = read_filter(filter_names[i], filter_dir=filter_dir) fnorm = int_tabulated(np.log(fdat.wl), fdat.response) vega_Fbar_nu_filter \ = int_tabulated2(vega_logwl, vega_F_nu, np.log(fdat.wl), fdat.response) / fnorm vega_mag_filter = -2.5*np.log10(vega_Fbar_nu_filter) vega_mag[i] = vega_mag_filter - vega_mag_V + 0.02 # If our input data are in Vega mag, convert them to AB mag first for i in range(len(units)): if units[i] == 'Vega mag': if filter_last: phot[:,i] = phot[:,i] + vega_mag[i] else: phot[i,:] = phot[i,:] + vega_mag[i] units[i] = 'AB mag' # Loop over photometric values for i in range(len(units)): # If units are Lsun or phot/s, this is a special value # that should not be converted. Also do nothing if target # units match current units. if (units[i] == 'Lsun') or (units[i] == 'phot/s') or \ (units[i] == target_units): continue # Do unit conversion try: # L_lambda to other units if (units[i] == 'erg/s/A') and \ (target_units == 'erg/s/Hz'): if filter_last: phot[:,i] = phot[:,i] * wl_cen[i]**2 * Angstrom / c else: phot[i,:] = phot[i,:] * wl_cen[i]**2 * Angstrom / c elif (units[i] == 'erg/s/A') and \ (target_units == 'AB mag'): if filter_last: L_nu = phot[:,i] * wl_cen[i]**2 * Angstrom / c F_nu = L_nu / (4.0*np.pi*(10.0*pc)**2) phot[:,i] = -2.5*np.log10(F_nu) - 48.6 else: L_nu = phot[i,:] * wl_cen[i]**2 * Angstrom / c F_nu = L_nu / (4.0*np.pi*(10.0*pc)**2) phot[i,:] = -2.5*np.log10(F_nu) - 48.6 elif (units[i] == 'erg/s/A') and \ (target_units == 'Vega mag'): if filter_last: L_nu = phot[:,i] * wl_cen[i]**2 * Angstrom / c F_nu = L_nu / (4.0*np.pi*(10.0*pc)**2) phot[:,i] = -2.5*np.log10(F_nu) - 48.6 - vega_mag[i] else: L_nu = phot[i,:] * wl_cen[i]**2 * Angstrom / c F_nu = L_nu / (4.0*np.pi*(10.0*pc)**2) phot[i,:] = -2.5*np.log10(F_nu) - 48.6 - vega_mag[i] elif (units[i] == 'erg/s/A') and \ (target_units == 'ST mag'): if filter_last: F_lambda = phot[:,i] / (4.0*np.pi*(10.0*pc)**2) phot[:,i] = -2.5*np.log10(F_lambda) - 21.1 else: F_lambda = phot[i,:] / (4.0*np.pi*(10.0*pc)**2) phot[i,:] = -2.5*np.log10(F_lambda) - 21.1 # L_nu to other units elif (units[i] == 'erg/s/Hz') and \ (target_units == 'erg/s/A'): if filter_last: phot[:,i] = phot[:,i] * c / \ (wl_cen[i]**2 * Angstrom) else: phot[i,:] = phot[i,:] * c / \ (wl_cen[i]**2 * Angstrom) elif (units[i] == 'erg/s/Hz') and \ (target_units == 'AB mag'): if filter_last: F_nu = phot[:,i] / (4.0*np.pi*(10.0*pc)**2) phot[:,i] = -2.5*np.log10(F_nu) - 48.6 else: F_nu = phot[i,:] / (4.0*np.pi*(10.0*pc)**2) phot[i,:] = -2.5*np.log10(F_nu) - 48.6 elif (units[i] == 'erg/s/Hz') and \ (target_units == 'Vega mag'): if filter_last: F_nu = phot[:,i] / (4.0*np.pi*(10.0*pc)**2) phot[:,i] = -2.5*np.log10(F_nu) - 48.6 - vega_mag[i] else: F_nu = phot[i,:] / (4.0*np.pi*(10.0*pc)**2) phot[i,:] = -2.5*np.log10(F_nu) - 48.6 - vega_mag[i] elif (units[i] == 'erg/s/Hz') and \ (target_units == 'ST mag'): if filter_last: L_lambda = phot[:,i] * c / (wl_cen[i]**2 * Angstrom) F_lambda = L_lambda / (4.0*np.pi*(10.0*pc)**2) phot[:,i] = -2.5*np.log10(F_lambda) - 21.1 else: L_lambda = phot[i,:] * c / (wl_cen[i]**2 * Angstrom) F_lambda = L_lambda / (4.0*np.pi*(10.0*pc)**2) phot[i,:] = -2.5*np.log10(F_lambda) - 21.1 # AB mag to other units elif (units[i] == 'AB mag') and \ (target_units == 'erg/s/A'): if filter_last: F_nu = 10.0**(-(phot[:,i] + 48.6)/2.5) F_lambda = F_nu * c / (wl_cen[i]**2 * Angstrom) phot[:,i] = F_lambda * 4.0*np.pi*(10.0*pc)**2 else: F_nu = 10.0**(-(phot[i,:] + 48.6)/2.5) F_lambda = F_nu * c / (wl_cen[i]**2 * Angstrom) phot[i,:] = F_lambda * 4.0*np.pi*(10.0*pc)**2 elif (units[i] == 'AB mag') and \ (target_units == 'erg/s/Hz'): if filter_last: F_nu = 10.0**(-(phot[:,i] + 48.6)/2.5) phot[:,i] = F_nu * 4.0*np.pi*(10.0*pc)**2 else: F_nu = 10.0**(-(phot[i,:] + 48.6)/2.5) phot[i,:] = F_nu * 4.0*np.pi*(10.0*pc)**2 elif (units[i] == 'AB mag') and \ (target_units == 'ST mag'): if filter_last: F_nu = 10.0**(-(phot[:,i] + 48.6)/2.5) F_lambda = F_nu * c / (wl_cen[i]**2 * Angstrom) phot[:,i] = -2.5*np.log10(F_lambda) - 21.1 else: F_nu = 10.0**(-(phot[i,:] + 48.6)/2.5) F_lambda = F_nu * c / (wl_cen[i]**2 * Angstrom) phot[i,:] = -2.5*np.log10(F_lambda) - 21.1 elif (units[i] == 'AB mag') and \ (target_units == 'Vega mag'): if filter_last: phot[:,i] = phot[:,i] - vega_mag[i] else: phot[i,:] = phot[i,:] - vega_mag[i] # ST mag to other units elif (units[i] == 'ST mag') and \ (target_units == 'erg/s/A'): if filter_last: F_lambda = 10.0**(-(phot[:,i] + 21.1)/2.5) phot[:,i] = F_lambda * 4.0*np.pi*(10.0*pc)**2 else: F_lambda = 10.0**(-(phot[i,:] + 21.1)/2.5) phot[i,:] = F_lambda * 4.0*np.pi*(10.0*pc)**2 elif (units[i] == 'ST mag') and \ (target_units == 'erg/s/Hz'): if filter_last: F_lambda = 10.0**(-(phot[:,i] + 21.1)/2.5) F_nu = F_lambda * wl_cen[i]**2 * Angstrom / c phot[:,i] = F_nu * 4.0*np.pi*(10.0*pc)**2 else: F_lambda = 10.0**(-(phot[i,:] + 21.1)/2.5) F_nu = F_lambda * wl_cen[i]**2 * Angstrom / c phot[i,:] = F_nu * 4.0*np.pi*(10.0*pc)**2 elif (units[i] == 'ST mag') and \ (target_units == 'AB mag'): if filter_last: F_lambda = 10.0**(-(phot[:,i] + 21.1)/2.5) F_nu = F_lambda * wl_cen[i]**2 * Angstrom / c phot[:,i] = -2.5*np.log10(F_nu) - 48.6 else: F_lambda = 10.0**(-(phot[i,:] + 21.1)/2.5) F_nu = F_lambda * wl_cen[i]**2 * Angstrom / c phot[i,:] = -2.5*np.log10(F_nu) - 48.6 elif (units[i] == 'ST mag') and \ (target_units == 'Vega mag'): if filter_last: F_lambda = 10.0**(-(phot[:,i] + 21.1)/2.5) F_nu = F_lambda * wl_cen[i]**2 * Angstrom / c phot[:,i] = -2.5*np.log10(F_nu) - 48.6 - vega_mag[i] else: F_lambda = 10.0**(-(phot[i,:] + 21.1)/2.5) F_nu = F_lambda * wl_cen[i]**2 * Angstrom / c phot[i,:] = -2.5*np.log10(F_nu) - 48.6 - vega_mag[i] except NameError: raise ValueError("Requested photometric " + "conversion requires " + "filter wavelengths") # Store new units units[i] = target_units