Skip to content
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 refactor #2139

Merged
merged 7 commits into from
Jul 7, 2023
Merged
Changes from all commits
Commits
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
142 changes: 142 additions & 0 deletions tardis/energy_input/GXPacket.py
Original file line number Diff line number Diff line change
@@ -3,6 +3,14 @@
from numba import int64, float64
from numba.experimental import jitclass

from tardis.energy_input.samplers import sample_decay_time, sample_energy
from tardis.energy_input.util import (
get_index,
get_random_unit_vector,
doppler_factor_3d,
H_CGS_KEV,
)


class GXPacketStatus(IntEnum):
BETA_DECAY = -1
@@ -68,3 +76,137 @@ def get_location_r(self):
+ self.location[1] ** 2.0
+ self.location[2] ** 2.0
)


# @njit(**njit_dict_no_parallel)
def initialize_packet_properties(
isotope_energy,
isotope_intensity,
positronium_energy,
positronium_intensity,
positronium_fraction,
packet_energy,
k,
tau_start,
tau_end,
initial_radius,
times,
effective_times,
inventory,
average_power_per_mass,
uniform_packet_energies=True,
):
"""Initialize the properties of an individual packet

Parameters
----------
isotope_energy : numpy array
_description_
isotope_intensity : numpy array
_description_
positronium_energy : numpy array
_description_
positronium_intensity : numpy array
_description_
positronium_fraction : float
_description_
packet_energy : float
_description_
k : int
_description_
tau_start : float
_description_
tau_end : float
_description_
initial_radius : float
_description_
times : numpy array
_description_
effective_times : numpy array
_description_

Returns
-------
_type_
_description_
"""
decay_time = np.inf

# packets sampled uniformly in time
if not uniform_packet_energies:
z = np.random.random()
decay_time = z * times[0] + (1 - z) * times[-1]
else:
decay_time = sample_decay_time(
tau_start,
end_tau=tau_end,
decay_time_min=0,
decay_time_max=times[-1],
)

decay_time_index = get_index(decay_time, times)

cmf_energy = sample_energy(isotope_energy, isotope_intensity)

z = np.random.random()
if z < positronium_fraction:
z = np.random.random()
if cmf_energy == 511 and z > 0.25:
cmf_energy = sample_energy(
positronium_energy, positronium_intensity
)

energy_factor = 1.0
if decay_time < times[0]:
energy_factor = decay_time / times[0]
decay_time = times[0]

# compute the new energy factor based on the ratio of the end-chain decay power
# to the average decay power
if not uniform_packet_energies:
energy_factor = (
get_chain_decay_power_per_ejectamass(
inventory,
decay_time,
times[0],
isotope_energy * isotope_intensity / 100,
positronium_energy * positronium_intensity / 100,
tau_end,
)
/ average_power_per_mass
)

scaled_r = initial_radius * effective_times[decay_time_index]

initial_location = scaled_r * get_random_unit_vector()
initial_direction = get_random_unit_vector()
initial_energy = packet_energy * energy_factor

# draw a random gamma-ray in shell
packet = GXPacket(
initial_location,
initial_direction,
1.0,
initial_energy,
0.0,
0.0,
GXPacketStatus.IN_PROCESS,
k,
decay_time,
)

packet.energy_rf = packet.energy_cmf / doppler_factor_3d(
packet.direction,
packet.location,
packet.time_current,
)

packet.nu_cmf = cmf_energy / H_CGS_KEV

packet.nu_rf = packet.nu_cmf / doppler_factor_3d(
packet.direction,
packet.location,
packet.time_current,
)

return packet, decay_time_index
Loading