import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
import math
# Lead
Energies = np.linspace(1,100,199)
LAC = np.array([
5.87E+04, 2.65E+04, 1.44E+04, 1.65E+04, 2.21E+04,
1.74E+04, 1.40E+04, 1.05E+04, 8.16E+03, 6.45E+03,
5.20E+03, 4.25E+03, 3.52E+03, 2.96E+03, 2.52E+03,
2.15E+03, 1.86E+03, 1.62E+03, 1.42E+03, 1.25E+03,
1.11E+03, 9.88E+02, 8.85E+02, 7.96E+02, 7.19E+02,
1.63E+03, 1.47E+03, 1.34E+03, 1.22E+03, 1.57E+03,
1.68E+03, 1.55E+03, 1.44E+03, 1.33E+03, 1.24E+03,
1.16E+03, 1.08E+03, 1.01E+03, 9.50E+02, 8.90E+02,
8.35E+02, 7.85E+02, 7.39E+02, 6.97E+02, 6.57E+02,
6.21E+02, 5.88E+02, 5.57E+02, 5.28E+02, 5.01E+02,
4.76E+02, 4.53E+02, 4.31E+02, 4.11E+02, 3.92E+02,
3.74E+02, 3.57E+02, 3.42E+02, 3.27E+02, 3.13E+02,
3.00E+02, 2.87E+02, 2.75E+02, 2.64E+02, 2.54E+02,
2.44E+02, 2.34E+02, 2.25E+02, 2.17E+02, 2.09E+02,
2.01E+02, 1.94E+02, 1.87E+02, 1.80E+02, 1.74E+02,
1.68E+02, 1.62E+02, 1.57E+02, 1.52E+02, 1.47E+02,
1.42E+02, 1.38E+02, 1.33E+02, 1.29E+02, 1.25E+02,
1.21E+02, 1.18E+02, 1.14E+02, 1.11E+02, 1.07E+02,
1.04E+02, 1.01E+02, 9.85E+01, 9.57E+01, 9.31E+01,
9.06E+01, 8.81E+01, 8.57E+01, 8.35E+01, 8.13E+01,
7.91E+01, 7.71E+01, 7.51E+01, 7.32E+01, 7.14E+01,
6.96E+01, 6.79E+01, 6.62E+01, 6.46E+01, 6.31E+01,
6.16E+01, 6.01E+01, 5.87E+01, 5.73E+01, 5.60E+01,
5.48E+01, 5.35E+01, 5.23E+01, 5.12E+01, 5.01E+01,
4.90E+01, 4.79E+01, 4.69E+01, 4.59E+01, 4.49E+01,
4.40E+01, 4.31E+01, 4.22E+01, 4.13E+01, 4.05E+01,
3.97E+01, 3.89E+01, 3.81E+01, 3.74E+01, 3.66E+01,
3.59E+01, 3.52E+01, 3.46E+01, 3.39E+01, 3.33E+01,
3.27E+01, 3.21E+01, 3.15E+01, 3.09E+01, 3.04E+01,
2.98E+01, 2.93E+01, 2.88E+01, 2.83E+01, 2.78E+01,
2.73E+01, 2.68E+01, 2.64E+01, 2.59E+01, 2.55E+01,
2.51E+01, 2.47E+01, 2.43E+01, 2.39E+01, 2.35E+01,
2.31E+01, 2.27E+01, 2.24E+01, 2.20E+01, 2.17E+01,
2.13E+01, 2.10E+01, 2.07E+01, 2.04E+01, 2.01E+01,
1.98E+01, 1.95E+01, 1.92E+01, 1.89E+01, 1.86E+01,
8.26E+01, 8.15E+01, 8.03E+01, 7.91E+01, 7.80E+01,
7.69E+01, 7.58E+01, 7.47E+01, 7.37E+01, 7.27E+01,
7.17E+01, 7.07E+01, 6.98E+01, 6.88E+01, 6.79E+01,
6.70E+01, 6.61E+01, 6.52E+01, 6.44E+01, 6.35E+01,
6.27E+01, 6.19E+01, 6.11E+01, 6.03E+01
])
I_0 = np.linspace(1,100,199)
P_0 = np.array([
0,
0,
5.11347E-06,
0.000121555,
3.08604E-05,
0,
0,
8.31027E-06,
8.7771E-05,
0.000359376,
0.000946079,
0.001877169,
0.003548272,
0.005150231,
0.006742632,
0.016790123,
0.014374974,
0.012801486,
0.024634899,
0.011637658,
0.015204991,
0.0107191,
0.008945965,
0.009214431,
0.015634477,
0.009179067,
0.009152242,
0.00943935,
0.010478518,
0.010235272,
0.010345676,
0.010531003,
0.010234591,
0.01073557,
0.010859254,
0.011212683,
0.011140717,
0.011089523,
0.010867856,
0.010987345,
0.011155824,
0.010809004,
0.011037909,
0.010752565,
0.010937409,
0.010729171,
0.010745011,
0.010583666,
0.010481057,
0.010277613,
0.010390617,
0.010158208,
0.010171552,
0.009960702,
0.010035846,
0.009819078,
0.009828614,
0.009517127,
0.009333436,
0.009406346,
0.009269727,
0.008867548,
0.008940091,
0.008521914,
0.008432146,
0.008507657,
0.008371625,
0.008354725,
0.00794921,
0.007963582,
0.007768016,
0.007666614,
0.00774671,
0.007514888,
0.007287714,
0.007189879,
0.007027001,
0.006843111,
0.006814986,
0.006970499,
0.006564386,
0.006477514,
0.006288494,
0.006329008,
0.006155955,
0.006046507,
0.005918102,
0.005944737,
0.00578743,
0.005658994,
0.00563903,
0.005555042,
0.005472806,
0.005113881,
0.005375275,
0.005027333,
0.005019937,
0.005089112,
0.004891952,
0.004766316,
0.004586465,
0.004527969,
0.004351737,
0.004400257,
0.004311621,
0.004322532,
0.004149311,
0.004076936,
0.003974694,
0.004101935,
0.003906201,
0.003932123,
0.003898879,
0.003951311,
0.003648835,
0.008549116,
0.003395487,
0.012448897,
0.003362305,
0.003395445,
0.003248954,
0.003229609,
0.003147635,
0.002916915,
0.002987517,
0.002910967,
0.002877974,
0.00273077,
0.002851213,
0.002961867,
0.002615342,
0.002670858,
0.002596836,
0.005511989,
0.002374152,
0.002347264,
0.002388104,
0.002970019,
0.002251517,
0.001893168,
0.001914275,
0.002046383,
0.001923139,
0.001900081,
0.00176021,
0.002097179,
0.001679893,
0.001536802,
0.001690877,
0.002112264,
0.001642358,
0.001356783,
0.00140141,
0.001483782,
0.001486027,
0.001390132,
0.001416076,
0.001205739,
0.001301917,
0.001234389,
0.001154041,
0.001219586,
0.00121325,
0.001145995,
0.001201994,
0.001040079,
0.001012389,
0.000876481,
0.001004057,
0.000925759,
0.000918481,
0.000809436,
0.000821415,
0.000789026,
0.000813591,
0.00068724,
0.000703664,
0.000756499,
0.000682395,
0.000586809,
0.000578912,
0.00049449,
0.000509462,
0.000464834,
0.00045288,
0.000404728,
0.000394017,
0.000356244,
0.000352201,
0.000281882,
0.000271873,
0.00022954,
0.000237268,
0.000188778,
0.000170614,
0.000129775,
9.31903E-05,
6.61896E-05,
2.2146E-05
])
thickness = 0.08
P_final = np.zeros(199)
for i in range(len(P_0)):
P_final[i] = P_0[i] * np.exp(-LAC[i] * thickness)
print(P_final)
# Calculate mean energies based on weighted averages
MeanFinalE = np.sum(Energies * P_final) / np.sum(P_final) if np.sum(P_final) > 0 else 0
MeanInitialE = np.sum(Energies * P_0) / np.sum(P_0) if np.sum(P_0) > 0 else 0
# Output results
print(f"Mean Initial Energy: {MeanInitialE:.2f} keV")
print(f"Mean Final Energy: {MeanFinalE:.2f} keV")