import numpy as np
from scipy import constants as sc
from scipy.interpolate import interp1d
from pathlib import Path
from scipy.special import erf as Erf
import pandas as pd
import sys
import os
import csv

class OpticalMedium():

	available_media = list()
	available_media.append("Rhodamine6G")

	def __init__(self, optical_medium):

		"""
			Initiazies an optical medium object.

			Parameters:

				optical_medium (str):		Optical medium

		"""

		if not type(optical_medium) == str:
			raise Exception("optical_medium is expected to be a string")

		if not optical_medium in self.available_media:
			raise Exception(optical_medium+" is an unknown optical medium")

		if optical_medium == "Rhodamine6G":
			self.medium = Rhodamine6G()


	def get_rates(self, lambdas, **kwargs):

		"""
			Calculates the rates of absorption and emission, for a specific optical medium.

			Parameters:

				lambdas (list, or other iterable): Wavelength points where the rates are to be calculated. Wavelength is in meters
				other medium specific arguments
				
		"""

		return self.medium.get_rates(lambdas=lambdas, **kwargs)




class Rhodamine6G(OpticalMedium):

	def __init__(self):
		pass


	def get_rates(self, lambdas, dye_concentration, n):

		"""
			Rates for Rhodamine 6G

			Parameters:

				lambdas (list, or other iterable):  	Wavelength points where the rates are to be calculated. Wavelength is in meters
				dye_concentration (float):				In mM (milimolar) 1 mM = 1 mol / m^3			
				n (float): 								index of refraction

		"""

		# absorption data
		min_wavelength = 480
		max_wavelength = 650
		absorption_spectrum_datafile = Path("data") / 'absorption_cross_sections_R6G_in_EthyleneGlycol_corrected.csv'
		absorption_spectrum_datafile = Path(os.path.dirname(os.path.abspath(__file__))) / absorption_spectrum_datafile
		raw_data2 = pd.read_csv(absorption_spectrum_datafile)
		initial_index = raw_data2.iloc[(raw_data2['wavelength (nm)']-min_wavelength).abs().argsort()].index[0]
		raw_data2 = raw_data2.iloc[initial_index:].reset_index(drop=True)
		final_index = raw_data2.iloc[(raw_data2['wavelength (nm)']-max_wavelength).abs().argsort()].index[0]
		raw_data2 = raw_data2.iloc[:final_index].reset_index(drop=True)
		absorption_data = raw_data2
		absorption_data_normalized = absorption_data['absorption cross-section (m^2)'].values / np.max(absorption_data['absorption cross-section (m^2)'].values)
		absorption_spectrum = np.squeeze(np.array([[absorption_data['wavelength (nm)'].values], [absorption_data_normalized]], dtype=float))
		interpolated_absorption_spectrum = interp1d(absorption_spectrum[0,:], absorption_spectrum[1,:], kind='cubic')
		
		# emission data
		fluorescence_spectrum_datafile =  Path("data") / 'fluorescence_spectrum_R6G_in_EthyleneGlycol_corrected.csv'
		fluorescence_spectrum_datafile = Path(os.path.dirname(os.path.abspath(__file__))) / fluorescence_spectrum_datafile
		raw_data = pd.read_csv(fluorescence_spectrum_datafile)
		initial_index = raw_data.iloc[(raw_data['wavelength (nm)']-min_wavelength).abs().argsort()].index[0]
		raw_data = raw_data.iloc[initial_index:].reset_index(drop=True)
		final_index = raw_data.iloc[(raw_data['wavelength (nm)']-max_wavelength).abs().argsort()].index[0]
		raw_data = raw_data.iloc[:final_index].reset_index(drop=True)
		fluorescence_data = raw_data
		fluorescence_data_normalized = fluorescence_data['fluorescence (arb. units)'].values / np.max(fluorescence_data['fluorescence (arb. units)'].values)
		emission_spectrum = np.squeeze(np.array([[fluorescence_data['wavelength (nm)'].values], [fluorescence_data_normalized]], dtype=float))
		interpolated_emission_spectrum = interp1d(emission_spectrum[0,:], emission_spectrum[1,:], kind='cubic')

		# Uses both datasets
		if np.min(1e9*np.array(lambdas)) < 480 or np.max(1e9*np.array(lambdas)) > 650:
			raise Exception('*** Restrict wavelength to the range between 480 and 650 nm ***')

		temperature = 300
		lamZPL = 545e-9
		n_mol_per_vol= dye_concentration*sc.Avogadro
		peak_Xsectn = 2.45e-20*n_mol_per_vol*sc.c/n
		wpzl = 2*np.pi*sc.c/lamZPL/1e12

		def freq(wl):
			return 2*np.pi*sc.c/wl/1e12
		def single_exp_func(det):
			f_p = 2*np.pi*sc.c/(wpzl+det)*1e-3
			f_m = 2*np.pi*sc.c/(wpzl-det)*1e-3
			return (0.5*interpolated_absorption_spectrum(f_p)) + (0.5*interpolated_emission_spectrum(f_m))
		def Err(det):
			return Erf(det*1e12)
		def single_adjust_func(det):
			return ((1+Err(det))/2.0*single_exp_func(det)) + ((1-Err(det))/2.0*single_exp_func(-1.0*det)*np.exp(sc.h/(2*np.pi*sc.k*temperature)*det*1e12)) 
			
		emission_rates = np.array([single_adjust_func(-1.0*freq(a_l)+wpzl) for a_l in lambdas])*peak_Xsectn
		absorption_rates = np.array([single_adjust_func(freq(a_l)-wpzl) for a_l in lambdas])*peak_Xsectn

		return absorption_rates, emission_rates
[((78, 14, 78, 55), 'pandas.read_csv', 'pd.read_csv', ({(78, 26, 78, 54): 'absorption_spectrum_datafile'}, {}), '(absorption_spectrum_datafile)', True, 'import pandas as pd\n'), ((86, 37, 86, 111), 'scipy.interpolate.interp1d', 'interp1d', (), '', False, 'from scipy.interpolate import interp1d\n'), ((91, 13, 91, 56), 'pandas.read_csv', 'pd.read_csv', ({(91, 25, 91, 55): 'fluorescence_spectrum_datafile'}, {}), '(fluorescence_spectrum_datafile)', True, 'import pandas as pd\n'), ((99, 35, 99, 105), 'scipy.interpolate.interp1d', 'interp1d', (), '', False, 'from scipy.interpolate import interp1d\n'), ((76, 33, 76, 45), 'pathlib.Path', 'Path', ({(76, 38, 76, 44): '"""data"""'}, {}), "('data')", False, 'from pathlib import Path\n'), ((84, 90, 84, 154), 'numpy.max', 'np.max', ({(84, 97, 84, 153): "absorption_data['absorption cross-section (m^2)'].values"}, {}), "(absorption_data['absorption cross-section (m^2)'].values)", True, 'import numpy as np\n'), ((85, 35, 85, 133), 'numpy.array', 'np.array', (), '', True, 'import numpy as np\n'), ((89, 36, 89, 48), 'pathlib.Path', 'Path', ({(89, 41, 89, 47): '"""data"""'}, {}), "('data')", False, 'from pathlib import Path\n'), ((97, 89, 97, 150), 'numpy.max', 'np.max', ({(97, 96, 97, 149): "fluorescence_data['fluorescence (arb. units)'].values"}, {}), "(fluorescence_data['fluorescence (arb. units)'].values)", True, 'import numpy as np\n'), ((98, 33, 98, 135), 'numpy.array', 'np.array', (), '', True, 'import numpy as np\n'), ((118, 10, 118, 23), 'scipy.special.erf', 'Erf', ({(118, 14, 118, 22): '(det * 1000000000000.0)'}, {}), '(det * 1000000000000.0)', True, 'from scipy.special import erf as Erf\n'), ((77, 54, 77, 79), 'os.path.abspath', 'os.path.abspath', ({(77, 70, 77, 78): '__file__'}, {}), '(__file__)', False, 'import os\n'), ((90, 56, 90, 81), 'os.path.abspath', 'os.path.abspath', ({(90, 72, 90, 80): '__file__'}, {}), '(__file__)', False, 'import os\n'), ((120, 96, 120, 144), 'numpy.exp', 'np.exp', ({(120, 103, 120, 143): '(sc.h / (2 * np.pi * sc.k * temperature) * det * 1000000000000.0)'}, {}), '(sc.h / (2 * np.pi * sc.k * temperature) * det * 1000000000000.0)', True, 'import numpy as np\n'), ((102, 16, 102, 33), 'numpy.array', 'np.array', ({(102, 25, 102, 32): 'lambdas'}, {}), '(lambdas)', True, 'import numpy as np\n'), ((102, 55, 102, 72), 'numpy.array', 'np.array', ({(102, 64, 102, 71): 'lambdas'}, {}), '(lambdas)', True, 'import numpy as np\n')]