#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Nov  5 12:58:59 2025

@author: christian
"""


"""
Problem:

- Import xls data into Python
- Plot the Intensity [a.u.]  vs pixels
- Interpolate and convert x-axis from pixels to nm (true wavelength) using Hg lamp data (using data in file: **Lampa_Calibrare_Mercur.xlsx**)
- Find response function of the spectrometer using the tungsten lamp data from file: "**Calibrare Intensitate Oxigen.xlsx**)": $R=\frac{I_{measured}}{I_{true}}$ (where True is computed by Planck's law of radiation (see notes in the pptx above)
- Convert y-axis from Intensity [a.u.] into Intensity in [W/(cm^2*sr*nm)] by dividing the measured Oxygen spectrum with the response function: $I_{oxygen, true}=\frac{I_{oxygen, measured}}{R}$
- Once the spectra is in real units: compute the density of one of the oxygen lines by integrating underneath one of the peaks (see equation from Slide 39 - bottom). We will give all of the constants that are in this equation (see the "I**ntensity_Calibration_Oxygen_Discharge_Solution.xlsx**")
"""

import numpy as np
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'

df_Hg = pd.read_excel(Hg_data, header=11, engine='openpyxl')
df_Ox = pd.read_excel(Ox_data, header=3, engine='openpyxl')
#df_Ar = pd.read_excel(Ar_data, header=14, engine='openpyxl')


# Trims mercury data to the first 2000 pixels
df_Hg = df_Hg[df_Hg['Pixels']<2000]

# %% 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 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)

# 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(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.]']

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