def direct_ionization_cross_section(self, energy: u.erg):
"""
Calculate direct ionization cross-section.
The cross-sections are calculated one of two ways:
- Using the method of [1]_ for hydrogenic and He-like ions
- Using the scaled cross-sections of [2]_ for all other ions
References
----------
.. [1] Fontes, C. J., et al., 1999, Phys. Rev. A., `59 1329 <https://journals.aps.org/pra/abstract/10.1103/PhysRevA.59.1329>`_
.. [2] Dere, K. P., 2007, A&A, `466, 771 <http://adsabs.harvard.edu/abs/2007A%26A...466..771D>`_
"""
is_hydrogenic = (self.atomic_number - self.charge_state == 1) and (self.atomic_number >= 6)
is_he_like = (self.atomic_number - self.charge_state == 2) and (self.atomic_number >= 10)
if is_hydrogenic or is_he_like:
return self._fontes_cross_section(energy)
else:
return self._dere_cross_section(energy)
python类erg()的实例源码
def _dere_cross_section(self, energy: u.erg):
"""
Calculate direct ionization cross-sections according to [1]_.
References
----------
.. [1] Dere, K. P., 2007, A&A, `466, 771 <http://adsabs.harvard.edu/abs/2007A%26A...466..771D>`_
"""
# Cross-sections from diparams file
cross_section_total = np.zeros(energy.shape)
for ip,bt_c,bt_e,bt_cross_section in zip(self._diparams['ip'], self._diparams['bt_c'], self._diparams['bt_e'],
self._diparams['bt_cross_section']):
U = energy/(ip.to(u.erg))
scaled_energy = 1. - np.log(bt_c)/np.log(U - 1. + bt_c)
f_interp = interp1d(bt_e.value, bt_cross_section.value, kind='cubic', fill_value='extrapolate')
scaled_cross_section = f_interp(scaled_energy.value)*bt_cross_section.unit
# Only nonzero at energies above the ionization potential
scaled_cross_section *= (U.value > 1.)
cross_section = scaled_cross_section*(np.log(U) + 1.)/U/(ip**2)
if not hasattr(cross_section_total, 'unit'):
cross_section_total = cross_section_total*cross_section.unit
cross_section_total += cross_section
return cross_section_total
def _fontes_cross_section(self, energy: u.erg):
"""
Calculate direct ionization cross-section according to [1]_.
References
----------
.. [1] Fontes, C. J., et al., 1999, Phys. Rev. A., `59 1329 <https://journals.aps.org/pra/abstract/10.1103/PhysRevA.59.1329>`_
"""
is_hydrogenic = (self.atomic_number - self.charge_state == 1) and (self.atomic_number >= 6)
U = energy/self.ip
A = 1.13
B = 1 if is_hydrogenic else 2
F = 1 if self.atomic_number < 20 else (140 + (self.atomic_number/20)**3.2)/141
if self.atomic_number >= 16:
c, d, C, D = -0.28394, 1.95270, 0.20594, 3.70590
if self.atomic_number > 20:
C += ((self.atomic_number - 20)/50.5)**1.11
else:
c, d, C, D = -0.80414, 2.32431, 0.14424, 3.82652
Qrp = 1./U*(A*np.log(U) + D*(1. - 1./U)**2 + C*U*(1. - 1./U)**4 + (c/U + d/U**2)*(1. - 1./U))
return B*(np.pi*const.a0.cgs**2)*F*Qrp/(self.ip.to(u.Ry).value**2)
def excitation_autoionization_rate(self):
"""
Calculate ionization rate due to excitation autoionization
"""
# Collision constant
c = (const.h.cgs**2)/((2. * np.pi * const.m_e.cgs)**(1.5) * np.sqrt(const.k_B.cgs))
kBTE = u.Quantity([(const.k_B.cgs * self.temperature) / (de.to(u.erg))
for de in self._easplups['delta_energy']])
# Descale upsilon
shape = self._easplups['bt_upsilon'].shape
xs = np.tile(np.linspace(0, 1, shape[1]), shape[0]).reshape(shape)
args = [xs, self._easplups['bt_upsilon'].value, kBTE.value, self._easplups['bt_c'].value,
self._easplups['bt_type']]
upsilon = u.Quantity(list(map(self.burgess_tully_descale, *args)))
# Rate coefficient
rate = c * upsilon * np.exp(-1 / kBTE) / np.sqrt(self.temperature[np.newaxis,:])
return rate.sum(axis=0)
def proton_collision_rate(self):
"""
Calculates the collision rate for de-exciting and exciting collisions for protons
"""
# Create scaled temperature--these are not stored in the file
bt_t = np.vectorize(np.linspace, excluded=[0, 1], otypes='O')(0, 1, [ups.shape[0]
for ups in self._psplups['bt_rate']])
# Get excitation rates directly from scaled data
energy_ratio = np.outer(const.k_B.cgs*self.temperature, 1.0/self._psplups['delta_energy'].to(u.erg))
ex_rate = np.array(list(map(self.burgess_tully_descale, bt_t, self._psplups['bt_rate'], energy_ratio.T,
self._psplups['bt_c'], self._psplups['bt_type'])))
ex_rate = u.Quantity(np.where(ex_rate > 0., ex_rate, 0.), u.cm**3/u.s).T
# Calculation de-excitation rates from excitation rate
omega_upper = 2.*self._elvlc['J'][self._psplups['upper_level'] - 1] + 1.
omega_lower = 2.*self._elvlc['J'][self._psplups['lower_level'] - 1] + 1.
dex_rate = (omega_lower/omega_upper)*ex_rate*np.exp(1./energy_ratio)
return dex_rate, ex_rate
def get_stellar_mass(self, z=0, rnorm=1., cosmo=Planck15):
"""
Get M/L, L, stellar mass, as measured in R-band
Returns:
M/Lr, log10(Lr), log10(stellar_mass)
"""
from sedpy.observate import Filter
import astropy.units as u
import astropy.constants as const
MLr, MLi, MLk = self.get_M2L()
#plt.plot(self.wave, self.spec); print(MLr)
rfilt = Filter('sdss_r0')
norm_spec = self.get_model(in_place=False)*rnorm
#rband_Jy = rfilt.obj_counts(self.wave, norm_spec)/rfilt.ab_zero_counts*3631
#rband_flam = rband_Jy/(3.34e4*rfilt.wave_pivot**2)#*u.erg/u.second/u.cm**2/u.AA
#dL = Planck15.luminosity_distance(z)
Mr = rfilt.ab_mag(self.wave, norm_spec) - cosmo.distmod(z).value
Lr = 10**(-0.4*(Mr-rfilt.solar_ab_mag))
#Lr = (rband_flam*4*np.pi*dL.to(u.cm).value**2)/3.828e+33*rfilt.wave_pivot*(1+z)
stellar_mass = Lr*MLr
return MLr, np.log10(Lr), np.log10(stellar_mass)
def test_creation():
wvln = np.linspace(1000., 4000., 1024)
flux = np.random.uniform(0., 1., wvln.size)
# try with list, array, Quantity input
s = Spectrum(list(wvln), list(flux))
s = Spectrum(wvln, flux)
s = Spectrum(wvln*u.angstrom, flux*u.erg/u.cm**2/u.angstrom)
# ------------------------------------------------------------------------
# Check that creation fails as expected if:
# 1) shapes don't match
with pytest.raises(ValueError):
s = Spectrum(wvln[:-1], flux)
with pytest.raises(ValueError):
s = Spectrum(wvln, flux[:-1])
# 2) object can't be coerced to a Quantity
with pytest.raises(TypeError):
s = Spectrum(wvln, None)
with pytest.raises(TypeError):
s = Spectrum(None, flux)
with pytest.raises(TypeError):
s = Spectrum(None, None)
# 3) wavelength goes negative
wvln2 = wvln.copy()
wvln2[:-10] *= -1.
with pytest.raises(ValueError):
s = Spectrum(wvln2, flux)
def test_integrate():
subslice = slice(100,200)
wvln = np.linspace(1000., 4000., 1024)
flux = np.zeros_like(wvln)
flux[subslice] = 1./np.ptp(wvln[subslice]) # so the integral is 1
s = Spectrum(wvln*u.angstrom, flux*u.erg/u.cm**2/u.angstrom)
# the integration grid is a sub-section of the full wavelength array
wvln_grid = s.wavelength[subslice]
i_flux = s.integrate(wvln_grid)
assert np.allclose(i_flux.value, 1.) # "close" because this is float comparison
def stellar_info(chain, data):
"""
computes stellar masses and SFRs
"""
gal_do, irlum_dict, nh_dict, BBebv_dict, SFRdict = data.dictkey_arrays #call dictionary info
#relevanta parameters form the MCMC chain
tau_mcmc = chain[:,0]
age_mcmc = chain[:,1]
GA = chain[:,6] - 18. #1e18 is the common normalization factor used in parspace.ymodel in order to have comparable NORMfactors
z = data.z
distance = z2Dlum(z)
#constants
solarlum = const.L_sun.to(u.erg/u.second) #3.839e33
solarmass = const.M_sun
Mstar_list=[]
SFR_list=[]
for i in range (len (tau_mcmc)):
N = 10**GA[i]* 4* pi* distance**2 / (solarlum.value)/ (1+z)
gal_do.nearest_par2dict(tau_mcmc[i], 10**age_mcmc[i], 0.)
tau_dct, age_dct, ebvg_dct=gal_do.t, gal_do.a,gal_do.e
SFR_mcmc =SFRdict[tau_dct, age_dct]
# Calculate Mstar. BC03 templates are normalized to M* = 1 M_sun.
# Thanks to Kenneth Duncan, and his python version of BC03, smpy
Mstar = np.log10(N * 1)
#Calculate SFR. output is in [Msun/yr].
SFR = N * SFR_mcmc
SFR_list.append(SFR.value)
Mstar_list.append(Mstar)
return np.array(Mstar_list) , np.array(SFR_list)
def effective_collision_strength(self):
"""
Calculate the effective collision strength or the Maxwellian-averaged collision
strength, typically denoted by upsilon.
Note
----
Need a more efficient way of calculating upsilon for all transitions. Current method is slow ions with
many transitions, e.g. Fe IX and Fe XI
"""
energy_ratio = np.outer(const.k_B.cgs*self.temperature, 1.0/self._scups['delta_energy'].to(u.erg))
upsilon = np.array(list(map(self.burgess_tully_descale, self._scups['bt_t'], self._scups['bt_upsilon'],
energy_ratio.T, self._scups['bt_c'], self._scups['bt_type'])))
upsilon = u.Quantity(np.where(upsilon > 0., upsilon, 0.))
return upsilon.T
def electron_collision_rate(self):
"""
Calculates the collision rate for de-exciting and exciting collisions for electrons
"""
c = (const.h.cgs**2)/((2. * np.pi * const.m_e.cgs)**(1.5) * np.sqrt(const.k_B.cgs))
upsilon = self.effective_collision_strength()
omega_upper = 2.*self._elvlc['J'][self._scups['upper_level'] - 1] + 1.
omega_lower = 2.*self._elvlc['J'][self._scups['lower_level'] - 1] + 1.
dex_rate = c*upsilon/np.sqrt(self.temperature[:, np.newaxis])/omega_upper
energy_ratio = np.outer(1./const.k_B.cgs/self.temperature, self._scups['delta_energy'].to(u.erg))
ex_rate = omega_upper/omega_lower*dex_rate*np.exp(-energy_ratio)
return dex_rate, ex_rate
def __init__(self, filename_cube, filename_white=None, pixelsize=0.2 * u.arcsec, n_fig=1,
flux_units=1E-20 * u.erg / u.s / u.cm ** 2 / u.angstrom, vmin=None, vmax=None, wave_cal='air'):
"""
Parameters
----------
filename_cube: string
Name of the MUSE datacube .fits file
filename_white: string
Name of the MUSE white image .fits file
pixel_size : float or Quantity, optional
Pixel size of the datacube, if float it assumes arcsecs.
Default is 0.2 arcsec
n_fig : int, optional
XXXXXXXX
flux_units : Quantity
XXXXXXXXXX
"""
# init
self.color = False
self.cmap = ""
self.flux_units = flux_units
self.n = n_fig
plt.close(self.n)
self.wave_cal = wave_cal
self.filename = filename_cube
self.filename_white = filename_white
self.load_data()
self.white_data = fits.open(self.filename_white)[1].data
self.hdulist_white = fits.open(self.filename_white)
self.white_data = np.where(self.white_data < 0, 0, self.white_data)
if not vmin:
self.vmin=np.nanpercentile(self.white_data,0.25)
else:
self.vmin = vmin
if not vmax:
self.vmax=np.nanpercentile(self.white_data,98.)
else:
self.vmax = vmax
self.gc2 = aplpy.FITSFigure(self.filename_white, figure=plt.figure(self.n))
self.gc2.show_grayscale(vmin=self.vmin, vmax=self.vmax)
# self.gc = aplpy.FITSFigure(self.filename, slices=[1], figure=plt.figure(20))
self.pixelsize = pixelsize
gc.enable()
# plt.close(20)
print("MuseCube: Ready!")
def __init__(self, filename_cube, filename_white=None, pixelsize=0.2 * u.arcsec, n_fig=1,
flux_units=1E-20 * u.erg / u.s / u.cm ** 2 / u.angstrom, vmin=None, vmax=None, wave_cal='air'):
"""
Parameters
----------
filename_cube: string
Name of the MUSE datacube .fits file
filename_white: string
Name of the MUSE white image .fits file
pixel_size : float or Quantity, optional
Pixel size of the datacube, if float it assumes arcsecs.
Default is 0.2 arcsec
n_fig : int, optional
XXXXXXXX
flux_units : Quantity
XXXXXXXXXX
"""
# init
self.color = False
self.cmap = ""
self.flux_units = flux_units
self.n = n_fig
plt.close(self.n)
self.wave_cal = wave_cal
self.filename = filename_cube
self.filename_white = filename_white
self.load_data()
self.white_data = fits.open(self.filename_white)[1].data
self.hdulist_white = fits.open(self.filename_white)
self.white_data = np.where(self.white_data < 0, 0, self.white_data)
if not vmin:
self.vmin=np.nanpercentile(self.white_data,0.25)
else:
self.vmin = vmin
if not vmax:
self.vmax=np.nanpercentile(self.white_data,98.)
else:
self.vmax = vmax
self.gc2 = aplpy.FITSFigure(self.filename_white, figure=plt.figure(self.n))
self.gc2.show_grayscale(vmin=self.vmin, vmax=self.vmax)
# self.gc = aplpy.FITSFigure(self.filename, slices=[1], figure=plt.figure(20))
self.pixelsize = pixelsize
gc.enable()
# plt.close(20)
print("MuseCube: Ready!")
