In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import scipy.optimize as opt
from scipy.integrate import quad, dblquad
from IPython.html.widgets import interact, interactive, fixed
from IPython.display import display
In [2]:
# 1D ~ 3.3e-30 C.m
DebyeToSI = (1/(299792458))*1e-21
kb = 1.3806488e-23
In [3]:
thetamin = 17.8 * (np.pi/180)
thetamax = 33.4 * (np.pi/180)
In [4]:
df = pd.read_excel('ShenData.xlsx', sheetname=None)
In [5]:
del df['Intro']
In [6]:
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']
In [7]:
df['40C']
Out[7]:
In [8]:
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])
In [9]:
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);
In [32]:
#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
In [33]:
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
In [34]:
compute_psi(E,p0k,alpha)
Out[34]:
In [39]:
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);
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.
Now that compute_psi is functioning, the goal is to reverse engineer the same phenomonological parameters (alpha, p0) that were used above.
In [46]:
guess = (p0k+100,alpha+100)
In [47]:
theta_best, theta_cov = opt.curve_fit(compute_psi, E_array, PSI, guess)
In [48]:
E,alpha,p0k
Out[48]:
In [18]:
E*1e-4
Out[18]:
In [19]:
alpha*1e9
Out[19]:
In [20]:
p0k*1e6
Out[20]: