def norm_apec(self, z):
"""
The normalization factor of the XSPEC APEC model assuming
EM = 1 (i.e., n_e = n_H = 1 cm^-3, and V = 1 cm^3)
norm = 1e-14 / (4*pi* (D_A * (1+z))^2) * int(n_e * n_H) dV
unit: [cm^-5]
This value will be used to calculate the cooling function values.
References
----------
* XSPEC: APEC model
https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XSmodelApec.html
"""
da = self.angular_diameter_distance(z, unit="cm")
norm = 1e-14 / (4*math.pi * (da * (1+z))**2)
return norm
python类cm()的实例源码
def ionization_rate(self):
"""
Total ionization rate.
Includes contributions from both direct ionization and excitation-autoionization
See Also
--------
direct_ionization_rate
excitation_autoionization_rate
"""
di_rate = self.direct_ionization_rate()
di_rate = np.zeros(self.temperature.shape)*u.cm**3/u.s if di_rate is None else di_rate
ea_rate = self.excitation_autoionization_rate()
ea_rate = np.zeros(self.temperature.shape)*u.cm**3/u.s if ea_rate is None else ea_rate
return di_rate + ea_rate
def recombination_rate(self):
"""
Total recombination rate.
Includes contributions from dielectronic recombination and radiative recombination.
See Also
--------
radiative_recombination_rate
dielectronic_recombination_rate
"""
rr_rate = self.radiative_recombination_rate()
rr_rate = np.zeros(self.temperature.shape)*u.cm**3/u.s if rr_rate is None else rr_rate
dr_rate = self.dielectronic_recombination_rate()
dr_rate = np.zeros(self.temperature.shape)*u.cm**3/u.s if dr_rate is None else dr_rate
return rr_rate + dr_rate
def GetFluxRatio(sptlist, Tsec, xgrid):
"""
Returns the flux ratio between the secondary star of temperature Tsec
and the (possibly multiple) primary star(s) given in the
'sptlist' list (given as spectral types)
xgrid is a np.ndarray containing the x-coordinates to find the
flux ratio at (in nm)
"""
prim_flux = np.zeros(xgrid.size)
# Determine the flux from the primary star(s)
for spt in sptlist:
end = search("[0-9]", spt).end()
T = MS.Interpolate(MS.Temperature, spt[:end])
R = MS.Interpolate(MS.Radius, spt[:end])
prim_flux += Planck(xgrid * units.nm.to(units.cm), T) * R ** 2
# Determine the secondary star flux
s_spt = MS.GetSpectralType(MS.Temperature, Tsec)
R = MS.Interpolate(MS.Radius, s_spt)
sec_flux = Planck(xgrid * units.nm.to(units.cm), Tsec) * R ** 2
return sec_flux / prim_flux
def _filter_streamlines(self, streamline, close_threshold=0.05,
loop_length_range: u.cm =[2.e+9, 5.e+10]*u.cm, **kwargs):
"""
Check extracted loop to make sure it fits given criteria. Return True if it passes.
Parameters
----------
streamline : yt streamline object
close_threshold : `float`
percentage of domain width allowed between loop endpoints
loop_length_range : `~astropy.Quantity`
minimum and maximum allowed loop lengths (in centimeters)
"""
streamline = streamline[np.all(streamline != 0.0, axis=1)]
loop_length = np.sum(np.linalg.norm(np.diff(streamline, axis=0), axis=1))
if np.fabs(streamline[0, 2] - streamline[-1, 2]) > close_threshold*self.extrapolated_3d_field.domain_width[2]:
return False
elif loop_length > loop_length_range[1].to(u.cm).value or loop_length < loop_length_range[0].to(u.cm).value:
return False
else:
return True
def peek(self, figsize=(10, 10), color='b', alpha=0.75,
print_to_file=None, **kwargs):
"""
Show extracted fieldlines overlaid on HMI image.
"""
fig = plt.figure(figsize=figsize)
ax = fig.gca(projection=self.hmi_map)
self.hmi_map.plot()
ax.set_autoscale_on(False)
for stream, _ in self.streamlines:
ax.plot(self._convert_angle_to_length(stream[:, 0]*u.cm,
working_units=u.arcsec).to(u.deg),
self._convert_angle_to_length(stream[:, 1]*u.cm,
working_units=u.arcsec).to(u.deg),
alpha=alpha, color=color, transform=ax.get_transform('world'))
if print_to_file is not None:
plt.savefig(print_to_file, **kwargs)
plt.show()
def flatten_emissivities(channel, emission_model):
"""
Compute product between wavelength response and emissivity for all ions
"""
flattened_emissivities = []
for ion in emission_model:
wavelength, emissivity = emission_model.get_emissivity(ion)
if wavelength is None or emissivity is None:
flattened_emissivities.append(None)
continue
interpolated_response = splev(wavelength.value, channel['wavelength_response_spline'], ext=1)
em_summed = np.dot(emissivity.value, interpolated_response)
unit = emissivity.unit*u.count/u.photon*u.steradian/u.pixel*u.cm**2
flattened_emissivities.append(u.Quantity(em_summed, unit))
return flattened_emissivities
def __init__(self, *args, **kwargs):
if len(args) == 1 and os.path.exists(args[0]):
data, header, wavelength = self._restore_from_file(args[0], **kwargs)
elif all([k in kwargs for k in ['data', 'header', 'wavelength']]):
data = kwargs.get('data')
header = kwargs.get('header')
wavelength = kwargs.get('wavelength')
else:
raise ValueError('''EISCube can only be initialized with a valid FITS file or NumPy
array with an associated wavelength and header.''')
# check dimensions
if data.shape[-1] != wavelength.shape[0]:
raise ValueError('''Third dimension of data cube must have the same length as
wavelength.''')
self.meta = header.copy()
self.wavelength = wavelength
self.data = data
self.cmap = kwargs.get('cmap', sunpy.cm.get_cmap('hinodexrt'))
self._fix_header()
def load_results(self, loop):
"""
Load EBTEL output for a given loop object.
Parameters
----------
loop : `synthesizAR.Loop` object
"""
# load text
N_s = len(loop.field_aligned_coordinate)
_tmp = np.loadtxt(loop.hydro_configuration['output_filename'])
# reshape into a 1D loop structure with units
time = _tmp[:, 0]*u.s
electron_temperature = np.outer(_tmp[:, 1], np.ones(N_s))*u.K
ion_temperature = np.outer(_tmp[:, 2], np.ones(N_s))*u.K
density = np.outer(_tmp[:, 3], np.ones(N_s))*(u.cm**(-3))
velocity = np.outer(_tmp[:, -2], np.ones(N_s))*u.cm/u.s
# flip sign of velocity at apex
i_mirror = np.where(np.diff(loop.coordinates.value[:, 2]) > 0)[0][-1] + 2
velocity[:, i_mirror:] = -velocity[:, i_mirror:]
return time, electron_temperature, ion_temperature, density, velocity
def non_equilibrium_ionization(self, time: u.s, temperature: u.K, density: u.cm**(-3), rate_matrix=None,
initial_condition=None):
"""
Compute the ionization fraction in non-equilibrium for a given temperature and density
timeseries.
"""
if rate_matrix is None:
rate_matrix = self._rate_matrix()
if initial_condition is None:
initial_condition = self.equilibrium_ionization(rate_matrix=rate_matrix)
interpolate_indices = [np.abs(self.temperature - t).argmin() for t in temperature]
y = np.zeros(time.shape + (self.atomic_number+1,))
y[0, :] = initial_condition[interpolate_indices[0], :]
identity = u.Quantity(np.eye(self.atomic_number + 1))
for i in range(1, time.shape[0]):
dt = time[i] - time[i-1]
term1 = identity - density[i]*dt/2.*rate_matrix[interpolate_indices[i], :, :]
term2 = identity + density[i-1]*dt/2.*rate_matrix[interpolate_indices[i-1], :, :]
y[i, :] = np.linalg.inv(term1) @ term2 @ y[i-1, :]
y[i, :] = np.fabs(y[i, :])
y[i, :] /= y[i, :].sum()
return u.Quantity(y)
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 emissivity(self, density: u.cm**(-3), include_energy=False, **kwargs):
"""
Calculate emissivity for all lines as a function of temperature and density
"""
populations = self.level_populations(density, include_protons=kwargs.get('include_protons', True))
if populations is None:
return (None, None)
wavelengths = np.fabs(self._wgfa['wavelength'])
# Exclude 0 wavelengths which correspond to two-photon decays
upper_levels = self._wgfa['upper_level'][wavelengths != 0*u.angstrom]
a_values = self._wgfa['A'][wavelengths != 0*u.angstrom]
wavelengths = wavelengths[wavelengths != 0*u.angstrom]
if include_energy:
energy = const.h.cgs*const.c.cgs/wavelengths.to(u.cm)
else:
energy = 1.*u.photon
emissivity = populations[:, :, upper_levels - 1]*(a_values*energy)
return wavelengths, emissivity
def cm_per_pix(self, z):
"""
Calculate the transversal length (unit: cm) corresponding to
1 ACIS pixel (i.e., 0.492 arcsec) at the *angular diameter distance*
of z.
"""
return self.kpc_per_pix(z) * au.kpc.to(au.cm)
def preprocessor(self, table, line, index):
line = line.strip().split()
if index == 0:
filetype = int(line[0])
table.append([filetype])
if filetype == 1:
self.dtypes = [int,float,float,float,float]
self.units = [None,(u.cm**3)/u.s,None,u.K,u.K]
self.headings = ['fit_type', 'A_fit', 'B_fit', 'T0_fit', 'T1_fit']
self.descriptions = ['fit type','A fit parameter','B fit parameter',
'T0 fit parameter','T1 fit parameter']
elif filetype == 2:
self.dtypes = [int,float,float,float,float,float,float]
self.units = [None,(u.cm**3)/u.s,None,u.K,u.K,None,u.K]
self.headings = ['fit_type', 'A_fit', 'B_fit', 'T0_fit', 'T1_fit', 'C_fit', 'T2_fit']
self.descriptions = ['fit type','A fit parameter','B fit parameter',
'T0 fit parameter','T1 fit parameter','C fit parameter',
'T2 fit parameter']
elif filetype == 3:
self.dtypes = [int,float,float]
self.units = [None,(u.cm**3)/u.s,None]
self.headings = ['fit_type', 'A_fit', 'eta_fit']
self.descriptions = ['fit type','A rad fit parameter','eta fit parameter']
else:
raise ValueError('Unrecognized .rrparams filetype {}'.format(filetype))
else:
if table[0][0] == 1 or table[0][0] == 2:
table[0] += line[3:]
else:
table[0] += line[2:]
def preprocessor(self, table, line, index):
line = line.strip().split()
if index == 0:
self._drparams_filetype = int(line[0])
if self._drparams_filetype == 1:
# Badnell type
self.dtypes = [int,float,float]
self.units = [None,u.K,(u.cm**3)/u.s*(u.K**(3/2))]
self.headings = ['fit_type', 'E_fit', 'c_fit']
self.descriptions = ['fit type', 'E fit parameter','c fit parameter']
elif self._drparams_filetype == 2:
# Shull type
self.dtypes = [int,float,float,float,float]
self.units = [None,(u.cm**3)/u.s*(u.K**(3/2)),u.dimensionless_unscaled,u.K,u.K]
self.headings = ['fit_type', 'A_fit', 'B_fit', 'T0_fit', 'T1_fit']
self.descriptions = ['fit type','A fit coefficient','B fit coefficient',
'T0 fit coefficient','T1 fit coefficient']
else:
raise ValueError('Unrecognized drparams filetype {}'.format(self._drparams_filetype))
else:
if self._drparams_filetype == 1:
tmp = np.array(line[2:],dtype=float)
if index%2 == 0:
tmp_col = table[-1]
for i in range(tmp.shape[0]):
table.append([self._drparams_filetype,tmp_col[i],tmp[i]])
del table[0]
else:
table.append(tmp)
else:
table.append([self._drparams_filetype]+line[2:])
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 __init__(self, name, z, x, y):
self.name = name
self.z = z
self.x = x
self.y = y
self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc
self.radToKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc
# setting up final xl file
def __init__(self, name, z, x, y):
self.name = name
self.z = z
self.x = x
self.y = y
self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc
self.radToKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc
# We grab our galaxy data and make a list of galaxy objects
def __init__(self, name, z, x, y):
self.name = name
self.z = z
self.x = x
self.y = y
self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc
self.radToKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc
# Grabbing and filling galaxy data
def __init__(self, name, z, x, y):
self.name = name
self.z = z
self.x = x
self.y = y
self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc
self.radToKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc
# We grab our galaxy data and make a list of galaxy objects
def __init__(self, name, z, x, y):
self.name = name
self.z = z
self.x = x
self.y = y
self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc
self.radToKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc
# setting up final xl file
def __init__(self, name, z, x, y):
self.name = name
self.z = z
self.x = x
self.y = y
self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc
self.radToKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc
# Grabbing and filling galaxy data
def __init__(self, name, z, x, y):
self.name = name
self.z = z
self.x = x
self.y = y
self.lDcm = cosmo.luminosity_distance(self.z)*u.Mpc.to(u.cm) / u.Mpc
self.radToKpc = conv.arcsec_per_kpc_proper(self.z)*0.05/u.arcsec*u.kpc
# Grabbing and filling galaxy data
def test_init(self):
afrho = Afrho(1000 * u.cm)
assert afrho.value == 1000
assert afrho.unit == u.cm
def test_scaler_ops(self):
afrho = Afrho(1000 * u.cm)
afrho = afrho / 2
assert afrho == 500 * u.cm
def test_quantity_ops(self):
afrho = Afrho(1000 * u.cm)
v = afrho * 2 * u.cm
assert v == 2000 * u.cm**2
assert not isinstance(v, Afrho)
def test_from_flam(self):
fluxd = 6.730018324465526e-14 * u.W / u.m**2 / u.um
aper = 1 * u.arcsec
eph = dict(rh=1.5 * u.au, delta=1.0 * u.au)
S = 1869 * u.W / u.m**2 / u.um
afrho = Afrho.from_fluxd(None, fluxd, aper, eph, S=S)
assert np.isclose(afrho.cm, 1000)
def test_from_fnu(self):
fluxd = 6.161081515869728 * u.mJy
nu = 2.998e14 / 11.7 * u.Hz
aper = 1 * u.arcsec
eph = dict(rh=1.5 * u.au, delta=1.0 * u.au)
S = 1.711e14 * u.Jy
afrho = Afrho.from_fluxd(nu, fluxd, aper, eph, S=S)
assert np.isclose(afrho.cm, 1000.0)
def test_to_phase(self):
afrho = Afrho(10 * u.cm).to_phase(15 * u.deg, 0 * u.deg)
assert np.isclose(afrho.cm, 5.8720)
def test_init(self):
efrho = Efrho(1000 * u.cm)
assert efrho.value == 1000
assert efrho.unit == u.cm
