Work on getting $E$ vs $\theta$ plot parameters w/ 3 extrema



In [1]:

    
%matplotlib inline
from sympy import *
import matplotlib.pyplot as plt
import numpy as np
init_printing(use_unicode=True)



In [2]:

    
r, u, v, c, r_c, u_c, v_c, E, p, r_p, u_p, v_p, e, a, b, q, b_0, b_1, b_2, b_3, q_0, q_1, q_2, q_3, q_4, q_5, beta, rho, epsilon, delta, d, K_3, Omega, Lambda, lamda, C, mu, Gamma, tau, nu, xi, P_x, eta, varep, gamma, P_0, theta_0, z, a_0, alpha_0, alpha, T_p, T_c, T = symbols('r u v c r_c u_c v_c E p r_p u_p v_p e a b q b_0 b_1 b_2 b_3 q_0 q_1 q_2 q_3 q_4 q_5 beta rho epsilon delta d K_3 Omega Lambda lamda C mu Gamma tau nu xi_2 P_x eta varepsilon gamma P_0 theta_0 z a_0 alpha_0 alpha T_p T_c T')



In [3]:

    
eptil, atil, btil, ctil, Ctil, gtil, thetatil, Ptil, gprm, thetaprm, Pprm, Tprm, chiprm = symbols('epsilontilde atilde btilde ctilde Ctilde gtilde_{0} thetatilde_{0} Ptilde_{0} gprm thetaprm Pprm Tprm chiprm')

Generalized Landau Model of Ferroelectric Liquid Crystals



In [12]:

    
def g1(a,b,q,xi,P_x,varep,eta,C,Omega):
    return (a*xi**2)/2+(b*xi**4)/4+(q*xi**6)/6+P_x**2/(2*varep)+(eta*P_x**4)/4+C*P_x*xi-(Omega*(P_x*xi)**2)/2



In [13]:

    
g1(a,b,q,xi,P_x,varep,eta,C,Omega)









    Out[13]:





$$C P_{x} \xi_{2} - \frac{\Omega P_{x}^{2}}{2} \xi_{2}^{2} + \frac{P_{x}^{4} \eta}{4} + \frac{P_{x}^{2}}{2 \varepsilon} + \frac{a \xi_{2}^{2}}{2} + \frac{b \xi_{2}^{4}}{4} + \frac{q \xi_{2}^{6}}{6}$$



In [14]:

    
g = (a*xi**2)/2+(b*xi**4)/4+(q*xi**6)/6+P_x**2/(2*epsilon)+(eta*P_x**4)/4+C*P_x*xi-(Omega*(P_x*xi)**2)/2



In [15]:

    
g









    Out[15]:





$$C P_{x} \xi_{2} - \frac{\Omega P_{x}^{2}}{2} \xi_{2}^{2} + \frac{P_{x}^{4} \eta}{4} + \frac{P_{x}^{2}}{2 \epsilon} + \frac{a \xi_{2}^{2}}{2} + \frac{b \xi_{2}^{4}}{4} + \frac{q \xi_{2}^{6}}{6}$$

$f(c,p) = \dfrac{1}{2}r_{c}c^{2}+\dfrac{1}{4}u_{c}c^{4}+\dfrac{1}{6}v_{c}c^{6}+\dfrac{1}{2}r_{p}p^{2}+\dfrac{1}{4}u_{p}p^{4}+\dfrac{1}{6}v_{p}p^{6}-\gamma cp-\dfrac{1}{2}ec^{2}p^{2}-Ep$

'Our' version of the equation

$P_x \neq 0 = p$

$P_y = 0$

$\xi_1 = 0$

$\xi_2 \neq 0 = c$

$r_p = \dfrac{1}{\epsilon}$

$u_p = \eta$

$\gamma = -C$

$e = \dfrac{1}{2}\Omega$

Missing: $p^6$ and $E$ terms



In [16]:

    
g.subs([(epsilon,1/r_p),(eta,u_p),(C,-gamma),(Omega,2*e),(xi,c),(P_x,p)])









    Out[16]:





$$\frac{a c^{2}}{2} + \frac{b c^{4}}{4} + \frac{c^{6} q}{6} - c^{2} e p^{2} - c \gamma p + \frac{p^{4} u_{p}}{4} + \frac{p^{2} r_{p}}{2}$$

Electrically Induced Tilt in Achiral Bent-Core Liquid Crystals



In [4]:

    
fc = a*c**2+b*c**4+q*c**6
fp = alpha*p**2+beta*p**4+gamma*p**6-E*p
fcp = -Omega*(p*c)**2



In [5]:

    
fc









    Out[5]:





$$a c^{2} + b c^{4} + c^{6} q$$



In [6]:

    
fp









    Out[6]:





$$- E p + \alpha p^{2} + \beta p^{4} + \gamma p^{6}$$



In [7]:

    
fcp









    Out[7]:





$$- \Omega c^{2} p^{2}$$

Mapping:

$\dfrac{1}{2}r_c = a = a_{0}(T-T_c),\ a_{0} = 6\times10^4$ $\text{ N/K m}^2$

$\dfrac{1}{4}u_c = b = 10^4 \text{ N/m}^2$

$\dfrac{1}{6}v_c = q = 9\times10^7 \text{ N/m}^2$

$\dfrac{1}{2}r_p = \alpha = \alpha_{0}(T-T_p) = 0.6\times10^9 \text{ J m/C}^2$

$\dfrac{1}{4}u_p = \beta \approx 10^4 \dfrac{ J m^5}{C^4}$

$\dfrac{1}{6}v_p = \gamma = 0$

$e = \Omega \approx 10^{10}-10^{11} \text{ N/C}$



In [8]:

    
# fc = fc.subs([(a_0*(T-T_c),a)])#,6e4),(b,10**4),(q,9e7)])
# fc



In [8]:

    
fp = fp.subs([(gamma,0)])#,(alpha_0*(T-T_p),alpha),(beta,10**4)])
fp









    Out[8]:





$$- E p + \alpha p^{2} + \beta p^{4}$$

$\dfrac{\partial f}{\partial c} = 0$



In [9]:

    
collect((fc+fp+fcp).diff(c),c)









    Out[9]:





$$4 b c^{3} + 6 c^{5} q + c \left(- 2 \Omega p^{2} + 2 a\right)$$

Solve for $p$



In [10]:

    
pc = solve((fc+fp+fcp).diff(c),p)[1]
pc









    Out[10]:





$$\sqrt{\frac{1}{\Omega} \left(a + 2 b c^{2} + 3 c^{4} q\right)}$$

Solve $\dfrac{\partial f}{\partial p} = 0$ for $E$



In [11]:

    
solve((fc+fp+fcp).diff(p),E)[0]









    Out[11]:





$$2 p \left(- \Omega c^{2} + \alpha + 2 \beta p^{2}\right)$$

Sub $p(c)$ into $E(c,p)$



In [12]:

    
(solve((fc+fp+fcp).diff(p),E)[0]).subs(p,pc)









    Out[12]:





$$2 \sqrt{\frac{1}{\Omega} \left(a + 2 b c^{2} + 3 c^{4} q\right)} \left(- \Omega c^{2} + \alpha + \frac{2 \beta}{\Omega} \left(a + 2 b c^{2} + 3 c^{4} q\right)\right)$$



In [13]:

    
simplify(series((solve((fc+fp+fcp).diff(p),E)[0]).subs(p,pc),c,n=7))









    Out[13]:





$$\frac{1}{\Omega a^{3}} \left(\Omega a^{3} \left(4 \beta \left(\frac{a}{\Omega}\right)^{\frac{3}{2}} + 2 \alpha \sqrt{\frac{a}{\Omega}} + \mathcal{O}\left(c^{7}\right)\right) + 2 a^{2} c^{2} \sqrt{\frac{a}{\Omega}} \left(- \Omega^{2} a + \Omega \alpha b + 6 a b \beta\right) + a c^{4} \sqrt{\frac{a}{\Omega}} \left(\Omega a \left(- 2 \Omega b + 3 \alpha q\right) - \Omega \alpha b^{2} + 18 a^{2} \beta q + 6 a b^{2} \beta\right) + c^{6} \sqrt{\frac{a}{\Omega}} \left(- 3 \Omega^{2} a^{2} q + \Omega a b \left(\Omega b - 3 \alpha q\right) + \Omega \alpha b^{3} + 18 a^{2} b \beta q - 2 a b^{3} \beta\right)\right)$$

Define above equation for plotting



In [14]:

    
def Ecc(a,b,q,c,Omega,alpha,beta):
    return [2*np.sqrt((a+2*b*c**2+3*q*c**4)/Omega)*(alpha-Omega*c**2+(2*beta/Omega)*(a+2*b*c**2+3*q*c**4)),
            (Omega*a**3*(4*beta*(a/Omega)**(3/2)+2*alpha*np.sqrt(a/Omega))
            +2*(a*c)**2*np.sqrt(a/Omega)*(-Omega**2*a+Omega*alpha*b+6*a*b*beta)
            +a*c**4*np.sqrt(a/Omega)*(Omega*a*(-2*Omega*b+3*alpha*q)-Omega*alpha*b**2+18*a**2*beta*q+6*a*b**2*beta)
            +c**60*np.sqrt(a/Omega)*(-3*q*(Omega*a)**2+Omega*a*b*(Omega*b-3*alpha*q)+Omega*alpha*b**3+18*a**2*b*beta*q-2*a*b**3*beta))/(Omega*a**3)]



In [15]:

    
def coeff0(a,alpha,beta,Omega):
    return 4*beta*(a/Omega)**(3/2)+2*alpha*np.sqrt(a/Omega)

def coeff2(a,Omega,alpha,b,beta):
    return 2*a**2*np.sqrt(a/Omega)*(-Omega**2*a+Omega*alpha*b+6*a*b*beta)/(Omega*a**3)

def coeff4(a, Omega, b, alpha, q, beta):
    return a*np.sqrt(a/Omega)*(Omega*a*(-2*Omega*b+3*alpha*q)-Omega*alpha*b**2+18*a**2*beta*q+6*a*b**2*beta)/(Omega*a**3)

def coeff6(a, Omega, b, alpha, q, beta):
    return np.sqrt(a/Omega)*(-3*q*(Omega*a)**2+Omega*a*b*(Omega*b-3*alpha*q)+Omega*alpha*b**3+18*a**2*b*beta*q-2*a*b**3*beta)/(Omega*a**3)



In [16]:

    
a = 6e5*(116-114.1)
Omega = 8e4
q = 1
b = 1.01e3
alpha = 0.5e3
beta = 1.05e6



In [18]:

    
coeff0(a,alpha,beta,Omega)









    Out[18]:





$$225932570.393$$



In [17]:

    
coeff2(a,Omega,alpha,b,beta)









    Out[17]:





$$-3488.45398211$$



In [488]:

    
coeff4(a, Omega, b, alpha, q, beta)









    Out[488]:





$$622.725336562$$



In [489]:

    
coeff6(a, Omega, b, alpha, q, beta)









    Out[489]:





$$0.153889123565$$



In [29]:

    
plt.figure(figsize=(11,8))
plt.plot(np.linspace(-17.5,17.5,201),Ecc(a,b,q,np.linspace(-17.5,17.5,201),Omega,alpha,beta)[0],label='$\mathregular{E}$')
# plt.plot(np.linspace(-17.5,17.5,201),Ecc(a,b,q,np.linspace(-17.5,17.5,201),Omega,alpha,beta)[1])
plt.scatter(0,Ecc(a,b,q,0,Omega,alpha,beta)[0],label='$\mathregular{E_{th}}$')
plt.ylim(2.259e8,2.2595e8),plt.xlim(-4,4)
plt.xlabel(c,fontsize=18)
plt.ylabel(E,fontsize=18,rotation='horizontal',labelpad=25)
plt.legend(loc='lower right',fontsize=18);



In [497]:

    
150000+2.258e8









    Out[497]:





$$225950000.0$$



In [441]:

    
plt.figure(figsize=(11,8))
plt.scatter(np.linspace(-17.5,17.5,201),Ecc(6e4*(116-114.1),10**4,9e7,np.linspace(-17.5,17.5,201),10**10,0.6e9,10**4)[1],color='r',label='expand')
plt.ylim(-0.5e14,0.25e14),plt.xlim(-5,5)
plt.xlabel('c',fontsize=18)
plt.ylabel('E',rotation='horizontal',labelpad=25,fontsize=18);



In [20]:

    
def quadterm(alpha,Omega,Eth,E,a):
    return [np.sqrt((Eth-E)/np.sqrt(4*Omega*a)),-np.sqrt((Eth-E)/np.sqrt(4*Omega*a))]



In [21]:

    
# np.linspace(0,60e6,1000)



In [23]:

    
plt.plot(np.linspace(0,60e6,1000),quadterm(0.6e9,10**10,23663663.66366366,np.linspace(0,60e6,1000),6e4*(116-114.1))[0],'b')
plt.plot(np.linspace(0,60e6,1000),quadterm(0.6e9,10**10,23663663.66366366,np.linspace(0,60e6,1000),6e4*(116-114.1))[1],'r');









    



/usr/local/lib/python3.4/dist-packages/IPython/kernel/__main__.py:2: RuntimeWarning: invalid value encountered in sqrt
  from IPython.kernel.zmq import kernelapp as app



In [50]:

    
plt.figure(figsize=(11,8))
plt.scatter(Ecc(6e4*(116-114.1),10**4,9e7,np.linspace(7.5,17.5),10**10,0.6e9,10**4)[0],np.linspace(7.5,17.5),color='r',label='no expand')
plt.ylim(0,35)
plt.xlabel('E',fontsize=18)
plt.ylabel('c',rotation='horizontal',labelpad=25,fontsize=18);

$\dfrac{\partial f}{\partial p} = 0$



In [17]:

    
(fc+fp+fcp).diff(p)









    Out[17]:





$$- E - 2 \Omega c^{2} p + 2 \alpha p + 4 \beta p^{3}$$

Solve for $c$



In [40]:

    
solve((fc+fp+fcp).diff(p),c)[1]









    Out[40]:





$$\frac{\sqrt{2}}{2} \sqrt{\frac{1}{\Omega p} \left(- E + 2 \alpha p + 4 \beta p^{3}\right)}$$

Sub $c(p)$ into $p(c)$



In [50]:

    
expand(simplify(pc.subs(c,solve((fc+fp+fcp).diff(p),c)[1])),p)









    Out[50]:





$$\frac{1}{2} \sqrt{\frac{3 E^{2} q}{\Omega^{3} p^{2}} - \frac{4 E b}{\Omega^{2} p} - \frac{12 E \alpha q}{\Omega^{3} p} - \frac{24 E}{\Omega^{3}} \beta p q + \frac{4 a}{\Omega} + \frac{8 \alpha}{\Omega^{2}} b + \frac{16 b}{\Omega^{2}} \beta p^{2} + \frac{12 q}{\Omega^{3}} \alpha^{2} + \frac{48 \alpha}{\Omega^{3}} \beta p^{2} q + \frac{48 q}{\Omega^{3}} \beta^{2} p^{4}}$$

Solve for $E(p)$



In [52]:

    
simplify(solve(expand(simplify(pc.subs(c,solve((fc+fp+fcp).diff(p),c)[1])),p)-p,E)[1])









    Out[52]:





$$\frac{1}{3 q} \left(2 p \left(\Omega b + 3 \alpha q + 6 \beta p^{2} q\right) + 2 \sqrt{\Omega^{2} p^{2} \left(3 \Omega p^{2} q - 3 a q + b^{2}\right)}\right)$$

Define above equation for plotting



In [46]:

    
def Epp(a,b,q,alpha,beta,Omega,p):
    return (2*p*(Omega*b+3*alpha*q+6*beta*q*p**2)+2*Omega*p*np.sqrt(3*Omega*q*p**2-3*a*q+b**2))/(3*q)



In [55]:

    
plt.figure(figsize=(11,8))
plt.scatter(Epp(6e4*(116-114.1),10**4,9e7,0.6e9,10**4,10**11,np.linspace(0,0.005))/10**7,np.linspace(0,0.005),color='r',label='T = 116')
plt.scatter(Epp(6e4*(110-114.1),10**4,9e7,0.6e9,10**4,10**11,np.linspace(0,0.005))/10**7,np.linspace(0,0.005),color='g',label='T = 110')
plt.scatter(Epp(6e4*(105-114.1),10**4,9e7,0.6e9,10**4,10**11,np.linspace(0,0.005))/10**7,np.linspace(0,0.005),color='b',label='T = 105')
plt.ylim(0,0.006),plt.xlim(0,20)
plt.xlabel('E',fontsize=18)
plt.ylabel('p',rotation='horizontal',labelpad=25,fontsize=18)
plt.legend(loc='upper right',fontsize=16);









    



/usr/local/lib/python3.4/dist-packages/IPython/kernel/__main__.py:2: RuntimeWarning: invalid value encountered in sqrt
  from IPython.kernel.zmq import kernelapp as app



In [19]:

    
def ppp(a,b,q,Omega,c):
    return np.sqrt((a+2*b*c**2+3*q*c**4)/Omega)



In [20]:

    
plt.figure(figsize=(11,8))
plt.scatter(np.linspace(0,35),ppp(6e4*(116-114.1),10**4,9e7,10**10,np.linspace(0,35)))
plt.ylim(-2,220),plt.xlim(-1,36)
plt.xlabel('c',fontsize=18)
plt.ylabel('p',rotation='horizontal',labelpad=25,fontsize=18);



In [9]:

    
simplify((solve((fp+fcp).diff(p),E)[0]).subs(p,solve((fc+fp+fcp).diff(c),p)[1]))









    Out[9]:





$$\frac{2}{\Omega} \sqrt{\frac{1}{\Omega} \left(a + 2 b c^{2} + 3 c^{4} q\right)} \left(\Omega \left(- \Omega c^{2} + \alpha\right) + 2 \beta \left(a + 2 b c^{2} + 3 c^{4} q\right)\right)$$



In [64]:

    
def Ec(om, a, b, c, q, alp, beta):
    return (2/om)*np.sqrt((a+2*b*c**2+3*q*c**4)/om)*(alp*om-(c*om)**2+2*beta*a+4*b*beta*c**2+6*beta*q*c**4)



In [86]:

    
expand((alpha-E/(2*p)+2*beta*p**2)**2)









    Out[86]:





$$\frac{E^{2}}{4 p^{2}} - \frac{E \alpha}{p} - 2 E \beta p + \alpha^{2} + 4 \alpha \beta p^{2} + 4 \beta^{2} p^{4}$$



In [65]:

    
plt.scatter(Ec(10**11,6e4,10**4,np.linspace(0,15,200),9e7,0.6e9,10**4),np.linspace(0,15,200),marker='s');



In [48]:

    
def Eth(T):
    E = []
    for i in T:
        if i > 114.1:
            E.append(0.64e6*(i-100)*np.sqrt(i-114.1))
        else:
            E.append(0)
    return E



In [58]:

    
plt.plot(np.linspace(113,119,10),Eth(np.linspace(113,119,10)))
plt.xlim(113,122),plt.ylim(0,3e7);



In [ ]: