TO DO:
Need to be able to scatter plot measured values of Psi on top of the current Psi plot.
Alpha and rho LaTeX not working in plots.
Legend needs to be move in the Psi plot.
Consider also applying the scatter of the data and curve fit style of plotting for the parameters
Plots in general need to look better
LaTeX in all of the equations
In [38]:
import numpy as np
from scipy.integrate import quad, dblquad
%matplotlib inline
import matplotlib.pyplot as plt
import scipy.optimize as opt
This notebook calculates and plots the theoretical tilt angles. It will also plot the alpha and p0 factors vs temperature that are given in the cell below this.
In [39]:
thetamin = 25.6*np.pi/180
thetamax = 33.7*np.pi/180
t = 1*10**-6 #Cell Thickness
In [40]:
tempsC = np.array([26, 27, 29, 31, 33, 35, 37])
voltages = np.array([2,3,6,7,9,11,12.5,14,16,18,20,22,23.5,26,27.5,29,31,32.5,34,36])
alpha_micro = np.array([.2575,.2475,.2275,.209,.189,.176,.15])
p0Debye = np.array([650,475,300,225,160,125,100]) #Temperature Increases to the right
In [41]:
#This Block just converts units
fields = np.array([entry/t for entry in voltages])
debye = 3.33564e-30
p0_array = np.array([entry*debye for entry in p0Debye]) #debye units to SI units
k = 1.3806488e-23
p0k_array = np.array([entry/k for entry in p0_array]) #p0k is used because it helps with the integration
KC = 273.15
tempsK = np.array([entry+KC for entry in tempsC]) #Celsius to Kelvin
alpha_array = np.array([entry*1e-6 for entry in alpha_micro])
In [88]:
PSIdata = np.array([11.4056,20.4615,25.4056,27.9021,29.028,29.6154,30.2517,30.8392,31.1329,31.5245,31.8671,32.014,32.3077,32.5034,32.7972,32.9929,33.1399,33.3357,33.4336,33.6783])
Edata = fields
T = tempsK[0]
In [43]:
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+alpha*E*np.cos(phi)))*np.sin(theta)
In [45]:
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
In [46]:
def denominator(theta,phi,T,p0k,alpha,E):
boltz = Boltz(theta,phi,T,p0k,alpha,E)
return ((np.cos(theta)**2) - ((np.sin(theta)**2) * (np.cos(phi)**2)))*boltz
In [83]:
def compute_psi(E,p0k,alpha):
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+alpha*E*np.cos(phi)))*np.sin(theta)
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)**2) * (np.cos(phi)**2)))*boltz
"""Computes the tilt angle(psi) by use of our tan(2psi) equation
Returns
-------
Float:
The statistical tilt angle with conditions T,p0k,alpha,E
"""
avg_numerator, avg_numerator_error = dblquad(numerator, 0, 2*np.pi, lambda theta: thetamin, lambda theta: thetamax,args=(T,p0k,alpha,E))
avg_denominator, avg_denominator_error = dblquad(denominator, 0, 2*np.pi, lambda theta: thetamin, lambda theta: thetamax,args=(T,p0k,alpha,E))
psi = np.arctan(avg_numerator / (avg_denominator)) * (180 /(2*np.pi)) #Converting to degrees from radians and divide by two
return psi
In [84]:
compute_psi(tdata[0],p0k_array[0]*1e7,alpha_array[0]*1e10)
Out[84]:
In [91]:
PSIdata[0]
Out[91]:
In [85]:
guess = [p0k_array[0]*1e7,alpha_array[0]*1e10]
In [86]:
guess
Out[86]:
In [90]:
theta_best, theta_cov = opt.curve_fit(compute_psi, Edata, PSIdata,guess,absolute_sigma=True)
In [ ]: