<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Work on getting $E$ vs $\theta$ plot parameters w/ 3 extrema
Step1:  Generalized Landau Model of Ferroelectric Liquid Crystals
Step2:  $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
Step3:  Electrically Induced Tilt in Achiral Bent-Core Liquid Crystals
Step4:  Mapping
Step5:  $\dfrac{\partial f}{\partial c} = 0$
Step6:  Solve for $p$
Step7:  Solve $\dfrac{\partial f}{\partial p} = 0$ for $E$
Step8:  Sub $p(c)$ into $E(c,p)$
Step9:  Define above equation for plotting
Step10:  $\dfrac{\partial f}{\partial p} = 0$
Step11:  Solve for $c$
Step12:  Sub $c(p)$ into $p(c)$
Step13:  Solve for $E(p)$
Step14:  Define above equation for plotting
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from sympy import *
import matplotlib.pyplot as plt
import numpy as np
init_printing(use_unicode=True)
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')
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')
Explanation: Work on getting $E$ vs $\theta$ plot parameters w/ 3 extrema
End of explanation
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
g1(a,b,q,xi,P_x,varep,eta,C,Omega)
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
g
Explanation: Generalized Landau Model of Ferroelectric Liquid Crystals
End of explanation
g.subs([(epsilon,1/r_p),(eta,u_p),(C,-gamma),(Omega,2*e),(xi,c),(P_x,p)])
Explanation: $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
End of explanation
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
fc
fp
fcp
Explanation: Electrically Induced Tilt in Achiral Bent-Core Liquid Crystals
End of explanation
# fc = fc.subs([(a_0*(T-T_c),a)])#,6e4),(b,10**4),(q,9e7)])
# fc
fp = fp.subs([(gamma,0)])#,(alpha_0*(T-T_p),alpha),(beta,10**4)])
fp
Explanation: 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}$
End of explanation
collect((fc+fp+fcp).diff(c),c)
Explanation: $\dfrac{\partial f}{\partial c} = 0$
End of explanation
pc = solve((fc+fp+fcp).diff(c),p)[1]
pc
Explanation: Solve for $p$
End of explanation
solve((fc+fp+fcp).diff(p),E)[0]
Explanation: Solve $\dfrac{\partial f}{\partial p} = 0$ for $E$
End of explanation
(solve((fc+fp+fcp).diff(p),E)[0]).subs(p,pc)
simplify(series((solve((fc+fp+fcp).diff(p),E)[0]).subs(p,pc),c,n=7))
Explanation: Sub $p(c)$ into $E(c,p)$
End of explanation
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)]
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)
a = 6e5*(116-114.1)
Omega = 8e4
q = 1
b = 1.01e3
alpha = 0.5e3
beta = 1.05e6
coeff0(a,alpha,beta,Omega)
coeff2(a,Omega,alpha,b,beta)
coeff4(a, Omega, b, alpha, q, beta)
coeff6(a, Omega, b, alpha, q, beta)
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);
150000+2.258e8
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);
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))]
# np.linspace(0,60e6,1000)
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');
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);
Explanation: Define above equation for plotting
End of explanation
(fc+fp+fcp).diff(p)
Explanation: $\dfrac{\partial f}{\partial p} = 0$
End of explanation
solve((fc+fp+fcp).diff(p),c)[1]
Explanation: Solve for $c$
End of explanation
expand(simplify(pc.subs(c,solve((fc+fp+fcp).diff(p),c)[1])),p)
Explanation: Sub $c(p)$ into $p(c)$
End of explanation
simplify(solve(expand(simplify(pc.subs(c,solve((fc+fp+fcp).diff(p),c)[1])),p)-p,E)[1])
Explanation: Solve for $E(p)$
End of explanation
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)
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);
def ppp(a,b,q,Omega,c):
    return np.sqrt((a+2*b*c**2+3*q*c**4)/Omega)
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);
simplify((solve((fp+fcp).diff(p),E)[0]).subs(p,solve((fc+fp+fcp).diff(c),p)[1]))
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)
expand((alpha-E/(2*p)+2*beta*p**2)**2)
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');
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
plt.plot(np.linspace(113,119,10),Eth(np.linspace(113,119,10)))
plt.xlim(113,122),plt.ylim(0,3e7);
Explanation: Define above equation for plotting
End of explanation
<END_TASK>