Skip to content

Commit e36a597

Browse files
arguellesclaude
arguelles
and
claude
committed
Add EmittingEarthAtm class for atmospheric neutrino production
This implements PR #41 functionality with improvements based on code review: - EmittingEarthAtm class that injects neutrino flux during propagation - Uses TriCubicInterpolator for 3D interpolation in (coszen, height, energy) - Supports all 3 flavors (nu_e, nu_mu, nu_tau) plus antineutrinos - Production profile generated from MCEq atmospheric flux calculations - Added GetCosZenith() accessor to EarthAtm::Track - Added default constructor for TriCubicInterpolator New files: - data/atmos_prod/PROD_MODEL_MCEQ.dat - production profile data - data/atmos_prod/generate_production_profile.py - MCEq data generator - examples/Emitting_Atmosphere/ - example demonstrating usage Based on work by ConnorSponsler in PR #41. Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
1 parent 365539c commit e36a597

10 files changed

Lines changed: 76979 additions & 4 deletions

File tree

‎configure‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -930,7 +930,8 @@ EXAMPLES := examples/Single_energy/single_energy \
930930
examples/HDF5_Write_Read/read \
931931
examples/Constant_density_layers/const_dens_layers \
932932
examples/Astrophysical_neutrino_flavor_ratio/astro_flv_ratio \
933-
examples/Decoherence/decoherence
933+
examples/Decoherence/decoherence \
934+
examples/Emitting_Atmosphere/emitting_atmosphere
934935
935936
CXXFLAGS= -std=c++11
936937
@@ -1070,6 +1071,10 @@ examples/Constant_density_layers/const_dens_layers : $(LIBDIR)/$(DYN_PRODUCT) ex
10701071
@echo Compiling constant density layer example
10711072
@$(CXX) $(EXAMPLES_FLAGS) examples/Constant_density_layers/main.cpp -lnuSQuIDS $(LDFLAGS) -o $@
10721073
1074+
examples/Emitting_Atmosphere/emitting_atmosphere : $(LIBDIR)/$(DYN_PRODUCT) examples/Emitting_Atmosphere/main.cpp
1075+
@echo Compiling emitting atmosphere example
1076+
@$(CXX) $(EXAMPLES_FLAGS) examples/Emitting_Atmosphere/main.cpp -lnuSQuIDS $(LDFLAGS) -o $@
1077+
10731078
.PHONY: install uninstall clean test docs benchmark
10741079
clean:
10751080
@echo Erasing generated files

‎data/atmos_prod/PROD_MODEL_MCEQ.dat‎

Lines changed: 76236 additions & 0 deletions
Large diffs are not rendered by default.

‎data/atmos_prod/generate_mceq_flux.py‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
#!/usr/bin/env python3
2+
"""
3+
Generate atmospheric neutrino differential flux production profiles using MCEq.
4+
5+
This script produces the data file needed by EmittingEarthAtm in nuSQuIDS.
6+
Output format: coszen height energy nu_e nu_e_bar nu_mu nu_mu_bar
7+
8+
Requirements:
9+
pip install MCEq crflux
10+
11+
Usage:
12+
python generate_mceq_flux.py
13+
"""
14+
15+
from MCEq.core import config, MCEqRun
16+
import crflux.models as crf
17+
import numpy as np
18+
import sys
19+
20+
# Configuration
21+
OUTPUT_FILE = "PROD_MODEL_MCEQ.dat"
22+
MONTH = 'January'
23+
LOCATION = 'SouthPole'
24+
25+
# Grid parameters
26+
H_PTS = 100 # Number of height points
27+
COSZEN_PTS = 50 # Number of cosine zenith points
28+
H_MIN_CM = 1e2 # 1 meter in cm
29+
H_MAX_CM = 6e6 # 60 km in cm
30+
31+
def main():
32+
print(f"Generating atmospheric flux data...")
33+
print(f" Height points: {H_PTS}")
34+
print(f" Coszen points: {COSZEN_PTS}")
35+
print(f" Location: {LOCATION}, {MONTH}")
36+
37+
# Create grids
38+
coszen_grid = np.linspace(-1, 1, COSZEN_PTS)
39+
zen_grid = np.arccos(coszen_grid)
40+
zen_deg_grid = np.degrees(zen_grid)
41+
42+
# Height grid from 60km to 1m (in cm), logarithmically spaced
43+
h_grid = np.logspace(np.log10(H_MAX_CM), np.log10(H_MIN_CM), H_PTS)
44+
45+
# Initialize MCEq
46+
print("Initializing MCEq...")
47+
mceq = MCEqRun(
48+
interaction_model='SIBYLL23C',
49+
primary_model=(crf.HillasGaisser2012, 'H3a'),
50+
theta_deg=zen_deg_grid[0],
51+
density_model=('MSIS00_IC', (LOCATION, MONTH)),
52+
)
53+
54+
# Open output file
55+
with open(OUTPUT_FILE, 'w') as f:
56+
# Write header comment
57+
f.write("# Atmospheric neutrino differential flux production from MCEq\n")
58+
f.write("# Format: coszen height[cm] energy[GeV] nu_e nu_e_bar nu_mu nu_mu_bar\n")
59+
f.write(f"# Generated with: SIBYLL23C, HillasGaisser2012 H3a, {LOCATION} {MONTH}\n")
60+
f.write(f"# Grid: {COSZEN_PTS} coszen x {H_PTS} heights x {len(mceq.e_grid)} energies\n")
61+
62+
total_iterations = COSZEN_PTS * H_PTS
63+
current = 0
64+
65+
for izen, coszen in enumerate(coszen_grid):
66+
# Set zenith angle
67+
mceq.set_theta_deg(zen_deg_grid[izen])
68+
69+
# Convert heights to slant depths
70+
X_grid = mceq.density_model.h2X(h_grid)
71+
72+
# Solve cascade equations
73+
mceq.solve(int_grid=X_grid)
74+
75+
for idx in range(H_PTS):
76+
# Get fluxes for all flavors (neutrinos and anti-neutrinos)
77+
nue_flux = mceq.get_solution('nue', grid_idx=idx)
78+
nuebar_flux = mceq.get_solution('antinue', grid_idx=idx)
79+
numu_flux = mceq.get_solution('numu', grid_idx=idx)
80+
numubar_flux = mceq.get_solution('antinumu', grid_idx=idx)
81+
82+
for ien, energy in enumerate(mceq.e_grid):
83+
# Write: coszen height energy nu_e nu_e_bar nu_mu nu_mu_bar
84+
f.write(f"{coszen:.6e} {h_grid[idx]:.6e} {energy:.6e} "
85+
f"{nue_flux[ien]:.6e} {nuebar_flux[ien]:.6e} "
86+
f"{numu_flux[ien]:.6e} {numubar_flux[ien]:.6e}\n")
87+
88+
current += 1
89+
if current % 100 == 0:
90+
print(f" Progress: {current}/{total_iterations} ({100*current/total_iterations:.1f}%)")
91+
92+
print(f"Done! Output written to {OUTPUT_FILE}")
93+
print(f"File contains {COSZEN_PTS * H_PTS * len(mceq.e_grid)} data points")
94+
95+
if __name__ == "__main__":
96+
main()

‎data/atmos_prod/generate_production_profile.py‎

Lines changed: 191 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,191 @@
1+
#!/usr/bin/env python3
2+
"""
3+
Generate atmospheric neutrino differential production profiles using MCEq.
4+
5+
This script produces the data file needed by EmittingEarthAtm in nuSQuIDS.
6+
It computes the derivative of the flux with respect to height to get the
7+
production rate at each point in the atmosphere.
8+
9+
Output format: coszen height[cm] energy[GeV] nu_e nu_e_bar nu_mu nu_mu_bar
10+
Where the flux values are differential production rates (d(flux)/d(height)).
11+
12+
Requirements:
13+
pip install MCEq crflux numpy
14+
15+
Usage:
16+
python generate_production_profile.py [--output OUTPUT_FILE]
17+
"""
18+
19+
from MCEq.core import config, MCEqRun
20+
import crflux.models as crf
21+
import numpy as np
22+
import argparse
23+
import sys
24+
25+
# Default configuration
26+
DEFAULT_OUTPUT = "PROD_MODEL_MCEQ.dat"
27+
MONTH = 'January'
28+
LOCATION = 'SouthPole'
29+
30+
# Grid parameters
31+
H_PTS = 100 # Number of height points
32+
COSZEN_PTS = 50 # Number of cosine zenith points
33+
H_MIN_CM = 1e2 # 1 meter in cm
34+
H_MAX_CM = 6e6 # 60 km in cm
35+
36+
def compute_differential(flux_array, height_array):
37+
"""
38+
Compute the differential production rate from integrated flux.
39+
40+
Uses central differences for interior points and forward/backward
41+
differences at boundaries.
42+
43+
Parameters:
44+
flux_array: 1D array of flux values at each height
45+
height_array: 1D array of heights (in cm)
46+
47+
Returns:
48+
1D array of differential production rates
49+
"""
50+
n = len(flux_array)
51+
diff = np.zeros(n)
52+
53+
# Use log-spacing aware differentiation
54+
log_h = np.log(height_array)
55+
56+
for i in range(n):
57+
if i == 0:
58+
# Forward difference at top of atmosphere
59+
dflux = flux_array[i+1] - flux_array[i]
60+
dlogh = log_h[i+1] - log_h[i]
61+
# d(flux)/d(h) = d(flux)/d(log h) * d(log h)/d(h) = d(flux)/d(log h) / h
62+
diff[i] = dflux / (dlogh * height_array[i]) if dlogh != 0 else 0
63+
elif i == n - 1:
64+
# Backward difference at ground
65+
dflux = flux_array[i] - flux_array[i-1]
66+
dlogh = log_h[i] - log_h[i-1]
67+
diff[i] = dflux / (dlogh * height_array[i]) if dlogh != 0 else 0
68+
else:
69+
# Central difference for interior points
70+
dflux = flux_array[i+1] - flux_array[i-1]
71+
dlogh = log_h[i+1] - log_h[i-1]
72+
diff[i] = dflux / (dlogh * height_array[i]) if dlogh != 0 else 0
73+
74+
return diff
75+
76+
def main():
77+
parser = argparse.ArgumentParser(description='Generate MCEq atmospheric production profile')
78+
parser.add_argument('--output', '-o', default=DEFAULT_OUTPUT,
79+
help=f'Output filename (default: {DEFAULT_OUTPUT})')
80+
parser.add_argument('--coszen-pts', type=int, default=COSZEN_PTS,
81+
help=f'Number of cosine zenith points (default: {COSZEN_PTS})')
82+
parser.add_argument('--height-pts', type=int, default=H_PTS,
83+
help=f'Number of height points (default: {H_PTS})')
84+
args = parser.parse_args()
85+
86+
print(f"Generating atmospheric neutrino production profile...")
87+
print(f" Height points: {args.height_pts}")
88+
print(f" Coszen points: {args.coszen_pts}")
89+
print(f" Location: {LOCATION}, {MONTH}")
90+
print(f" Output: {args.output}")
91+
92+
# Create grids
93+
coszen_grid = np.linspace(-1, 1, args.coszen_pts)
94+
zen_grid = np.arccos(coszen_grid)
95+
zen_deg_grid = np.degrees(zen_grid)
96+
97+
# Height grid from 60km to ground (in cm), logarithmically spaced
98+
# Going from high altitude to low altitude
99+
h_grid = np.logspace(np.log10(H_MAX_CM), np.log10(H_MIN_CM), args.height_pts)
100+
101+
# Initialize MCEq
102+
print("Initializing MCEq...")
103+
mceq = MCEqRun(
104+
interaction_model='SIBYLL23C',
105+
primary_model=(crf.HillasGaisser2012, 'H3a'),
106+
theta_deg=zen_deg_grid[0],
107+
density_model=('MSIS00_IC', (LOCATION, MONTH)),
108+
)
109+
110+
e_grid = mceq.e_grid
111+
n_energies = len(e_grid)
112+
113+
print(f" Energy points: {n_energies} ({e_grid[0]:.2e} - {e_grid[-1]:.2e} GeV)")
114+
115+
# Storage for all data
116+
all_data = []
117+
118+
total_zenith = args.coszen_pts
119+
for izen, coszen in enumerate(coszen_grid):
120+
print(f" Processing zenith {izen+1}/{total_zenith} (cos(zen) = {coszen:.3f})...")
121+
122+
# Set zenith angle
123+
mceq.set_theta_deg(zen_deg_grid[izen])
124+
125+
# Convert heights to slant depths
126+
X_grid = mceq.density_model.h2X(h_grid)
127+
128+
# Solve cascade equations
129+
mceq.solve(int_grid=X_grid)
130+
131+
# Get fluxes at all heights for each flavor
132+
# Shape: (n_heights, n_energies)
133+
nue_flux = np.array([mceq.get_solution('nue', grid_idx=idx) for idx in range(args.height_pts)])
134+
nuebar_flux = np.array([mceq.get_solution('antinue', grid_idx=idx) for idx in range(args.height_pts)])
135+
numu_flux = np.array([mceq.get_solution('numu', grid_idx=idx) for idx in range(args.height_pts)])
136+
numubar_flux = np.array([mceq.get_solution('antinumu', grid_idx=idx) for idx in range(args.height_pts)])
137+
nutau_flux = np.array([mceq.get_solution('nutau', grid_idx=idx) for idx in range(args.height_pts)])
138+
nutaubar_flux = np.array([mceq.get_solution('antinutau', grid_idx=idx) for idx in range(args.height_pts)])
139+
140+
# Compute differential production for each energy
141+
for ie in range(n_energies):
142+
# Get flux column for this energy
143+
nue_col = nue_flux[:, ie]
144+
nuebar_col = nuebar_flux[:, ie]
145+
numu_col = numu_flux[:, ie]
146+
numubar_col = numubar_flux[:, ie]
147+
nutau_col = nutau_flux[:, ie]
148+
nutaubar_col = nutaubar_flux[:, ie]
149+
150+
# Compute differentials
151+
dnue = compute_differential(nue_col, h_grid)
152+
dnuebar = compute_differential(nuebar_col, h_grid)
153+
dnumu = compute_differential(numu_col, h_grid)
154+
dnumubar = compute_differential(numubar_col, h_grid)
155+
dnutau = compute_differential(nutau_col, h_grid)
156+
dnutaubar = compute_differential(nutaubar_col, h_grid)
157+
158+
# Store data for each height
159+
for ih in range(args.height_pts):
160+
all_data.append([
161+
coszen,
162+
h_grid[ih],
163+
e_grid[ie],
164+
dnue[ih],
165+
dnuebar[ih],
166+
dnumu[ih],
167+
dnumubar[ih],
168+
dnutau[ih],
169+
dnutaubar[ih]
170+
])
171+
172+
# Write output file
173+
print(f"Writing {len(all_data)} data points to {args.output}...")
174+
with open(args.output, 'w') as f:
175+
# Write header
176+
f.write("# Atmospheric neutrino differential production profile from MCEq\n")
177+
f.write("# Format: coszen height[cm] energy[GeV] nu_e nu_e_bar nu_mu nu_mu_bar nu_tau nu_tau_bar\n")
178+
f.write("# Values are differential production rates: d(flux)/d(height)\n")
179+
f.write(f"# Generated with: SIBYLL23C, HillasGaisser2012 H3a, {LOCATION} {MONTH}\n")
180+
f.write(f"# Grid: {args.coszen_pts} coszen x {args.height_pts} heights x {n_energies} energies\n")
181+
f.write("#\n")
182+
183+
for row in all_data:
184+
f.write(f"{row[0]:.6e} {row[1]:.6e} {row[2]:.6e} "
185+
f"{row[3]:.6e} {row[4]:.6e} {row[5]:.6e} {row[6]:.6e} "
186+
f"{row[7]:.6e} {row[8]:.6e}\n")
187+
188+
print(f"Done! Output written to {args.output}")
189+
190+
if __name__ == "__main__":
191+
main()

‎examples/Emitting_Atmosphere/emitting_atmosphere_fluxes.png‎

42.4 KB
Loading

0 commit comments

Comments
 (0)