from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
List of constants
Sr_mol = 1e-6
Avogadro = 6.02e23
Sr_nuclei = Sr_mol * Avogadro
print(Sr_nuclei)
nuclear_data_1 = Path(“Nuclear_data_1.csv”)
nuclear_data_2 = Path(“Nuclear_data_2.csv”)
nuclear_plot_filename = “nuclear_decay_fit.png”
def activity_rb(t_sec, lambda_Rb, lambda_Sr):
return Sr_nuclei * lambda_Rb * lambda_Sr * np.exp(-lambda_Sr * t_sec) / (lambda_Sr - lambda_Rb) * (np.exp(-lambda_Rb * t_sec) - np.exp(-lambda_Sr * t_sec))
def activity_Rb(t_sec, lambda_Rb, lambda_Sr):
return Sr_nuclei * lambda_Rb * lambda_Sr * np.exp(-lambda_Sr * t_sec) / (lambda_Sr - lambda_Rb) * (np.exp(-lambda_Rb * t_sec) - np.exp(-lambda_Sr * t_sec))
Data reading, cleaning, and removing of outliers
def read_and_clean_file(file_path: Path) → pd.DataFrame:
try:
df = pd.read_csv(
file_path,
sep=None,
engine=“python”,
comment=“%”,
header=None
)
except FileNotFoundError:
print(f"Error: File ‘{file_path}’ not found.“)
return pd.DataFrame()
except pd.errors.EmptyDataError:
print(f"Error: File ‘{file_path}’ is empty.”)
return pd.DataFrame()
except Exception as exc:
print(f"Error reading ‘{file_path}’: {exc}")
return pd.DataFrame()
if df.shape[1] < 3:
print(f"Error: File '{file_path}' does not have at least 3 columns.")
return pd.DataFrame()
df = df.iloc[:, :3]
df.columns = ["time", "activity", "uncertainty"]
df = df.apply(pd.to_numeric, errors="coerce")
df.dropna(inplace=True)
df = df[(df["time"] >= 0) & (df["activity"] > 0) & (df["uncertainty"] > 0)]
return df
def remove_outliers_iqr(df: pd.DataFrame) → pd.DataFrame:
Q1 = df[“activity”].quantile(0.25)
Q3 = df[“activity”].quantile(0.75)
IQR = Q3 - Q1
lower_quartile = Q1 - 1.5 * IQR
upper_quartile = Q3 + 1.5 * IQR
mask = (df["activity"] >= lower_quartile) & (df["activity"] <= upper_quartile)
return df[mask]
Fit and stats
def fit_parameters(times_h: np.ndarray, activities_tbq: np.ndarray, uncertainties_tbq: np.ndarray):
time_sec = times_h * 3600.0
activity_bq = activities_tbq * 1e12
activity_uncertainty_bq = uncertainties_tbq * 1e12
popt, pcov = curve_fit(
activity_rb,
time_sec,
activity_bq,
sigma=activity_uncertainty_bq,
absolute_sigma=True,
p0=[5e-4, 5e-3],
maxfev=20000,
)
perr = np.sqrt(np.diag(pcov))
model = activity_rb(time_sec, *popt)
chi2 = np.sum(((activity_bq - model) / activity_uncertainty_bq) ** 2)
chi2_red = chi2 / (len(activity_bq) - 2)
return popt, perr, pcov, chi2, chi2_red
def half_life(lambda_value: float) → float:
return np.log(2.0) / lambda_value
def half_life_unc(lambda_value: float, sigma_lambda: float) → float:
return (np.log(2.0) / lambda_value**2) * sigma_lambda
def predict_with_uncertainty(popt, pcov, t_min: float):
t_sec = t_min * 60.0
lambda_Rb, lambda_Sr = popt
eps_Rb = lambda_Rb * 1e-6
eps_Sr = lambda_Sr * 1e-6
df_dRb = (
activity_Rb([t_sec], lambda_Rb + eps_Rb, lambda_Sr)[0]
- activity_Rb([t_sec], lambda_Rb - eps_Rb, lambda_Sr)[0]
) / (2.0 * eps_Rb)
df_dSr = (
activity_Rb([t_sec], lambda_Rb, lambda_Sr + eps_Sr)[0]
- activity_Rb([t_sec], lambda_Rb, lambda_Sr - eps_Sr)[0]
) / (2.0 * eps_Sr)
return()
Plotting
def make_publication_plot(times_h: np.ndarray, activities_tbq: np.ndarray, uncertainties_tbq: np.ndarray, popt, chi2_red: float) None:
plt.style.use(“'Solarize_Light2”)
fig, ax = plt.subplots(figsize=(9, 6), dpi=350)
ax.errorbar(
times_h,
activities_tbq,
yerr=uncertainties_tbq,
fmt="x",
ms=4,
capsize=3,
color="darkblue",
ecolor="black",
elinewidth=1,
label="Data",
)
t_plot = np.linspace(times_h.min(), times_h.max(), 600)
model_tbq = activity_rb(t_plot * 3600.0, *popt) / 1e12
ax.plot(
t_plot,
model_tbq,
color="red",
lw=2.2,
label="Best-fit model",
)
ax.set_xlabel("Time (hours)", fontsize=12)
ax.set_ylabel("Activity of 79Rb (TBq)", fontsize=12)
ax.set_title("Nuclear Decay of 79Sr into 79Rb", fontsize=14)
lambda_rb, lambda_sr = popt
text = (
f"λ_Rb = {lambda_rb:.3e} s⁻¹\n"
f"λ_Sr = {lambda_sr:.3e} s⁻¹\n"
f"χ²_red = {chi2_red:.2f}"
)
ax.text(
0.97,
0.97,
text,
transform=ax.transAxes,
ha="right",
va="top",
fontsize=10,
bbox=dict(facecolor="white", alpha=0.9),
)
ax.legend()
fig.tight_layout()
fig.savefig(nuclear_plot_filename)
plt.close(fig)
plt.show()
Main
def main(): None:
“”"
Load, read and clean the data files
“”"
df1 = read_and_clean_file(nuclear_data_1)
df2 = read_and_clean_file(nuclear_data_2)
if df1.empty or df2.empty:
print("One or both files is not able to be processed.")
return
"""
Combine the two sets of data
"""
combined_df = pd.concat([df1, df2], ignore_index=True)
combined_df.sort_values(by="time", inplace=True)
print(f"Original combined data shape: {combined_df.shape}")
"""
Remove outliers
"""
cleaned_df = remove_outliers_iqr(combined_df)
cleaned_df.reset_index(drop=True, inplace=True)
print(f"Cleaned data shape (outliers removed): {cleaned_df.shape}")
"""
Extract the arrays
"""
times_h = cleaned_df["time"].to_numpy()
activities_tbq = cleaned_df["activity"].to_numpy()
uncertainties_tbq = cleaned_df["uncertainty"].to_numpy()
"""
Fit model (chi squared)
"""
popt, perr, pcov, chi2, chi2_red = fit_parameters(
times_h, activities_tbq, uncertainties_tbq
)
lambda_Rb, lambda_Sr = popt
sigma_Rb, sigma_Sr = perr
"""
Convert decay constants to minutes ^-1
"""
lambda_Rb_min = lambda_Rb * 60.0
lambda_Sr_min = lambda_Sr * 60.0
sigma_Rb_min = sigma_Rb * 60.0
sigma_Sr_min = sigma_Sr * 60.0
"""
Half-lives in seconds
"""
t_half_Rb_sec = half_life(lambda_Rb)
t_half_Sr_sec = half_life(lambda_Sr)
dt_half_Rb_sec = half_life_unc(lambda_Rb, sigma_Rb)
dt_half_Sr_sec = half_life_unc(lambda_Sr, sigma_Sr)
"""
Conversion of half-lives in seconds to minutes
"""
t_half_Rb_min = t_half_Rb_sec / 60.0
t_half_Sr_min = t_half_Sr_sec / 60.0
dt_half_Rb_min = dt_half_Rb_sec / 60.0
dt_half_Sr_min = dt_half_Sr_sec / 60.0
"""
Prediction at 90 minutes
"""
pred_bq, pred_sigma_bq = predict_with_uncertainty(popt, pcov, 90.0)
pred_tbq = pred_bq / 1e12
pred_sigma_tbq = pred_sigma_bq / 1e12
"""
Print results
"""
print("\n=== Best-fit decay constants (in minutes ^-1) ===")
print(f"λ_Rb = {lambda_Rb_min:.3g} ± {sigma_Rb_min:.3g} min⁻¹")
print(f"λ_Sr = {lambda_Sr_min:.3g} ± {sigma_Sr_min:.3g} min⁻¹")
print("Half-lives (in minutes)")
print(f"t_1/2,Rb = {t_half_Rb_min:.3g} ± {dt_half_Rb_min:.3g} min")
print(f"t_1/2,Sr = {t_half_Sr_min:.3g} ± {dt_half_Sr_min:.3g} min")
print("Chi-squared")
print(f"χ^2 = {chi2:.2f}")
print(f"χ^2_red = {chi2_red:.2f}")
print("\n=== Prediction at t = 90 minutes ===")
print(f"A(90 min) = {pred_tbq:.3g} ± {pred_sigma_tbq:.2g} TBq")
"""
Summary table of data
"""
print("\n==================== SUMMARY TABLE ====================")
print(f"{'Quantity':<30} {'Value':>15} {'Uncertainty':>15}")
print("--------------------------------------------------------")
print(f"{'λ_Rb (minutes ^-1)':<30} {lambda_Rb_min:.3g:>15} {sigma_Rb_min:.3g:>15}")
print(f"{'λ_Sr (minutes ^-1)':<30} {lambda_Sr_min:.3g:>15} {sigma_Sr_min:.3g:>15}")
print(f"{'t_1/2 Rb (min)':<30} {t_half_Rb_min:.3g:>15} {dt_half_Rb_min:.3g:>15}")
print(f"{'t_1/2 Sr (min)':<30} {t_half_Sr_min:.3g:>15} {dt_half_Sr_min:.3g:>15}")
print(f"{'χ^2_red':<30} {chi2_red:.2f:>15} {'—':>15}")
print("========================================================\n")
"""
Data analysis
"""
print("\n==================== DATA ANALYSIS ====================")
print(f"{'Quantity':<30} {'Rel. Uncertainty':>20} {'Status':>12} {'Explanation':>30}")
print("--------------------------------------------------------------------------")
def classify_uncertainty(value, sigma):
rel = abs(sigma / value)
if rel < 0.20:
status = "GOOD"
explanation = "uncertainty < 20%"
elif rel <= 0.25:
status = "ACCEPTABLE"
explanation = "uncertainty between 20-25%"
else:
status = "POOR"
explanation = "uncertainty > 25%"
return rel, status, explanation
# Compute classifications
rel_Rb_lambda, status_Rb_lambda, expl_Rb_lambda = classify_uncertainty(lambda_Rb_min, sigma_Rb_min)
rel_Sr_lambda, status_Sr_lambda, expl_Sr_lambda = classify_uncertainty(lambda_Sr_min, sigma_Sr_min)
rel_Rb_half, status_Rb_half, expl_Rb_half = classify_uncertainty(t_half_Rb_min, dt_half_Rb_min)
rel_Sr_half, status_Sr_half, expl_Sr_half = classify_uncertainty(t_half_Sr_min, dt_half_Sr_min)
# Print table
print(f"{'λ_Rb (min⁻¹)':<30} {rel_Rb_lambda:>20.3f} {status_Rb_lambda:>12} {expl_Rb_lambda:>30}")
print(f"{'λ_Sr (min⁻¹)':<30} {rel_Sr_lambda:>20.3f} {status_Sr_lambda:>12} {expl_Sr_lambda:>30}")
print(f"{'t_1/2 Rb (min)':<30} {rel_Rb_half:>20.3f} {status_Rb_half:>12} {expl_Rb_half:>30}")
print(f"{'t_1/2Sr (min)':<30} {rel_Sr_half:>20.3f} {status_Sr_half:>12} {expl_Sr_half:>30}")
print(f"{'χ^2_red':<30} {'N/A':>20} {'N/A':>12} {'no uncertainty':>30}")
print("==========================================================================\n")
"""
Plot
"""
make_publication_plot(times_h, activities_tbq, uncertainties_tbq, popt, chi2_red)
print(f"Plot saved as '{nuclear_plot_filename}'.")