Skip to content

21 cm part updated #128

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Anikbh11 wants to merge 3 commits into
base: main
Choose a base branch
from
Open

21 cm part updated #128

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5b04665
21 cm part updated
Anikbh11 Mar 27, 2020
c75c6e3
Added a new option to select line deposition method
Anikbh11 Apr 9, 2020
4edad87
Minor update ensuring lightray redshift is global for wavelength space
Anikbh11 Apr 9, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
120 changes: 114 additions & 6 deletions trident/absorption_spectrum/absorption_line.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,13 +18,24 @@
from yt.utilities.physical_constants import \
charge_proton_cgs, \
mass_electron_cgs, \
speed_of_light_cgs
speed_of_light_cgs, \
planck_constant_cgs, \
boltzmann_constant_cgs
from yt.utilities.on_demand_imports import _scipy, NotAModule
from yt.units.yt_array import YTArray, YTQuantity

special = _scipy.special
signal = _scipy.signal
tau_factor = None
_cs = None

def delta(lam1,lambda_bins):
idx = np.digitize(lam1,lambda_bins,right=True)
if idx == len(lambda_bins):
phi = np.zeros_like(lambda_bins)
else:
phi = signal.unit_impulse(lambda_bins.size,idx)
return phi

def voigt_scipy(a, u):
x = np.asarray(u).astype(np.float64)
Expand Down Expand Up @@ -143,10 +154,10 @@ def voigt_old(a, u):
k1 = k1.astype(np.float64).clip(0)
return k1


def tau_profile(lambda_0, f_value, gamma, v_doppler, column_density,
delta_v=None, delta_lambda=None,
lambda_bins=None, n_lambda=12000, dlambda=0.01):
lambda_bins=None, n_lambda=12000, dlambda=0.01,
deposition_method='voigt'):
r"""
Create an optical depth vs. wavelength profile for an
absorption line using a voigt profile.
Expand All @@ -159,7 +170,7 @@ def tau_profile(lambda_0, f_value, gamma, v_doppler, column_density,
f_value : float
absorption line f-value.
gamma : float
absorption line gamma value.
absorption line gamma value
v_doppler : float in cm/s
doppler b-parameter.
column_density : float in cm^-2
Expand All @@ -180,6 +191,13 @@ def tau_profile(lambda_0, f_value, gamma, v_doppler, column_density,
dlambda : float in angstroms
lambda bin width in angstroms if lambda_bins is None.
Default: 0.01.
:deposition_method: 'voigt' or 'delta'
Sets the line profile in which spectra are deposited. If set to
voigt, the resulting line profiles are deposited as voigt profiles
. If set to delta, the line profiles are set to delta. This is
useful for modelling the 21 cm Forest of neutral hydrogen and in cases
where thermal broadening is to be ignored.
Default: voigt

"""
global tau_factor
Expand Down Expand Up @@ -213,11 +231,101 @@ def tau_profile(lambda_0, f_value, gamma, v_doppler, column_density,
# tau_0
tau_X = tau_factor * column_density * f_value / v_doppler
tau0 = tau_X * lambda_0 * 1e-8

Copy link
Collaborator

@chummels chummels Jun 18, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think whitespace changes like these will fail the flake tests. In general, try to leave whitespace alone.

# dimensionless frequency offset in units of doppler freq
x = _cs / v_doppler * (lam1 / lambda_bins - 1.0)
a = gamma / (4.0 * np.pi * nudop) # damping parameter
phi = voigt(a, x) # line profile
if deposition_method == 'voigt':
phi = voigt(a, x) # line profile
else:
phi = delta(lam1,lambda_bins)
tauphi = tau0 * phi # profile scaled with tau0

return (lambda_bins, tauphi)

def tau_profile_21cm(lambda_0, f_value, gamma, temperature, number_density,
Copy link
Collaborator

@brittonsmith brittonsmith Apr 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think we need both a tau_profile_21cm function and a variable deposition method in the regular tau_profile function. I would say choose one or the other, perhaps whichever option has the least amount of duplicate code.

Copy link
Collaborator

@chummels chummels Jun 18, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed. I would say we could just do with tau_profile and have a keyword which triggers the 21 cm deposition to be different, and in turn, this would get triggered by using a line like: H I 21cm or something?

h_now, delta_v=None, delta_lambda=None,
lambda_bins=None, n_lambda=12000, dlambda=0.01,
deposition_method='voigt'):
r"""
Create an optical depth vs. wavelength profile for the
21 cm forest. The optical depth is calculated using eq. 1
in https://arxiv.org/abs/1510.02296.
At the moment in the implementation, we make a very reasonable
assumption that (1+ 1/H(z) dv/dr) ~ 1.

Parameters
----------

lambda_0 : float in angstroms
central wavelength.
f_value : float
absorption line f-value.
gamma : float
absorption line gamma value. For this case, represents the einstein coefficient A_10
temperature : float in Kelvin
Gas temperature. Assumption that T_K = T_S is made here.
number_density : float in cm^-3
neutral hydrogen number density
h_now ; float in s^-1
Hhubble constant at redshift z. H(z)
delta_v : float in cm/s
velocity offset from lambda_0.
Default: None (no shift).
delta_lambda : float in angstroms
wavelength offset.
Default: None (no shift).
lambda_bins : array in angstroms
wavelength array for line deposition. If None, one will be
created using n_lambda and dlambda.
Default: None.
n_lambda : int
size of lambda bins to create if lambda_bins is None.
Default: 12000.
dlambda : float in angstroms
lambda bin width in angstroms if lambda_bins is None.
Default: 0.01.
:deposition_method: 'voigt' or 'delta'
Sets the line profile in which spectra are deposited. If set to
voigt, the resulting line profiles are deposited as voigt profiles
. If set to delta, the line profiles are set to delta. This is
useful for modelling the 21 cm Forest of neutral hydrogen and in cases
where thermal broadening is to be ignored.
Default: voigt

"""
global tau_factor
if tau_factor is None:
lam0 = YTQuantity(lambda_0,'angstrom')
gam = YTQuantity(gamma,'1/s')
tau_factor = ((3 / 32 / np.pi) * planck_constant_cgs * gam *
speed_of_light_cgs * lam0**2 / boltzmann_constant_cgs
).in_cgs().d

global _cs
if _cs is None:
_cs = speed_of_light_cgs.d[()]

# shift lambda_0 by delta_v
if delta_v is not None:
lam1 = lambda_0 * (1 + delta_v / _cs)
elif delta_lambda is not None:
lam1 = lambda_0 + delta_lambda
else:
lam1 = lambda_0

# create wavelength
if lambda_bins is None:
lambda_bins = lam1 + \
np.arange(n_lambda, dtype=np.float) * dlambda - \
n_lambda * dlambda / 2 # wavelength vector (angstroms)

# tau_0
tau0 = (tau_factor * number_density) / (temperature * h_now)
if deposition_method == 'voigt':
raise RuntimeError('21 cm currently only supports delta profile deposition')
else:
phi = delta(lam1,lambda_bins) #line profile as a delta
tauphi = tau0 * phi # profile scaled with tau0

return (lambda_bins, tauphi)
Expand Down
77 changes: 59 additions & 18 deletions trident/absorption_spectrum/absorption_spectrum.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,9 @@

from trident.absorption_spectrum.absorption_line import \
tau_profile
from trident.absorption_spectrum.absorption_line import \
tau_profile_21cm
from yt.utilities.cosmology import Cosmology

pyfits = _astropy.pyfits

Expand Down Expand Up @@ -88,16 +91,27 @@ class AbsorptionSpectrum(object):
wavelength. If set to velocity, the spectra are flux vs.
velocity offset from the rest wavelength of the absorption line.
Default: wavelength

:deposition_method: 'voigt' or 'delta'

Sets the line profile in which spectra are deposited. If set to
voigt, the resulting line profiles are deposited as voigt profiles
. If set to delta, the line profiles are set to delta. This is
useful for modelling the 21 cm Forest of neutral hydrogen and in cases
where thermal broadening is to be ignored.
Default: voigt

"""

def __init__(self, lambda_min, lambda_max, n_lambda=None, dlambda=None,
bin_space='wavelength'):
bin_space='wavelength',deposition_method='voigt'):
Copy link
Collaborator

@chummels chummels Jun 18, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minor issue, but separate keyword arguments by a space.


if bin_space not in _bin_space_units:
raise RuntimeError(
'Invalid bin_space value: "%s". Valid values are: "%s".' %
(bin_space, '", "'.join(list(_bin_space_units))))
self.bin_space = bin_space
self.deposition_method = deposition_method
lunits = _bin_space_units[self.bin_space]

if dlambda is not None:
Expand Down Expand Up @@ -468,8 +482,10 @@ def make_spectrum(self, input_object, output_file=None,
input_ds.domain_left_edge = input_ds.domain_left_edge.to('code_length')
input_ds.domain_right_edge = input_ds.domain_right_edge.to('code_length')

if self.bin_space == 'velocity':
Copy link
Collaborator

@chummels chummels Jun 18, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So we won't be able to do in velocity space? Seems like leaving this in would allow for it, but maybe not?

self.zero_redshift = getattr(input_ds, 'current_redshift', 0)
self.hubble_constant = getattr(input_ds,'hubble_constant')
Copy link
Collaborator

@chummels chummels Jun 18, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I fear that this may break functionality for simulations which aren't cosmological. Notice for the redshift, there is a final argument that sets it to 0 in the case that there is no current_redshift parameter. This is to ensure that the code works on non-cosmological datasets by just setting the current redshift value to 0.

Could you group these parameter instantiations along with your co = Cosmology and hnow code below, and maybe put it inside an if/then block for deposition_method='delta', or a try/except block and give a default value for the z=0 case, so non-cosmological datasets will still run?

self.omega_matter = getattr(input_ds,'omega_matter')
self.omega_lambda = getattr(input_ds,'omega_lambda')
self.zero_redshift = getattr(input_ds, 'current_redshift', 0)

# temperature field required to calculate voigt profile widths
if ('temperature' not in input_ds.derived_field_list) and \
Expand Down Expand Up @@ -707,6 +723,7 @@ def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity,
for store, line in parallel_objects(self.line_list, njobs=njobs,
storage=self.line_observables_dict):
column_density = field_data[line['field_name']] * field_data['dl']
number_density = field_data[line['field_name']]
if (column_density < 0).any():
mylog.warning(
"Setting negative densities for field %s to 0! Bad!" % line['field_name'])
Expand Down Expand Up @@ -738,19 +755,23 @@ def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity,

# the total number of absorbers per transition
n_absorbers = len(lambda_obs)
temperature = field_data['temperature'].d
Copy link
Collaborator

@chummels chummels Jun 18, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this necessary to strip off the Kelvin unit?


# thermal broadening b parameter
thermal_b = np.sqrt((2 * boltzmann_constant_cgs *
field_data['temperature']) /
line['atomic_mass'])

field_data['temperature']) /
line['atomic_mass'])
# the actual thermal width of the lines
thermal_width = (lambda_obs * thermal_b /
c_kms).to('angstrom')


co = Cosmology(self.hubble_constant,self.omega_matter,self.omega_lambda,0.0)
h_now = co.hubble_parameter(self.zero_redshift).d #H(z) s^-1
# Sanitize units for faster runtime of the tau_profile machinery.
lambda_0 = line['wavelength'].d # line's rest frame; angstroms
cdens = column_density.in_units("cm**-2").d # cm**-2
ndens = number_density.in_units("cm**-3").d # cm**-2
thermb = thermal_b.to('cm/s').d # thermal b coefficient; cm / s
dlambda = delta_lambda.d # lambda offset; angstroms
# Array to store sum of the tau values for each index in the
Expand All @@ -777,7 +798,6 @@ def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity,
# 3) a bin width will be divisible by vbin_width times a power of
# 10; this will assure we don't get spikes in the deposited
# spectra from uneven numbers of vbins per bin

if self.bin_space == 'wavelength':
my_width = thermal_width
elif self.bin_space == 'velocity':
Expand All @@ -786,10 +806,9 @@ def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity,
raise RuntimeError('What bin space is this?')

resolution = my_width / self.bin_width
n_vbins_per_bin = (10 ** (np.ceil( np.log10( subgrid_resolution /
resolution) ).clip(0, np.inf) ) ).astype('int')
n_vbins_per_bin = (10 ** (np.ceil( np.log10(subgrid_resolution/
resolution)).clip(0, np.inf))).astype('int')
vbin_width = self.bin_width.d / n_vbins_per_bin

# a note to the user about which lines components are unresolved
if (my_width < self.bin_width).any():
mylog.info("%d out of %d line components will be " +
Expand Down Expand Up @@ -851,7 +870,6 @@ def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity,
# may be <0 or >the size of the array. In the end, we deposit
# the bins that actually overlap with the AbsorptionSpectrum's
# range in lambda.

left_index, center_index, right_index = \
self._get_bin_indices(
my_lambda, self.bin_width,
Expand All @@ -872,6 +890,7 @@ def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity,
my_vbins = vbins * \
wavelength_zero_point.d / c_kms.d + \
wavelength_zero_point.d

else:
raise RuntimeError('What bin_space is this?')

Expand All @@ -882,11 +901,20 @@ def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity,
'increasing the bin size.') % my_vbins.size)

# the virtual bins and their corresponding opacities
my_vbins, vtau = \
tau_profile(
lambda_0, line['f_value'], line['gamma'],
thermb[i], cdens[i],
delta_lambda=dlambda[i], lambda_bins=my_vbins)
if (line['label'] == '21 cm'):
Copy link
Collaborator

@brittonsmith brittonsmith Apr 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Similar to my comment above, it's somewhat redundant to be able to switch the deposition method and change the behavior for one specific line. Because of the huge range of wavelengths between the hyperfine transition and the electron energy level transitions, I don't think anyone should be making a spectrum that includes both the 21 cm and something from the Lyman series. I'd rather see a workflow where the user specifies that they want only the 21 cm line and the delta deposition method manually. This would remove the need to have separate function calls here, which is going to be more difficult to maintain and sets a precedent for adding even more.

Copy link
Collaborator

@chummels chummels Jun 18, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agree with Britton.

my_vbins, vtau = \
tau_profile_21cm(
lambda_0, line['f_value'], line['gamma'],
temperature[i], ndens[i], h_now,
delta_lambda=dlambda[i], lambda_bins=my_vbins,
deposition_method=self.deposition_method)
else:
my_vbins, vtau = \
tau_profile(
lambda_0, line['f_value'], line['gamma'],
thermb[i], cdens[i],
delta_lambda=dlambda[i], lambda_bins=my_vbins,
deposition_method=self.deposition_method)

# If tau has not dropped below min tau threshold by the
# edges (ie the wings), then widen the wavelength
Expand Down Expand Up @@ -925,13 +953,26 @@ def _add_lines_to_spectrum(self, field_data, use_peculiar_velocity,
# normal use of the word by observers. It is an equivalent
# with in tau, not in flux, and is only used internally in
# this subgrid deposition as EW_tau.
vEW_tau = vtau * vbin_width[i]
# The different behavior for the 21 cm line here stems from
# the fact that line profiles are treated as deltas instead
# of voigt profiles.
if self.deposition_method == 'voigt':
vEW_tau = vtau * vbin_width[i]
elif self.deposition_method == 'delta':
vEW_tau = vtau
else:
raise RuntimeError('Unknown line deposition method')
EW_tau = np.zeros(right_index - left_index)
EW_tau_indices = np.arange(left_index, right_index)
for k, val in enumerate(EW_tau_indices):
EW_tau[k] = vEW_tau[n_vbins_per_bin[i] * k:
n_vbins_per_bin[i] * (k + 1)].sum()
EW_tau = EW_tau/self.bin_width.d
if self.deposition_method == 'voigt':
EW_tau = EW_tau/self.bin_width.d
elif self.deposition_method == 'delta':
EW_tau = EW_tau
else:
raise RuntimeError('Unknown line deposition method')

# only deposit EW_tau bins that actually intersect the original
# spectral wavelength range (i.e. lambda_field)
Expand Down
6 changes: 4 additions & 2 deletions trident/data/line_lists/lines.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#Ion Wavelength [A] gamma f_value alt. name
H I 1215.670000 4.690000e+08 4.160000e-01 Ly a
H I 2.1000000e9 2.850000e-15 4.160000e-01 21 cm
Copy link
Collaborator

@brittonsmith brittonsmith Apr 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is quite pedantic, but the rest of the H lines are in order of decreasing wavelength. I would put this line above Lya.

H I 1025.722200 5.570000e+07 7.910000e-02 Ly b
H I 972.536740 1.280000e+07 2.900000e-02 Ly c
H I 949.742980 4.120000e+06 1.390000e-02 Ly d
Expand Down Expand Up @@ -49,7 +50,7 @@ C I 1193.031000 6.390000e+07 4.090000e-02
C I 1188.833000 1.950000e+07 1.240000e-02
C I 1157.910000 3.520000e+07 2.120000e-02
C II 1335.663000 4.760000e+07 1.270000e-02 C II* 1336
C II 1334.532000 2.410000e+08 1.290000e-01
C II 1334.532000 2.400000e+08 1.280000e-01
Copy link
Collaborator

@brittonsmith brittonsmith Apr 21, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@chummels, can you comment on these changes to the line data?

Copy link
Collaborator

@chummels chummels Jun 18, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These changes are in line with the Morton data, but it just seems odd to include these minor changes in a PR that's about depositing HI as 21 cm lines. But they can stay, as I don't think they'll break the answer tests. I don't think so anyway.

C II 1037.018000 1.460000e+09 1.180000e-01 C II* 1037
C II 1036.337000 7.380000e+08 1.190000e-01
C II 903.962000 2.700000e+09 3.310000e-01
Expand Down Expand Up @@ -79,7 +80,7 @@ N V 1242.804000 3.370000e+08 7.800000e-02
N V 1238.821000 3.400000e+08 1.560000e-01
O I 1306.029000 6.760000e+07 5.190000e-02 O I* 1306
O I 1304.858000 2.030000e+08 5.180000e-02 O I* 1305
O I 1302.168000 3.410000e+08 5.200000e-02
O I 1302.168000 3.150000e+08 4.800000e-02
O I 1039.230000 9.430000e+07 9.160000e-03
O I 988.773000 2.260000e+08 4.640000e-02
O I 988.655000 5.770000e+07 8.460000e-03
Expand Down Expand Up @@ -196,6 +197,7 @@ Ar VII 585.748000 7.830000e+09 1.210000e+00
Ca X 574.010000 3.200000e+09 1.600000e-01
Ca X 557.765000 3.500000e+09 3.260000e-01
Fe II 1608.450830 1.910000e+08 5.910000e-02
Fe II 2382.765200 3.130000e+08 0.320000e-01
Fe II 1143.225730 1.000000e+08 1.900000e-02
Fe II 1127.098400 6.000000e+06 1.100000e-03
Fe II 1125.447630 1.000000e+08 1.600000e-02
Expand Down
Loading