``````

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 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.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.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.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.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_angles = 50

angles_single = np.linspace(0,2*np.pi,n_angles,endpoint=False)

``````
``````

In [5]:

``````
``````

In [6]:

``````
``````

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)

``````
``````

In [8]:

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

``````
``````

In [9]:

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();

``````
``````

``````