## Copyright (C) 2013 APC CNRS Universite Paris Diderot
## <cecile.cavet@apc.univ-paris7.fr>
## <lejeune@apc.univ-paris7.fr>
## <cecile.portello-roucelle@apc.univ-paris7.fr>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, see http://www.gnu.org/licenses/gpl.html
import numpy as np
from scipy.integrate import trapz, quad
from astropy.constants import c
from astropy import cosmology
#from cosmolopy import magnitudes
# Physical constants in cgs system
PARSEC = 3.086e18 # in cm
CELERITY = 2.99792458e10 # in cm/s ##warning: c in A/s in the original code
FLUX_CONV = 1./(4. * np.pi * 100. * PARSEC**2) # in erg/s/A/cm^2
M_AB_0 = -2.5 * np.log10(3631.e-23) # zero point of AB apparent magnitude
COSMO = 'default'
[docs]def set_cosmo(key):
""" Set global cosmology variable. Possible values are 'default' (H0=70, Om0=0.3) or 'WMAP7'.
"""
global COSMO
if key in ['default', 'WMAP7']:
COSMO = key
[docs]def get_cosmo(key=None):
""" Get astropy cosmology object from global variable. See set_cosmo for more details.
"""
if key is None:
key = COSMO
if key=='default':
return cosmology.FlatLambdaCDM(H0=70.0, Om0=0.3)
elif key=='WMAP7':
return cosmology.WMAP7
[docs]def mag_AB_abs(wave, flux, trans):
""" Compute the absolute magnitude in AB system
See Eq. 3.1 of O. Ilbert thesis
:math:`M_{AB} = -2.5 * \log_{10}( \\frac{\int(F(\lambda,z)T(\lambda)d\lambda)} {\int(T(\lambda)(c/\lambda^2)d\lambda)} ) - M_{AB,0}`
:math:`M_{AB} = -2.5 * \log_{10}(F(\\nu) / \\text{area}) - M_{AB,0}`
**Examples**
>>> import numpy as np
>>> print mag_AB_abs(np.array([5000, 5010, 5020]), np.array([0.0064, 0.0063, 0.0065]), np.array([1, 1, 1]))
-15.418265234
"""
## check size of each array
assert wave.shape==flux.shape
assert wave.shape==trans.shape
## replace nan by 0
flux[np.isnan(flux)] = 0
trans[np.isnan(trans)]= 0
arean = trapz(trans * CELERITY * 1.e8 / wave**2, x=None, dx=1.0, axis=-1) # c in A/s
fmel = trapz(flux * trans, x=None, dx=1.0, axis=-1)
mab = -2.5 * np.log10(fmel / arean) - M_AB_0
return mab
[docs]def mag_AB_to_flux(mab):
""" Compute the flux from the AB magnitude
**Examples**
>>> print mag_AB_to_flux(-15.41)
5.29695457906e-14
"""
return 10**(-0.4 * (mab + M_AB_0))
[docs]def flux_to_mag_AB(flux, wave, system='cgs'):
""" Compute the AB magnitude from the flux
**Examples**
>>> print flux_to_mag_AB(5.29695457906e-14)
-15.41
"""
assert wave.shape==flux.shape
if system == 'SI':
flux = 1.e23 * CELERITY * flux / wave**2
return -2.5 * np.log10(flux) - M_AB_0
[docs]def do_redshift(z, wave, flux):
""" Redshift the SED spectra
**Examples**
>>> do_redshift(0.1, np.array([5000, 5010, 5020]), np.array([0.0064, 0.0063, 0.0065]))
(array([ 5500., 5511., 5522.]), array([ 0.00581818, 0.00572727, 0.00590909]))
"""
assert wave.size==flux.size
wave_z = wave * (1. + z)
flux_z = flux / (1. + z)
return wave_z, flux_z
[docs]def apply_extinction_law(flux, k, ebv=0):
""" Apply inter-galactic extinction laws given E(B-V) for an SED
**Examples**
>>> apply_extinction_law(np.array([0.0064, 0.0063, 0.0065]), np.array([ 4.45, 4.43, 4.41]), ebv=1)
array([ 0.00010621, 0.0001065 , 0.00011192])
"""
## The extinction curve are expressed in k(lambda[A]) vs lambda[A] [k(l) vs l]
## flux_attenuated=flux_intrinsic * 10^[-0.4*k(l)*E(B-V)]
## to do : consider ext laws normalization to confirm E(B-V) to give as an input for user
assert k.shape==flux.shape
flux_ext = flux * 10.**(-0.4 * ebv * k)
## replace nan by previous flux value
flux_ext[np.isnan(flux_ext)] = flux[np.isnan(flux_ext)]
return flux_ext
[docs]def compute_int_mag_AB(wave_res, magc_res, filt_res) :
flux_res = 10**(magc_res/-2.5) *3631 ## flux en Jy
fluxn = trapz(flux_res* filt_res* c.value / (wave_res*10**(-10))**2, x=None, dx=1.0, axis=-1)
arean = trapz(filt_res * c.value / (wave_res*10**(-10))**2, x=None, dx=1.0, axis=-1) ## c en m/s ## wave en A
mag_int = - 2.5*np.log10(fluxn*1./arean ) + 8.9 # compute AB mag from flux in Jy
# no error on photometry accounted for here.
return mag_int
[docs]def compute_pivot_wave(wave_res, filt_res) :
totconv = trapz(wave_res* filt_res, x=None, dx=1.0, axis=-1)
sensi = trapz(filt_res*1./wave_res, x=None, dx=1.0, axis=-1)
pivot_wave = np.sqrt(totconv*1./sensi)
return pivot_wave