-
-
Notifications
You must be signed in to change notification settings - Fork 421
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.
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
Gamma ray spectra dep #2793
base: master
Are you sure you want to change the base?
Gamma ray spectra dep #2793
Changes from 20 commits
ddc060c
cc6a1e5
1be314c
82dc7de
cf0258d
05f09f5
6e5f459
8d2179d
0ab1571
6dd1c54
ca8c054
9e8c19a
101c037
8484acd
1d0398e
76115ad
906935e
099f7a7
bc5c0c6
9f83760
8d39eb0
0a5ac59
8fb6f48
cd6e099
112f98a
984b95e
0f9976e
923b98c
00a34a0
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -47,6 +47,7 @@ Atharwa Kharkar <[email protected]> atharwa_24 <[email protected] | |
Atul Kumar <[email protected]> | ||
|
||
Anirban Dutta <[email protected]> | ||
Anirban Dutta <[email protected]> | ||
|
||
Barnabás Barna <[email protected]> | ||
|
||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,7 +1,6 @@ | ||
import numpy as np | ||
from numba import njit | ||
|
||
from tardis.energy_input.gamma_ray_estimators import deposition_estimator_kasen | ||
from tardis.energy_input.gamma_ray_grid import ( | ||
distance_trace, | ||
move_packet, | ||
|
@@ -38,17 +37,16 @@ | |
pair_creation_opacity_type, | ||
electron_number_density_time, | ||
mass_density_time, | ||
inv_volume_time, | ||
iron_group_fraction_per_shell, | ||
inner_velocities, | ||
outer_velocities, | ||
times, | ||
dt_array, | ||
times, | ||
effective_time_array, | ||
energy_bins, | ||
energy_df_rows, | ||
energy_plot_df_rows, | ||
energy_out, | ||
energy_deposited_gamma, | ||
energy_deposited_positron, | ||
packets_info_array, | ||
): | ||
"""Propagates packets through the simulation | ||
|
@@ -103,21 +101,23 @@ | |
escaped_packets = 0 | ||
scattered_packets = 0 | ||
packet_count = len(packets) | ||
# Logging does not work with numba | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Remove the "print statement". |
||
print("Entering gamma ray loop for " + str(packet_count) + " packets") | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I would much prefer using logging, and just doing it right above the function call rather than inside the function. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yep now this is moving from being in-development, we should formalize some of these logs (which has already been done in some of the other major functions) |
||
|
||
deposition_estimator = np.zeros_like(energy_df_rows) | ||
|
||
for i in range(packet_count): | ||
packet = packets[i] | ||
time_index = get_index(packet.time_current, times) | ||
energy_deposited_positron[ | ||
packet.shell, time_index | ||
] += packet.positron_energy | ||
|
||
andrewfullard marked this conversation as resolved.
Show resolved
Hide resolved
|
||
if time_index < 0: | ||
print(packet.time_current, time_index) | ||
raise ValueError("Packet time index less than 0!") | ||
|
||
scattered = False | ||
|
||
initial_energy = packet.energy_cmf | ||
# Not used now. Useful for the deposition estimator. | ||
# initial_energy = packet.energy_cmf | ||
|
||
while packet.status == GXPacketStatus.IN_PROCESS: | ||
# Get delta-time value for this step | ||
|
@@ -214,17 +214,6 @@ | |
|
||
packet = move_packet(packet, distance) | ||
|
||
deposition_estimator[packet.shell, time_index] += ( | ||
(initial_energy * 1000) | ||
* distance | ||
* (packet.energy_cmf / initial_energy) | ||
* deposition_estimator_kasen( | ||
comoving_energy, | ||
mass_density_time[packet.shell, time_index], | ||
iron_group_fraction_per_shell[packet.shell], | ||
) | ||
) | ||
|
||
if distance == distance_time: | ||
time_index += 1 | ||
|
||
|
@@ -246,26 +235,9 @@ | |
|
||
packet, ejecta_energy_gained = process_packet_path(packet) | ||
|
||
# Save packets to dataframe rows | ||
# convert KeV to eV / s / cm^3 | ||
energy_df_rows[packet.shell, time_index] += ( | ||
ejecta_energy_gained * 1000 | ||
) | ||
|
||
energy_plot_df_rows[i] = np.array( | ||
[ | ||
i, | ||
ejecta_energy_gained * 1000 | ||
# * inv_volume_time[packet.shell, time_index] | ||
/ dt, | ||
packet.get_location_r(), | ||
packet.time_current, | ||
packet.shell, | ||
compton_opacity, | ||
photoabsorption_opacity, | ||
pair_creation_opacity, | ||
] | ||
) | ||
energy_deposited_gamma[ | ||
packet.shell, time_index | ||
] += ejecta_energy_gained | ||
|
||
if packet.status == GXPacketStatus.PHOTOABSORPTION: | ||
# Packet destroyed, go to the next packet | ||
|
@@ -279,14 +251,18 @@ | |
|
||
if packet.shell > len(mass_density_time[:, 0]) - 1: | ||
rest_energy = packet.nu_rf * H_CGS_KEV | ||
lum_rf = (packet.energy_rf * 1.6022e-9) / dt | ||
bin_index = get_index(rest_energy, energy_bins) | ||
bin_width = ( | ||
energy_bins[bin_index + 1] - energy_bins[bin_index] | ||
) | ||
energy_out[bin_index, time_index] += rest_energy / ( | ||
bin_width * dt | ||
freq_bin_width = bin_width / H_CGS_KEV | ||
energy_out[bin_index, time_index] += ( | ||
packet.energy_rf | ||
/ dt | ||
/ freq_bin_width # Take light crossing time into account | ||
) | ||
|
||
luminosity = packet.energy_rf / dt | ||
packet.status = GXPacketStatus.ESCAPED | ||
escaped_packets += 1 | ||
if scattered: | ||
|
@@ -303,7 +279,7 @@ | |
packet.nu_cmf, | ||
packet.nu_rf, | ||
packet.energy_cmf, | ||
lum_rf, | ||
luminosity, | ||
packet.energy_rf, | ||
packet.shell, | ||
] | ||
|
@@ -313,12 +289,10 @@ | |
print("Scattered packets:", scattered_packets) | ||
|
||
return ( | ||
energy_df_rows, | ||
energy_plot_df_rows, | ||
energy_out, | ||
deposition_estimator, | ||
bin_width, | ||
packets_info_array, | ||
energy_deposited_gamma, | ||
energy_deposited_positron, | ||
) | ||
|
||
|
||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -5,6 +5,7 @@ | |
import radioactivedecay as rd | ||
|
||
from tardis.energy_input.util import KEV2ERG | ||
from tardis.model.matter.decay import IsotopicMassFraction | ||
|
||
logger = logging.getLogger(__name__) | ||
logging.basicConfig(level=logging.INFO) | ||
|
@@ -83,10 +84,10 @@ | |
cumulative_decay_df : pd.DataFrame | ||
total decays for x g of isotope for time 't' | ||
""" | ||
time_delta = u.Quantity(time_delta, u.s) | ||
time_delta = u.Quantity(time_delta, u.d) | ||
total_decays = {} | ||
for shell, inventory in inventories.items(): | ||
total_decays[shell] = inventory.cumulative_decays(time_delta.value) | ||
total_decays[shell] = inventory.cumulative_decays(time_delta.value, "d") | ||
|
||
flattened_dict = {} | ||
|
||
|
@@ -167,3 +168,91 @@ | |
) | ||
|
||
return isotope_decay_df | ||
|
||
|
||
def time_evolve_mass_fraction(raw_isotope_mass_fraction, time_array): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Maybe combine this with the time_evolve_cumilative_decay function below, and just have the combined function return both the new mass fractions and cumulative decays. |
||
""" | ||
Function to evolve the mass fraction of isotopes with time. | ||
|
||
Parameters | ||
---------- | ||
raw_isotope_mass_fraction : pd.DataFrame | ||
isotope mass fraction in mass fractions. | ||
time_array : np.array | ||
array of time in days. | ||
|
||
Returns | ||
------- | ||
time_evolved_isotope_mass_fraction : pd.DataFrame | ||
time evolved mass fraction of isotopes. | ||
""" | ||
|
||
initial_isotope_mass_fraction = raw_isotope_mass_fraction | ||
isotope_mass_fraction_list = [] | ||
|
||
for time in time_array: | ||
|
||
decayed_isotope_mass_fraction = IsotopicMassFraction( | ||
initial_isotope_mass_fraction | ||
).decay(time) | ||
isotope_mass_fraction_list.append(decayed_isotope_mass_fraction) | ||
initial_isotope_mass_fraction = decayed_isotope_mass_fraction | ||
|
||
time_evolved_isotope_mass_fraction = pd.concat( | ||
isotope_mass_fraction_list, keys=time_array, names=["time"] | ||
) | ||
|
||
return time_evolved_isotope_mass_fraction | ||
|
||
|
||
def time_evolve_cumulative_decay( | ||
raw_isotope_mass_fraction, shell_masses, gamma_ray_lines, time_array | ||
): | ||
""" | ||
Function to calculate the total decays for each isotope for each shell at each time step. | ||
|
||
Parameters | ||
---------- | ||
raw_isotope_mass_fraction : pd.DataFrame | ||
isotope abundance in mass fractions. | ||
shell_masses : numpy.ndarray | ||
shell masses in units of g | ||
gamma_ray_lines : pd.DataFrame | ||
gamma ray lines from nndc stored as a pandas dataframe. | ||
time_array : numpy.ndarray | ||
array of time steps in days. | ||
|
||
Returns | ||
------- | ||
time_evolve_decay_df : pd.DataFrame | ||
dataframe of isotopes for each shell with their decay mode, number of decays, radiation type, | ||
radiation energy and radiation intensity at each time step. | ||
|
||
""" | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This function duplicates some code from the function above that could be factored out into a different function |
||
|
||
isotope_decay_df_list = [] | ||
initial_isotope_mass_fraction = raw_isotope_mass_fraction | ||
|
||
dt = np.diff(time_array) | ||
# Do not append dataframes to empty lists | ||
for time in dt: | ||
# | ||
isotope_dict = create_isotope_dicts( | ||
initial_isotope_mass_fraction, shell_masses | ||
) | ||
inventories = create_inventories_dict(isotope_dict) | ||
total_decays = calculate_total_decays(inventories, time) | ||
isotope_df_time = create_isotope_decay_df(total_decays, gamma_ray_lines) | ||
isotope_decay_df_list.append(isotope_df_time) | ||
|
||
decayed_isotope_mass_fraction = IsotopicMassFraction( | ||
initial_isotope_mass_fraction | ||
).decay(time) | ||
|
||
initial_isotope_mass_fraction = decayed_isotope_mass_fraction | ||
|
||
time_evolved_decay_df = pd.concat( | ||
isotope_decay_df_list, keys=time_array, names=["time"] | ||
) | ||
|
||
return time_evolved_decay_df |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Change mailmap!