diff options
Diffstat (limited to 'tutorials/module_4/spectroscopy_problem/spectroscopy.py')
| -rw-r--r-- | tutorials/module_4/spectroscopy_problem/spectroscopy.py | 110 |
1 files changed, 94 insertions, 16 deletions
diff --git a/tutorials/module_4/spectroscopy_problem/spectroscopy.py b/tutorials/module_4/spectroscopy_problem/spectroscopy.py index f66734a..404e5ca 100644 --- a/tutorials/module_4/spectroscopy_problem/spectroscopy.py +++ b/tutorials/module_4/spectroscopy_problem/spectroscopy.py @@ -23,6 +23,7 @@ import pandas as pd import matplotlib.pyplot as plt + Hg_data = 'Lampa_Calibrare_Mercur.xlsx' #Ar_data = 'Spectru Descarcare Argon.xlsx' Ox_data = 'Calibrare Intensitate Oxigen.xlsx' @@ -35,45 +36,122 @@ df_Ox = pd.read_excel(Ox_data, header=3, engine='openpyxl') # Trims mercury data to the first 2000 pixels df_Hg = df_Hg[df_Hg['Pixels']<2000] +# %% Plot Intensity vs Pixels for mercury -"Plot Intensity vs Pixels for mercury" +# Add Vertical lines to help match pixel count +wavelength_true = np.array([365.015, 404.656, 407.783, 435.833, 546.074, 576.960, 579.006]) +bounds = np.array([635,1451]) +m = (bounds[-1]-bounds[0])/(wavelength_true[-1]-wavelength_true[0]) +c = bounds[0]-(m*wavelength_true[0]) +pixel_count = m * wavelength_true + c # Plot Hg intensity-pixel plot plt.figure(figsize=(8,5)) plt.plot(df_Hg['Pixels'], df_Hg['Intensity'], linestyle='-') plt.xlabel('Pixels') plt.ylabel('Intensity [a.u.]') +plt.vlines(pixel_count,0,17000,colors='r') plt.title('Pixel-Intensity (Hg)') plt.grid(True) plt.show() -"Calibrate length dimension" -C = np.polyfit(df_Hg['Pixels'],df_Hg['Intensity'],3) -lambda_Hg = lambda p : C[3] + C[2]*p+C[1]*p**2+C[0]*p**3 +# %% Calibrate length dimension using Mercury +# Arrays to be interpolated +wavelength_peak_nm = np.array([365.015, 404.656, 435.833, 546.074, 576.960, 579.006]) # From Literature +pixel_peak = np.array([635, 784, 902, 1324, 1443, 1451]) + +#Polynomial coefficients +C = np.polyfit(pixel_peak,wavelength_peak_nm,3) -lambda_calibrated = lambda_Hg(df_Hg['Pixel']) +# Calibration function +wavelength_Hg = lambda p : C[3] + C[2]*p+C[1]*p**2+C[0]*p**3 +# or +#def wavelength_Hg(p): +# return C[3] + C[2]*p+C[1]*p**2+C[0]*p**3 +wavelength_calibrated = wavelength_Hg(df_Hg['Pixels']) # Plot Hg intensity-wavelength plot plt.figure(figsize=(8,5)) -plt.plot(lambda_Hg, df_Hg['Intensity'], linestyle='-') +plt.plot(wavelength_calibrated, df_Hg['Intensity'], linestyle='-') plt.xlabel('Wavelength [nm]') plt.ylabel('Intensity [a.u.]') plt.title('Wavelength-Intensity (Hg)') plt.grid(True) plt.show() +# %% Calibrate intensity dimension using Tungsten +# Now that we have calibrated the wavelength dimension can can now use +# a tungsten source to calibrate intensity. +wavelength = df_Ox['Wavelength [nm]'] +I_W_measured = df_Ox['I_Tungsten [a.u.]']-df_Ox['I_Background [a.u.]'] -""" -"Calibrate intensity dimension" -epsilon = # Surface emissivity -lambd = - -I_W_true = epsilon*((2hc**2)/(lambd**5)*1/(e**(hc/kT)-1)) -I_W_meas= df_Ox['I_Tungsten [a.u.]'] -R=I_W_meas/I_W_true -I_plasma_meas=R*I_plasma_true -"""
\ No newline at end of file +# Calculate theoretical Intensity of Tungsten using Planks eqn. (TO BE CHECKED) +epsilon = 0.4 +h = 6.626e-34 # J/s +c = 2.99e8 # m/s +k = 1.381e-23 # J/K +T = 1800 # K + +I_planck = lambda wavelength0 : epsilon*((2*h*c**2)/(wavelength0**5)*1/(np.exp(h*c/(k*T))-1)) +I_W_true = I_planck(wavelength) + +#Response Function +R=I_W_measured/I_W_true +#I_plasma_meas=R*I_plasma_true + +# Plot Response function vs lambda +plt.figure(figsize=(8,5)) +plt.plot(wavelength, R, linestyle='-') +plt.xlabel('Wavelength [nm]') +plt.ylabel('R') +plt.title('Response function') +plt.grid(True) +plt.show() + +# Calculate Intensity of Oxygen +I_Ox_measured = df_Ox['I_Plasma Oxygen [a.u.]']-df_Ox['I_Background [a.u.]'] + +I_Ox_true = R * I_Ox_measured + + +# Plot Hg intensity-wavelength plot +plt.figure(figsize=(8,5)) +plt.plot(wavelength, I_Ox_measured, linestyle='-') +plt.xlabel('Wavelength [$nm$]') +plt.ylabel('Intensity [a.u,]') +plt.xlim(750,870) +plt.title('Wavelength-Intensity ($Ox_{measured}$)') +plt.grid(True) +plt.show() + +# Plot Hg intensity-wavelength plot +plt.figure(figsize=(8,5)) +plt.plot(wavelength, I_Ox_true, linestyle='-') +plt.xlabel('Wavelength [$nm$]') +plt.ylabel('Intensity [$W/m^2/sr/nm$]') +plt.xlim(750,870) +plt.title('Wavelength-Intensity ($Ox_{true}$)') +plt.grid(True) +plt.show() + + +# %% Calculate Dencity of Ox at upper state + +from scipy.interpolate import interp1d +from scipy.integrate import quad + +# Create an interpolated callable function +I_Ox_true = interp1d(wavelength, I_Ox_true, kind='linear', + fill_value=0, bounds_error=False) + +# Now integrate +delta = 1 +peak = 846.5 +bound_lower = peak - delta +bound_upper = peak + delta + +numerator, err = quad(I_Ox_true, bound_lower, bound_upper) |