In [18]:
import numpy as np
from scipy.integrate import quad, dblquad
%matplotlib inline
import matplotlib.pyplot as plt
from scipy.optimize import minimize
In [19]:
thetamin = 25.6*np.pi/180
thetamax = 33.7*np.pi/180
t = 1*10**-6 #Cell Thickness
In [20]:
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])
In [21]:
voltages = np.array([1.826,3.5652,5.3995,7.2368,9.0761,10.8711,12.7109,14.5508,16.3461,18.1414,19.9816,21.822,23.6174,25.4577,27.253,29.0935,30.889,32.7924,34.5699,35.8716])
measured_psi1 = 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]])
In [22]:
#This Block just converts units
fields = np.array([entry/t for entry in voltages])
KC = 273.15
tempsK = np.array([entry+KC for entry in tempsC]) #Celsius to Kelvin
In [23]:
# measured_psi1 = np.array([[11,20.5,25.5,27.5,29,30,30.5,31,31.25,31.5,31.75,32,32.25,32.5,32.75,33,33.25,33.5,33.75,34]])
# measured_psi2 = np.array([[7.6, 11.5, 22.3, 24.7, 27.8, 29.4, 30.1, 30.7, 31.2, 31.6, 31.9, 32.2, 32.4, 32.6, 32.7, 32.8, 32.9, 32.9, 33.0, 33.1]])
# measured_psi3 = np.array([[4.7, 7.3, 15.5, 18.1, 22.7, 25.9, 27.5, 28.6, 29.6, 30.3, 30.8, 31.2, 31.5, 31.8, 32.0, 32.1, 32.3, 32.4, 32.5, 32.6]])
# measured_psi4 = np.array([[3.5, 5.4, 11.5, 13.8, 18.1, 21.9, 24.1, 25.9, 27.5, 28.7, 29.5, 30.1,30.5, 31.0, 31.3, 31.5, 31.7, 31.9, 32.0, 32.2]])
# measured_psi5 = np.array([[2.5, 3.7, 8.0, 9.6, 12.9, 16.3, 18.7, 20.9, 23.4, 25.3, 26.8, 27.9, 28.5, 29.4, 29.8, 30.2, 30.6, 30.8, 31.1, 31.3]])
# measured_psi6 = np.array([[1.9, 2.9, 6.1, 7.3, 9.8, 12.6, 14.7, 16.8, 19.4, 21.7, 23.6, 25.2, 26.1, 27.4, 28.0, 28.6, 29.2, 29.5, 29.9, 30.3]])
# measured_psi7 = np.array([[1.5, 2.3, 4.7, 5.6, 7.5, 9.6, 11.2, 12.9, 15.2, 17.5, 19.6, 21.4, 22.7, 24.4, 25.37, 26.1, 27.02, 27.5, 28.0, 28.6]])
In [24]:
# AllPsi = np.concatenate((measured_psi1,measured_psi2,measured_psi3,measured_psi4,measured_psi5,measured_psi6,measured_psi7),axis=0)
In [25]:
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 [26]:
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 [27]:
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 [28]:
def compute_psi(T,p0k,alpha,E,thetamin,thetamax):
"""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 = (1/2)*np.arctan(avg_numerator / (avg_denominator)) * (180 /(np.pi)) #Converting to degrees from radians and divide by two
return psi
In [29]:
def compute_error(xo,fields,T,thetamin,thetamax,measured_psi):
"""Computes the squared error for a pair of parameters by comparing it to all measured tilt angles
at one temperature.
This will be used with the minimization function, xo is a point that the minimization checks.
Parameters/Conditions
----------
x0:
An array of the form [alpha^13,p0^33].
Returns
-------
Float: Error
"""
alpha = xo[0]/(1e10)
p0 = xo[1]/(1e30)
p0k = p0/1.3806488e-23
computed_psi = np.array([compute_psi(T,p0k,alpha,E,thetamin,thetamax) for E in fields])
Err = computed_psi - measured_psi
ErrSqr = np.array([i**2 for i in Err])
return np.sum(ErrSqr)*1e8 #Scaling the Squared Error up here seems to help with minimization precision.
It might be better to use the minimization function individually for each temperature range. The minimization function returns a minimization object, which gives extra information about the results. The two important entries are fun and x.
fun is the scalar value of the function that is being minimized. In our case fun is the squared error.
x is the solution array of the form [alpha^10,p0^30]
The reason it might be better to just minimze the squared error function, instead of using the minimize_func that I wrote below is because the minimize function is very picky about the initial guess. Also the minimization function tends to stop when the result of the function is one the order of 10^-3.
Right now everything below this might not work as well as manually guessing and checking. The idea for this section was to automate that process and just return our entire solution arrays at the end of the notebook.
In [30]:
def minimize_func(guess,fields,T,thetamin,thetamax,measured_psi,bnds):
"""A utility function that is will help me construct alpha and p0 arrays later.
Uses the imported minimize function and compute_error to best fit our parameters
at a temperature.
Parameters/Conditions
----------
guess:
The initial guess for minimize().
Returns
-------
Array: [alpha,p0]
"""
results = minimize(compute_error,guess,args=(fields,T,thetamin,thetamax,measured_psi),method = 'SLSQP',bounds = bnds)
xres = np.array(dict(results.items())['x'])
"""Minimize returns a special minimization object. That is similar to a dictionary but not quite.
xres is grabbing just the x result of the minimization object, which is the [alpha,p0] array that
we care about"""
alpha_results = xres[0]
p0_results = xres[1]
return np.array([alpha_results,p0_results])
In [31]:
guess = (2575,2168)
bnds = ((1000,2600),(200,2400))
In [32]:
results = minimize(compute_error,guess,args=(fields,tempsK[0],thetamin,thetamax,measured_psi1),method = 'TNC',bounds = bnds)
results
Out[32]:
In [33]:
res = np.array(dict(results.items())['x'])
alpha = res[0]
p0 = res[1]
alpha = alpha*1e-4
p0 = p0/3.33564
print("alpha micro: " + str(alpha))
print('p0 debye: ' + str(p0))
In [37]:
#Minimization claims that it did not succeed. But the results were pretty good. I think it believes that it did not succeed because I have the squared error scaled up very high.
In [19]:
def solution(initial_guess,fields,tempsK,thetamin,thetamax,AllPsi,initial_bnds):
"""Constructs Alpha and p0 arrays where each entry is the value of alpha,p0 at the corresponding temperature in
tempsK. Initial guess and initial bounds are changed each iteration of the loop to the previous values of alpha and p0.
Alpha and p0 decrease so this helps to cut down on the range.
Parameters/Conditions
----------
initial_guess:
The initial guess for minimize().
initial_bnds:
The initial bounds for minimize().
Returns
-------
Array,Array: Alpha Array in micro meters, p0 Array in debye
"""
alpha = np.array([])
p0 = np.array([])
guess = initial_guess
bnds = initial_bnds
for i in range(len(tempsK)):
res = minimize_func(guess,fields,tempsK[i],thetamin,thetamax,AllPsi[i],bnds)
alpha = np.append(alpha,res[0])
p0 = np.append(p0,res[1])
guess = (res[0]-10,res[1]-10)
bnds = ((initial_bnds[0][0],res[0]),(initial_bnds[1][0],res[1]))
alpha = alpha*1e-4
p0 = p0/(3.33564)
return alpha,p0
In [20]:
initial_guess = (2575,2168)
initial_bnds = ((1000,2600),(200,2300))
In [21]:
alpha_micro,p0Debye = solution(initial_guess,fields,tempsK,thetamin,thetamax,AllPsi,initial_bnds)