notebook.community

Edit and run 





In [1]:



    


# %matplotlib notebook
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.axis as ax
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import LightSource
from matplotlib.tri import Triangulation, UniformTriRefiner, CubicTriInterpolator
import numpy as np
from scipy.special import kn, k1, k0, jn
from scipy.integrate import quad
from matplotlib import cm



    











In [2]:



    


def r(x,y):
    return np.sqrt(x**2 + y**2)

def Nx(x,y,thetaS):
    return -y*thetaS/r(x,y)**2

def Ny(x,y,thetaS):
    return x*thetaS/r(x,y)**2

def n_perp(x, y, thetaS):
    return np.sqrt(Nx(x, y, thetaS)**2 + Ny(x, y, thetaS)**2)

def Nx_E(x, y, d, r1, thetaS):
    return -(thetaS*y*k1(r(x,y)*d**-1)/(r(x,y)*k1(r1*d**-1)))

def Ny_E(x, y, d, r1, thetaS):
    return (thetaS*x*k1(r(x,y)*d**-1)/(r(x,y)*k1(r1*d**-1)))

def nE_perp(x, y, d, r1, thetaS):
    return np.sqrt(Nx_E(x, y, d, r1, thetaS)**2 + Ny_E(x, y, d, r1, thetaS)**2)

def simple_f(x, y, d, r1, thetaS, K):
    return K*(thetaS/(d*k1(r1/d)))**2*(k0(r(x,y)/d)**2 + k1(r(x,y)/d)**2)/2

def simple_f2(r, d, r1, thetaS, K):
    return K*(thetaS/(d*k1(r1/d)))**2*(k0(r/d)**2 + k1(r/d)**2)/2



    











In [3]:



    


xaxis = np.zeros(100)
yaxis = np.zeros(100)
x = np.linspace(1,10, 100)
y = np.linspace(1,10, 100)

r1 = 1
dlist = [1,2,3,4]

colors = cm.gist_rainbow(np.linspace(0,1,len(dlist)))

dlabels = ['d = 1','d = 2','d = 3','d = 4']

No field





In [56]:



    


# This is the single cylinder case with no field
plt.figure(figsize=(10,8))
plt.plot(r(x,yaxis), n_perp(x,yaxis,0.2)) # plots the nx solution for each value of thetaS
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
plt.title('Zero Field $\mathregular{|n_{\perp}|}$ vs r', fontsize=18,y=1.01)
plt.xlabel('r', fontsize=20)
plt.ylabel('$\mathregular{|n_{\perp}|}$', fontsize=20, rotation='horizontal', labelpad=30)
plt.ylim(0,0.2),plt.xlim(0,7)
plt.tick_params(axis='both',which='major',labelsize=14)
plt.tick_params(axis='x', top='off')
plt.tick_params(axis='y', right='off');

E-Field





In [161]:



    


#This is the E-field case with varying d along y-axis
plt.figure(figsize=(10,8))
ax = plt.subplot(111)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.plot(r(xaxis,y), n_perp(xaxis,y,0.2), '--k', label='E = 0')
for color, dlabel, dval in zip(colors,dlabels,dlist):
    plt.plot(r(xaxis,y), nE_perp(xaxis,y,dval,r1,0.2), color=color)
    plt.plot(0,0, color=color, label=dlabel)
plt.title('E-field $\mathregular{|n_{\perp}|}$ vs r', fontsize=20,y=1.01)
plt.xlabel('r', fontsize=20)
plt.ylabel('$\mathregular{|n_{\perp}|}$', fontsize=20, rotation='horizontal', labelpad=30)
plt.legend(loc='best',fontsize=16)
plt.ylim(0,0.2),plt.xlim(0,7)
plt.tick_params(axis='both',which='major',labelsize=14)
plt.tick_params(axis='x', top='off')
plt.tick_params(axis='y', right='off');



    


















    
























In [64]:



    


# Free energy density of the E-field case w/ varying d-values, constant thetaS
plt.figure(figsize=(10,8))
ax = plt.subplot(111)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
for color, dlabel, dval in zip(colors,dlabels,dlist):
    plt.plot(np.linspace(1,10,200), simple_f2(np.linspace(1,10,200),dval,r1,0.2,1), color=color,label=dlabel)
plt.title('E-field f vs r', fontsize=20,y=1.01)
plt.xlabel('r', fontsize=20)
plt.ylabel('f', fontsize=20, rotation='horizontal', labelpad=25)
plt.xlim(0.5,2.5),plt.ylim(0,0.03)
# xticks = plt.gca().xaxis.get_major_ticks()
# xticks[0].label1.set_visible(False)
plt.tick_params(axis='both',which='major',labelsize=14)
plt.tick_params(axis='x', top='off')
plt.tick_params(axis='y', right='off')
plt.legend(loc='best',fontsize=16);



    


















    
























In [4]:



    


def int_f(d,theta_s,K,R):
    return K*R*np.pi*theta_s**2*k0(R/d)/(d*k1(R/d))



    











In [47]:



    


dvalues = np.linspace(0.1,10,100)

plt.figure(figsize=(11,8))
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
plt.plot(dvalues,int_f(dvalues,0.2,1,1))
plt.ylim(0,0.6),plt.xlim(0,10.1)
plt.xlabel('d',fontsize=20)
plt.ylabel('F/L',fontsize=20,labelpad=25,rotation='horizontal')
plt.title('Total Free Energy vs d',fontsize=20,y=1.01)
plt.tick_params(axis='both',which='major',labelsize=14)
plt.tick_params(axis='x', top='off')
plt.tick_params(axis='y', right='off');

2 Cyl E-Field

Line plot (maybe not)





In [79]:



    


ds = [1,2,3,10,50]
dslabels = ['d = 1','d = 2','d = 3','d = 10','d = 50']
colors2 = cm.gnuplot(np.linspace(0,1,5))



    











In [27]:



    


def twod(x,y,a,d,theta_1,theta_2):
    return surf_eqn(x,y,a,d,theta_1,theta_2) + n_eqn(x,y,a,d,theta_1,theta_2)



    











In [183]:



    


xxx = np.arange(1,10.02,0.02)
yyy = xxx.copy()



    











In [135]:



    


plt.plot(xxx,twod(xxx,yyy,a,1,0.2,0.2));

Single Cyl 3D





In [4]:



    


n_radii = 30
n_angles = 50

radii_single = np.linspace(1,10,n_radii)
angles_single = np.linspace(0,2*np.pi,n_angles,endpoint=False)



    











In [5]:



    


angles_single = np.repeat(angles_single[...,np.newaxis],n_radii,axis=1)



    











In [6]:



    


x_single = (radii_single*np.cos(angles_single)).flatten()
y_single = (radii_single*np.sin(angles_single)).flatten()



    











In [7]:



    


tri = Triangulation(x_single,y_single)
xmid = x_single[tri.triangles].mean(axis=1)
ymid = y_single[tri.triangles].mean(axis=1)
mask = np.where(xmid*xmid + ymid*ymid < 1, 1, 0)
tri.set_mask(mask)



    











In [8]:



    


work = n_perp(tri.x,tri.y,0.2)



    











In [9]:



    


xx = np.append(-radii_single[::-1],radii_single)
yy = xx.copy()



    











In [10]:



    


pos = np.where(abs(xx)==1)[0]



    











In [11]:



    


xx[pos] = np.nan
yy[pos] = np.nan



    











In [12]:



    


n_perp(xx,yy,0.2)



    


















    
Out[12]:








array([ 0.01414214,  0.01459509,  0.01507801,  0.01559399,  0.01614653,
        0.01673967,  0.01737805,  0.01806705,  0.01881293,  0.01962306,
        0.0205061 ,  0.02147235,  0.02253417,  0.02370647,  0.02500743,
        0.02645948,  0.02809054,  0.02993591,  0.03204078,  0.03446403,
        0.03728381,  0.04060613,  0.04457847,  0.04941228,  0.05542188,
        0.06309568,  0.07323606,  0.08725999,  0.10792682,         nan,
               nan,  0.10792682,  0.08725999,  0.07323606,  0.06309568,
        0.05542188,  0.04941228,  0.04457847,  0.04060613,  0.03728381,
        0.03446403,  0.03204078,  0.02993591,  0.02809054,  0.02645948,
        0.02500743,  0.02370647,  0.02253417,  0.02147235,  0.0205061 ,
        0.01962306,  0.01881293,  0.01806705,  0.01737805,  0.01673967,
        0.01614653,  0.01559399,  0.01507801,  0.01459509,  0.01414214])

















In [19]:



    


fignx = plt.figure(figsize=(11,8))
axnx = fignx.gca(projection='3d')
surfnx = axnx.plot_trisurf(tri,work, cmap=cm.jet,linewidth=0)
axnx.set_aspect('equal')
axnx.set_xlabel('x', fontsize=14)
axnx.set_ylabel('y', fontsize=14)
axnx.set_zlabel('nx', fontsize=14)
fignx.colorbar(surfnx,shrink=0.7);

Single Cyl E-field





In [14]:



    


def simple_ff(r, theta_s, d, R, K):
    return K*(theta_s/(d*k1(R/d)))**2*(k0(r/d)**2 + k1(r/d)**2)/2



    











In [15]:



    


Z = simple_ff(np.sqrt(tri.x**2+tri.y**2), 0.2, 1, 1, 1)



    











In [16]:



    


fig2 = plt.figure(figsize=(11,8))
ax2 = fig2.gca(projection='3d')
surf2 = ax2.plot_trisurf(tri,Z, cmap=cm.jet,linewidth=0)
ax2.set_aspect('equal')
ax2.set_xlabel('x', fontsize=14)
ax2.set_ylabel('y', fontsize=14)
ax2.set_zlabel('f', fontsize=14)
ax2.set_title('Electric Field Free Energy Density', fontsize=18)
fig2.colorbar(surf2,shrink=0.7)
fig2.tight_layout();

Content source: brettavedisian/Liquid-Crystals-Summer-2015

Similar notebooks:

notebook.community | gallery | about