<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 TO DO
Step1:  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.
Material Characteristics
Step2:  Data
Step4:  Langevin-Debye Model
$$ U(\theta,\phi) = -\rho_0E\
Step5:  Second, Calculate the Tilt Angle $\psi$
$$ numerator() \
Step6:  $$ denominator()\
Step9:  $$ tan(2\psi) =  \frac{\int_{\theta_{min}}^{\theta_{max}} \int_0^{2\pi} {sin\
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from scipy.integrate import quad, dblquad
%matplotlib inline
import matplotlib.pyplot as plt
import scipy.optimize as opt
Explanation: 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
End of explanation
thetamin = 25.6*np.pi/180
thetamax = 33.7*np.pi/180
t = 1*10**-6 #Cell Thickness
Explanation: 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.
Material Characteristics
End of explanation
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
#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])
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]
Explanation: Data
End of explanation
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)
Explanation: Langevin-Debye Model
$$ U(\theta,\phi) = -\rho_0E\:sin\:\theta\:cos\:\phi\:(1+\alpha E\:cos\:\phi) $$
First, Calculate the Boltzmann Factor and the Partition Function
$$ {Boltz() returns:}\:\: e^{\frac{-U}{k_bT}}\:sin\:{\theta}\ $$
End of explanation
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
Explanation: Second, Calculate the Tilt Angle $\psi$
$$ numerator() \:returns: {sin\:{2\theta}\:cos\:{\phi}}\:e^{\frac{-U}{k_bT}}\:sin\:{\theta} $$
End of explanation
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
Explanation: $$ denominator()\: returns: {({cos}^2{\theta} - {sin}^2{\theta}\:{cos}^2{\phi}})\:e^{\frac{-U}{k_bT}}\:sin\:{\theta} $$
End of explanation
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
compute_psi(tdata[0],p0k_array[0]*1e7,alpha_array[0]*1e10)
PSIdata[0]
guess = [p0k_array[0]*1e7,alpha_array[0]*1e10]
guess
theta_best, theta_cov = opt.curve_fit(compute_psi, Edata, PSIdata,guess,absolute_sigma=True)
Explanation: $$ 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} $$
End of explanation
<END_TASK>