# 1D ~ 3.3e-30 C.m
DebyeToSI = (1/(299792458))*1e-21
kb = 1.3806488e-23

"""
Explanation: Constants
End of explanation
"""


thetamin = 17.8 * (np.pi/180)
thetamax = 33.4 * (np.pi/180)

df = pd.read_excel('ShenData.xlsx', sheetname=None)

del df['Intro']

for key in df:
    del df[key]['Vp']
    del df[key]['θ1']
    del df[key]['θ2']
    del df[key]['Extinction +']
    del df[key]['Extinction -']
    del df[key]['Retardation']
    del df[key]['Birefringence']

df['40C']

Fields = np.array([df[key]['Field (V/μm)'] for key in df])
Angles = np.array([df[key]['Tilt Angle'] for key in df])
TempsC = np.array([df[key]['Temperature (Celsius)'][0] for key in df])
TempsK = np.array([entry + 273.15 for entry in TempsC])

plt.figure(figsize=(8,6));
colors = np.array(['red','green','blue','orange','brown','pink', 'purple','black','cyan','magenta','yellow','gray'])
for i in range(len(Fields)):
    plt.scatter(Fields[i],Angles[i],c=colors[i],label= str(int(TempsC[i])) + '$^\circ$C',marker = 'o',s=50, edgecolor = colors[i]);
plt.xlabel('Field (V/$\mu$m)',size=15);
plt.ylabel('Tilt Angle',size=15);
plt.title('Tilt Angle for Shen Data')
plt.xlim(0,70);
plt.ylim(0,40);
plt.legend(loc=4);

"""
Explanation: 1. Import Excel Data
End of explanation
"""


#Test values at 40C
E_array = np.array(Fields[-1])
T = TempsK[-1]
PSI = np.array(Angles[-1])

#Phenomenological points from Shen paper at 40C
p0 = 619 #Debye
alpha = .338 # Micron/Volt


#Convert to SI
p0 = p0 * DebyeToSI
alpha = .338 * 1e-6
E = E_array[1]*1e6


#Combine p0 and kb to a single more reasonable number
p0k = p0/kb


#Scaling try: Puts E, alpha, and p0k all on the same order of magnitude ~200.
E = E*1e-4
alpha = alpha*1e9
p0k = p0k*1e6
E_array = E_array*1e2

def compute_psi(E,p0k,alpha):
    """Computes the tilt angle(psi) via the tan(2psi) equation at one temperature.
    Returns
    -------
    Float:
    The statistical tilt angle with conditions T,p0k,alpha,E"""
    
    def Boltz(theta,phi,T,p0k,alpha,E):
        """Compute the integrand for the Boltzmann factor.
        Returns
        -------
        A function of theta,phi,T,p0k,alpha,E to be used within dblquad
        """
        return np.exp((1/T)*p0k*E*np.sin(theta)*np.cos(phi)*((1*1e-2)+(alpha*E*np.cos(phi)*1e-7)))*np.sin(theta)
        #Factors: 1e-2 and 1-e7 and to cancel out the rescaling, of E,alpha, and p0k, done in the above cell
    
    
    def numerator(theta,phi,T,p0k,alpha,E):
        boltz = Boltz(theta,phi,T,p0k,alpha,E)
        return np.sin(2*theta)*np.cos(phi)*boltz

    def denominator(theta,phi,T,p0k,alpha,E):
        boltz = Boltz(theta,phi,T,p0k,alpha,E)
        return (np.cos(theta)**2 - ((np.sin(theta)*np.cos(phi))**2))*boltz
    
    avg_numerator, avg_numerator_error = dblquad(numerator, 0, 2*np.pi, lambda theta1: thetamin, lambda theta2: thetamax,args=(T,p0k,alpha,E))
    
    avg_denominator, avg_denominator_error = dblquad(denominator, 0, 2*np.pi, lambda theta3: thetamin, lambda theta4: thetamax,args=(T,p0k,alpha,E))
    
    psi = np.arctan(avg_numerator / avg_denominator) * (180 /(2*np.pi)) #Convert to degrees from radians and divide by two
    return psi

compute_psi(E,p0k,alpha)

plt.plot(E_array*1e-2,[compute_psi(entry,p0k,alpha) for entry in E_array],label = "Compute Output")
plt.scatter(E_array*1e-2,PSI,label="Shen's Data",color = 'red');
#Factor: 1e-2 is added to cancel out the rescaling of E_array
plt.xlabel('Field (V/$\mu$m)',size=15);
plt.ylabel('Tilt Angle',size=15);
plt.title('Tilt Angle Comparison for ' + str(T-273.15) + "$^\circ$C");
plt.xlim(0,70);
plt.ylim(0,40);
plt.legend(loc=4);

"""
Explanation: 2. Curve Fitting
Langevin-Debye Model
$$ U(\theta, \phi, E, \rho_0) = -\rho_0E\:sin\:\theta\:cos\:\phi\:(1+\alpha E\:cos\:\phi) $$
$$ tan(2\psi) =  \frac{\langle sin\:2\theta\:cos\:{\phi}\rangle}{\langle{cos}^2{\theta} - {sin}^2{\theta}\:{cos}^2{\phi}\rangle} $$
$$ \langle X \rangle = \int_{\theta_{min}}^{\theta_{max}} \int_0^{2\pi} X(\theta,\phi) f(\theta, \phi) sin(\theta) d\theta d\phi   $$
$$ f(\theta,\phi) = \frac{e^{\frac{-U}{k_b T}}}{\int_{\theta_{min}}^{\theta_{max}} \int_0^{2\pi} e^{\frac{-U}{k_b T}} \:sin\:{\theta}\: d\theta d\phi}  $$
$$ tan(2\psi) = \frac{\int_{\theta_{min}}^{\theta_{max}} \int_0^{2\pi} {sin\:{2\theta}\:cos\:{\phi}}\:e^{\frac{-U}{k_bT}}\:sin\:{\theta}\: d\theta d\phi}{\int_{\theta_{min}}^{\theta_{max}} \int_0^{2\pi} ({{cos}^2{\theta} - {sin}^2{\theta}\:{cos}^2{\phi}})\:e^{\frac{-U}{k_bT}}\:sin\:{\theta}\: d\theta d\phi}  $$
Computing
End of explanation
"""


guess = (p0k+100,alpha+100)

theta_best, theta_cov = opt.curve_fit(compute_psi, E_array, PSI, guess)

E,alpha,p0k

E*1e-4

alpha*1e9

p0k*1e6

"""
Explanation: Compute_psi is working as intended. When supplied the correct parameters as taken from Shen's paper, Compute_psi outputs the same curve as seen in Shen's paper.
Scipy Curve Fit
Now that compute_psi is functioning, the goal is to reverse engineer the same phenomonological parameters (alpha, p0) that were used above.
End of explanation
"""


