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")