In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import patches
from sympy import *
from scipy.optimize import fsolve, fmin, root
from collections import Counter, defaultdict
init_printing(use_unicode=True)




In [2]:

r, c, u, v, E = symbols('r c u v E')



## Simple Case

$f = \dfrac{1}{2}rc^{2} - Ec$ with $r > 0$

$c = \dfrac{E}{r}$



In [2]:

def fsimple(r,c,E):
return (1/2)*r*c**2-E*c

def csimple(E,r):
return E/r




In [3]:

c = np.linspace(-5,5,1001)




In [4]:

plt.figure(figsize=(11,8))
plt.gca().spines['left'].set_position('zero')
plt.gca().spines['bottom'].set_position('zero')
plt.gca().spines['right'].set_color('none')
plt.gca().spines['top'].set_color('none')
plt.gca().xaxis.set_ticks_position('bottom')
plt.gca().yaxis.set_ticks_position('left')
plt.xlim(-5,5), plt.ylim(-5,15)
plt.xlabel('c', fontsize=14)
plt.ylabel('f', rotation='horizontal',verticalalignment='center',  fontsize=14)
plt.gca().yaxis.set_label_coords(0.53,1)
plt.gca().xaxis.set_label_coords(1.03,0.22)
yticks1 = plt.gca().yaxis.get_major_ticks()
yticks1[1].label1.set_visible(False)
plt.plot(c,fsimple(0.5,c,1),label='r = 0.5')
plt.plot(c,fsimple(1,c,1),label='r = 1')
plt.plot(c,fsimple(2,c,1),label='r = 2')
plt.legend(loc='best');







In [5]:

E = np.linspace(0,5)




In [6]:

plt.figure(figsize=(11,8))
plt.plot(E,csimple(E,0.5),label='r = 0.5')
plt.plot(E,csimple(E,1),label='r = 1')
plt.plot(E,csimple(E,2),label='r = 2')
plt.plot(E,csimple(E,4),label='r = 4')
plt.plot(E,csimple(E,8),label='r = 8')
plt.plot(E,csimple(E,16),label='r = 16')
plt.gca().spines['right'].set_color('none')
plt.gca().spines['top'].set_color('none')
plt.gca().xaxis.set_ticks_position('bottom')
plt.gca().yaxis.set_ticks_position('left')
plt.xlabel('E',fontsize=18)
plt.legend(loc='best');






### Not so simple case

$f = \dfrac{1}{2}rc^{2} + \dfrac{1}{4}|u|c^{4} - Ec$



In [7]:

def f2(r,c,u,E):
return (1/2)*r*c**2+(1/4)*u*c**4-E*c




In [8]:

plt.figure(figsize=(11,8))
plt.xlabel('c',fontsize=18)
plt.xlim(-3,3), plt.ylim(-5,15)
plt.gca().spines['left'].set_position('zero')
plt.gca().spines['bottom'].set_position('zero')
plt.gca().spines['right'].set_color('none')
plt.gca().spines['top'].set_color('none')
plt.gca().xaxis.set_ticks_position('bottom')
plt.gca().yaxis.set_ticks_position('left')
plt.gca().yaxis.set_label_coords(0.53,1)
plt.gca().xaxis.set_label_coords(1.03,0.27)
yticks2 = plt.gca().yaxis.get_major_ticks()
yticks2[1].label1.set_visible(False)
plt.plot(c,f2(1,c,1,1),label='r = 1')
plt.plot(c,f2(5,c,1,1),label='r = 5')
plt.legend(loc='best');







In [9]:

return ((np.sqrt(12*r**3*u**3+81*u**4*E**2)+9*u**2*E)/2)**(1/3)/(3**(2/3)*u) - (r*(2/3)**(1/3)/((np.sqrt(12*r**3*u**3+81*u**4*E**2)+9*u**2*E))**(1/3))




In [10]:

E2 = np.linspace(-1,5)




In [11]:

plt.figure(figsize=(11,8))
plt.xlim(0,5),plt.ylim(0,2)
plt.xlabel('E',fontsize=18)
plt.legend(loc='best');







In [12]:

def Eofc1(r,c,u):
return r*c+u*c**3




In [13]:

plt.figure(figsize=(11,8))
plt.xlim(-1,5),plt.ylim(0,2)
plt.gca().spines['left'].set_position('zero')
plt.gca().spines['bottom'].set_position('zero')
plt.gca().spines['right'].set_color('none')
plt.gca().spines['top'].set_color('none')
plt.gca().xaxis.set_ticks_position('bottom')
plt.gca().yaxis.set_ticks_position('left')
plt.xlim(-1,5), plt.ylim(0,2)
plt.xlabel('E', fontsize=18)
plt.ylabel('c', rotation='horizontal',verticalalignment='center',  fontsize=18)
plt.gca().yaxis.set_label_coords(0.2,1)
plt.gca().xaxis.set_label_coords(1.03,0.03)
yticks3 = plt.gca().yaxis.get_major_ticks()
yticks3[0].label1.set_visible(False)
plt.plot(Eofc1(-1,c,1),c,label='r = -1')
plt.plot(Eofc1(1,c,1),c,label='r = 1')
plt.plot(Eofc1(5,c,1),c,label='r = 5')
plt.legend(loc='best');






### Even less simple

$f = \dfrac{1}{2}rc^{2} - \dfrac{1}{4}|u|c^{4} + \dfrac{1}{6}|v|c^{6} - Ec$



In [14]:

def f3(c,r,u,v,E):
return (1/2)*r*c**2-(1/4)*abs(u)*c**4+(1/6)*abs(v)*c**6-E*c

def df3(c,r,u,v,E):
return r*c+u*c**3+v*c**5-E




In [16]:

f = (1/2)*r*c**2+(1/4)*u*c**4+(1/6)*v*c**6-E*c
f




---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-16-fd797fb4db5f> in <module>()
----> 1 f = (1/2)*r*c**2+(1/4)*u*c**4+(1/6)*v*c**6-E*c
2 f

ValueError: operands could not be broadcast together with shapes (50,) (1001,)




In [ ]:

E_c = expand(solve(f.diff(c),E)[0])
E_c




In [15]:

plt.plot(df3(np.linspace(0,1),0.3,-1,1,0),np.linspace(0,1))
# plt.scatter(Eu(0.3,-1,1),0.35608274213,c='g',label='$E_U,\ c_L$')
# plt.scatter(El(0.3,-1,1),0.687899033839,c='r',label='$E_L,\ c_U$')
plt.xlim(0,0.1),plt.ylim(0,1)
plt.legend(fontsize=16);




/opt/conda/lib/python3.5/site-packages/matplotlib/axes/_axes.py:531: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots.
warnings.warn("No labelled objects found. "




In [58]:

simplify(c_U.subs([(r,0.3),(u,-1),(v,1)])),simplify(c_L.subs([(r,0.3),(u,-1),(v,1)]))




Out[58]:

$$\left ( 0.217532774716107 \sqrt{10}, \quad 0.112603250061049 \sqrt{10}\right )$$




In [59]:

0.217532774716107*np.sqrt(10),0.112603250061049*np.sqrt(10)




Out[59]:

$$\left ( 0.687899033839, \quad 0.35608274213\right )$$




In [7]:

solve(simplify(E_c.diff(c)),c)




Out[7]:

$$\left [ - \frac{\sqrt{10}}{10} \sqrt{- \frac{3 u}{v} - \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}}, \quad \frac{\sqrt{10}}{10} \sqrt{- \frac{3 u}{v} - \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}}, \quad - \frac{\sqrt{10}}{10} \sqrt{- \frac{3 u}{v} + \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}}, \quad \frac{\sqrt{10}}{10} \sqrt{- \frac{3 u}{v} + \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}}\right ]$$




In [57]:

c_L = solve(simplify(E_c.diff(c)),c)[1]
c_U = solve(simplify(E_c.diff(c)),c)[3]




In [7]:

c_L,c_U




Out[7]:

$$\left ( \frac{\sqrt{10}}{10} \sqrt{- \frac{3 u}{v} - \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}}, \quad \frac{\sqrt{10}}{10} \sqrt{- \frac{3 u}{v} + \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}}\right )$$




In [40]:

E_c.subs(c,c_U)




Out[40]:

$$\frac{\sqrt{10} r}{10} \sqrt{- \frac{3 u}{v} + \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}} + \frac{\sqrt{10} u}{100} \left(- \frac{3 u}{v} + \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}\right)^{\frac{3}{2}} + \frac{\sqrt{10} v}{1000} \left(- \frac{3 u}{v} + \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}\right)^{\frac{5}{2}}$$




In [43]:

E_c.subs(c,c_L)




Out[43]:

$$\frac{\sqrt{10} r}{10} \sqrt{- \frac{3 u}{v} - \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}} + \frac{\sqrt{10} u}{100} \left(- \frac{3 u}{v} - \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}\right)^{\frac{3}{2}} + \frac{\sqrt{10} v}{1000} \left(- \frac{3 u}{v} - \frac{1}{v} \sqrt{- 20 r v + 9 u^{2}}\right)^{\frac{5}{2}}$$




In [13]:

E_U = simplify(E_c.subs(c,c_L))
E_L = simplify(E_c.subs(c,c_U))
E_U,E_L




Out[13]:

$$\left ( \frac{\sqrt{10}}{250 v} \sqrt{\frac{1}{v} \left(- 3 u - \sqrt{- 20 r v + 9 u^{2}}\right)} \left(20 r v - 3 u^{2} - u \sqrt{- 20 r v + 9 u^{2}}\right), \quad \frac{\sqrt{10}}{250 v} \sqrt{\frac{1}{v} \left(- 3 u + \sqrt{- 20 r v + 9 u^{2}}\right)} \left(20 r v - 3 u^{2} + u \sqrt{- 20 r v + 9 u^{2}}\right)\right )$$




In [16]:

r = np.linspace(0,0.5)




In [17]:

def E_U(r,u,v):
return (np.sqrt(10*v*(-3*u-np.sqrt(9*u**2-20*r*v)))*(20*r*v-3*u**2-u*np.sqrt(9*u**2-20*r*v))/(250*v**2))

def E_L(r,u,v):
return (np.sqrt(10*v*(-3*u+np.sqrt(9*u**2-20*r*v)))*(20*r*v-3*u**2+u*np.sqrt(9*u**2-20*r*v))/(250*v**2))




In [18]:

plt.figure(figsize=(11,8))
plt.xlabel('c',fontsize=18)
plt.xlim(-2,2), plt.ylim(-5,5)
plt.gca().spines['left'].set_position('zero')
plt.gca().spines['bottom'].set_position('zero')
plt.gca().spines['right'].set_color('none')
plt.gca().spines['top'].set_color('none')
plt.gca().xaxis.set_ticks_position('bottom')
plt.gca().yaxis.set_ticks_position('left')
plt.gca().yaxis.set_label_coords(0.53,0.97)
plt.gca().xaxis.set_label_coords(1.04,0.52)
yticks4 = plt.gca().yaxis.get_major_ticks()
yticks4[3].label1.set_visible(False)
plt.plot(c,f3(c,1,1,1,1),label='r = 1')
plt.plot(c,f3(c,3,1,1,1),label='r = 3')
plt.plot(c,f3(c,5,1,1,1),label='r = 5')
plt.legend(loc='best');







In [19]:

def Eofc2(r,c,u,v):
return r*c-u*c**3+v*c**5




In [44]:

plt.figure(figsize=(11,8))
plt.gca().spines['left'].set_position('zero')
plt.gca().spines['bottom'].set_position('zero')
plt.gca().spines['right'].set_color('none')
plt.gca().spines['top'].set_color('none')
plt.gca().xaxis.set_ticks_position('bottom')
plt.gca().yaxis.set_ticks_position('left')
plt.xlim(0,0.17), plt.ylim(0,1)
plt.yticks([0,0.2,0.4,0.6,0.8,1.0],('0','0.2','0.4','0.6','0.8','1.0'),fontsize=18)
plt.xticks([0,0.05,0.1,0.15],('0','0.05','0.1','0.15'),fontsize=18)
plt.xlabel('E-Field (E)', fontsize=18)
plt.ylabel('Tilt ($\mathregular{\\theta}$)', rotation='horizontal',verticalalignment='center', fontsize=18)
plt.title('Figure 2',y=-0.15,fontsize=20)
plt.gca().text(0.15, 0.85, 'High Tilt', fontsize=18)
plt.text(0.01, 0.79, 'Low Temp', fontsize=18,color='blue')
plt.text(0.15, 0.36, 'High Temp', fontsize=18,color='red')
plt.gca().text(0.02, 0.2, 'Low Tilt', fontsize=18)
plt.gca().yaxis.set_label_coords(0.07,1)
plt.gca().xaxis.set_label_coords(1.1,0.02)
yticks5 = plt.gca().yaxis.get_major_ticks()
yticks5[0].label1.set_visible(False)
# plt.plot(Eofc2(0.15,c,1,1),c,label='r(T) = 0.15')
# plt.plot(Eofc2(0.2,c,1,1),c,label='r(T) = 0.2')
plt.plot(Eofc2(0.3,c,1,1),c,label='r(T) = 0.3')
# plt.plot([0.13145,0.13145],[0,1],'--k')
# plt.plot(Eofc2(0.4,c,1,1),c,label='r(T) = 0.4')
plt.plot(Eofc2(0.45,c,1,1),c,label='r(T) = 0.45')
plt.plot(Eofc2(0.5,c,1,1),c,label='r(T) = 0.5')
plt.legend(loc='lower right',fontsize=14);







In [45]:

s = np.linspace(0,0.5,101)




In [46]:

r3 = 3/16
Ec = np.sqrt(3/10)*(6/25)
r_T, E_T = [],[]
u=-1
v=1
n = 100
for i in np.linspace(0,Ec,n): #E
for j in np.linspace(r3,9/20,n): #r
if (E_L(j,u,v) < i) and (i < E_U(j,u,v)):
c1 = fsolve(df3,np.array([0,1]),args=(j,u,v,i))

f1 = f3(c1[0],j,u,v,i)
f2 = f3(c1[1],j,u,v,i)

if f1 <= f2:
r_T.append(j)
E_T.append(i)
r3 = j
break




In [48]:

plt.figure(figsize=(11,8))
plt.plot(s,E_U(s,u,v),label='$\mathregular{E_{Upper}}$')
plt.plot(s,E_L(s,u,v),label='$\mathregular{E_{Lower}}$')
plt.fill_between(s,E_U(s,u,v),E_L(s,u,v), color='#D98880',alpha='0.5')
plt.plot(r_T,E_T,'--r',label='$\mathregular{E_{T}}$')
plt.scatter(0.45,E_L(s,-1,1)[90],c='k')
plt.plot([0,0.45],[E_L(s,-1,1)[90],E_L(s,-1,1)[90]],'--k')
plt.plot([0.45,0.45],[0,E_L(s,-1,1)[90]],'--k')
plt.gca().spines['left'].set_position('zero')
plt.gca().spines['bottom'].set_position('zero')
plt.gca().spines['right'].set_color('none')
plt.gca().spines['top'].set_color('none')
plt.gca().xaxis.set_ticks_position('bottom')
plt.gca().yaxis.set_ticks_position('left')
yticks6 = plt.gca().yaxis.get_major_ticks()
plt.title('Figure 3',y=-0.15,fontsize=20)
yticks6[0].label1.set_visible(False)
plt.gca().yaxis.set_label_coords(0.13,0.97)
plt.gca().xaxis.set_label_coords(1.05,0.02)
plt.xticks([0,0.1875,0.25,0.45,0.5],('$\mathregular{r_{2}}$','$\mathregular{r_{T}}$','$\mathregular{r_{1}}$','$\mathregular{r_{C}}$','0.5'),fontsize=18)
plt.yticks([0,0.05,0.10,0.13145,0.15,0.20],('0','0.05','0.10','$\mathregular{E_{C}}$','0.15','0.20'),fontsize=18)
plt.gca().annotate('phase boundary', xy=(0.37,0.09), xytext=(0.2, 0.12),arrowprops=dict(arrowstyle="->",facecolor='black'),fontsize=16)
plt.gca().text(0.1, 0.06, 'High Tilt', fontsize=18)
plt.gca().text(0.35, 0.04, 'Low Tilt', fontsize=18)
plt.gca().text(0.47, 0.142, 'B', fontsize=16)
plt.gca().text(0.1, 0.1, '$\mathregular{u_{\\theta} = -1}$', fontsize=14)
plt.gca().text(0.1, 0.09, '$\mathregular{v_{\\theta} = 1}$', fontsize=14)
plt.gca().annotate("A", xy=(0.385, 0.06), xytext=(0.245, 0.075), arrowprops=dict(arrowstyle="->",connectionstyle='arc3'),fontsize=16)
plt.gca().annotate("", xy=(0.4125,0.131), xytext=(0.4165,0.135),arrowprops=dict(arrowstyle="->",connectionstyle='arc3'))
plt.gca().annotate("", xy=(0.465,0.122), xytext=(0.48,0.126),arrowprops=dict(arrowstyle="->",connectionstyle='arc3'))
e1 = patches.Arc((0.45, 0.13145), 0.07, 0.02, angle=0, theta1=-60, theta2=165, linewidth=1, fill=False, zorder=2)
plt.xlabel('r(T)',fontsize=18)
plt.ylabel('E',fontsize=18,rotation='horizontal')
plt.ylim(0,0.15),plt.xlim(-0.05,0.6)
plt.legend(loc='lower right',fontsize=16);




/opt/conda/lib/python3.5/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in sqrt
from ipykernel import kernelapp as app
/opt/conda/lib/python3.5/site-packages/ipykernel/__main__.py:5: RuntimeWarning: invalid value encountered in sqrt



### $E_{C} = 0.13145$



In [23]:

Elower(0.45,-1,1),Eupper(0.45,-1,1)




Out[23]:

(0.13145341380123987, 0.13145341380123987)




In [24]:

for i in range(s.shape[0]):
if Elower(s,-1,1)[i]==Eupper(s,-1,1)[i]:
print(i,s[i])




90 0.45

/usr/local/lib/python3.4/dist-packages/IPython/kernel/__main__.py:9: RuntimeWarning: invalid value encountered in sqrt
/usr/local/lib/python3.4/dist-packages/IPython/kernel/__main__.py:4: RuntimeWarning: invalid value encountered in sqrt




In [25]:

for i in range(s.shape[0]):
if Eupper(s,-1,1)[i]==0:
print(i,s[i])




50 0.25

/usr/local/lib/python3.4/dist-packages/IPython/kernel/__main__.py:4: RuntimeWarning: invalid value encountered in sqrt