In [1]:

%matplotlib notebook
from sympy import *
from sympy import besselk, besselj
from sympy.utilities.lambdify import lambdify
from sympy.utilities.autowrap import ufuncify
from sympy.utilities.lambdify import implemented_function
from sympy.physics.hydrogen import R_nl
import numpy as np
from scipy.special import kn, jv, kv, iv, yv, k0, k1
from scipy.integrate import quad
import matplotlib.pyplot as plt
from scipy import misc
init_printing(use_unicode=True)




In [2]:

x, y, r, R, E, K, phi, theta_s, a, R_L, R_R, d, p, v, theta_1, theta_2 = symbols('x y r R E K phi theta_s a R_L R_R d p v theta_1 theta_2')



## Single Cylinder



In [22]:

def Nx(x,y,theta_s):
return -theta_s*(y/(x**2+y**2))

def Ny(x,y,theta_s):
return theta_s*(x/(x**2+y**2))




In [23]:

Nx(x,y,theta_s)




Out[23]:

$$- \frac{\theta_{s} y}{x^{2} + y^{2}}$$




In [24]:

Ny(x,y,theta_s)




Out[24]:

$$\frac{\theta_{s} x}{x^{2} + y^{2}}$$




In [25]:

DXNX = Nx(x,y,theta_s).diff(x)
DYNX = Nx(x,y,theta_s).diff(y)
DXNY = Ny(x,y,theta_s).diff(x)
DYNY = Ny(x,y,theta_s).diff(y)




In [26]:

def single_f(x, y, theta_s, dxnx, dynx, dxny, dyny):
return (dxnx**2 + dynx**2 + dxny**2 + dyny**2 + 2*dxnx*dyny - 2*dxny*dynx)/2




In [27]:

simplify(single_f(x, y, theta_s, DXNX, DYNX, DXNY, DYNY))




Out[27]:

$$0$$



## Single Cylinder E-Field



In [3]:

def nx(r, phi, theta_s, R, d):
return (-theta_s*sin(phi)*besselk(1, r/d)/besselk(1, R/d))

def ny(r, phi, theta_s, R, d):
return (theta_s*cos(phi)*besselk(1, r/d)/besselk(1, R/d))




In [4]:

nxx = nx(r, phi, theta_s, 1, d)
nxx




Out[4]:

$$- \frac{\theta_{s} K_{1}\left(\frac{r}{d}\right)}{K_{1}\left(\frac{1}{d}\right)} \sin{\left (\phi \right )}$$




In [5]:

nyy = ny(r, phi, theta_s, 1, d)
nyy




Out[5]:

$$\frac{\theta_{s} K_{1}\left(\frac{r}{d}\right)}{K_{1}\left(\frac{1}{d}\right)} \cos{\left (\phi \right )}$$




In [6]:

dxnx = cos(phi)*nx(r, phi, theta_s, 1, d).diff(r) - (sin(phi)*nx(r, phi, theta_s, 1, d).diff(phi))/r
dynx = sin(phi)*nx(r, phi, theta_s, 1, d).diff(r) + (cos(phi)*nx(r, phi, theta_s, 1, d).diff(phi))/r
dxny = cos(phi)*ny(r, phi, theta_s, 1, d).diff(r) - (sin(phi)*ny(r, phi, theta_s, 1, d).diff(phi))/r
dyny = sin(phi)*ny(r, phi, theta_s, 1, d).diff(r) + (cos(phi)*ny(r, phi, theta_s, 1, d).diff(phi))/r



### $(\partial_{x}n_{x}+\partial_{y}n_{y})^2$



In [23]:

(d*E)**2*(dxnx+dyny)**2




Out[23]:

$$0$$



### $\dfrac{E^2\theta_s^2}{2}(\partial_{x}n_{y}-\partial_{y}n_{x})^2$



In [22]:

simplify((d*E)**2*(dxny-dynx)**2/2)




Out[22]:

$$\frac{E^{2} \theta_{s}^{2} K^{2}_{0}\left(\frac{r}{d}\right)}{2 K^{2}_{1}\left(\frac{1}{d}\right)}$$



### $\dfrac{1}{2}(n_{x}^2 + n_{y}^2)$



In [20]:

simplify((nxx**2+nyy**2)/(2))




Out[20]:

$$\frac{\theta_{s}^{2} K^{2}_{1}\left(\frac{r}{d}\right)}{2 K^{2}_{1}\left(\frac{1}{d}\right)}$$




In [18]:

def f(r, phi, theta_s, R, d, E, dxnx, dynx, dxny, dyny):
return (E**2/2)*(d**2*((dxnx+dyny)**2 + (dxny-dynx)**2) + (nx(r, phi, theta_s, 1, d)**2 + ny(r, phi, theta_s, 1, d)**2))




In [19]:

simplify(f(r, phi, theta_s, 1, d, E, dxnx, dynx, dxny, dyny)).




Out[19]:

$$\frac{E^{2} \theta_{s}^{2}}{2 K^{2}_{1}\left(\frac{1}{d}\right)} \left(K^{2}_{0}\left(\frac{r}{d}\right) + K^{2}_{1}\left(\frac{r}{d}\right)\right)$$




In [17]:

Ef = lambda r,theta_s,d: (theta_s**2*(k0(r/d)**2 + k1(r/d)**2) - k1(1/d)**2)/(2*(d*k1(1/d))**2)




In [48]:

splay = lambda r: 0




In [34]:

twist = lambda r,theta_s,d: (theta_s*k0(r/d)/(d*k1(1/d)))**2/2




In [43]:

field = lambda r,theta_s,d: ((theta_s*k1(r/d))**2-k1(1/d)**2)/(2*(d*k1(1/d))**2)



### $d = 2$, $\theta_s = 0.2$



In [51]:

Is, es = quad(splay,1,np.inf)
Is




Out[51]:

$$0.0$$




In [52]:

It, et = quad(twist,1,np.inf,args=(0.2,2))
It




Out[52]:

$$0.0010043517423015916$$




In [53]:

If, ef = quad(field,1,np.inf,args=(0.2,2))
If




Out[53]:

$$0.12743840739612153$$




In [54]:

IE, e = quad(Ef, 1, np.inf, args=(0.2, 2))
IE




Out[54]:

$$0.128442759133668$$




In [34]:

integrate((0.2/(2*besselk(1,1/2)))**2*(besselk(0,r/2)**2+besselk(1,r/2)**2), (r,1,oo))




Out[34]:

$$\frac{1}{K^{2}_{1}\left(0.5\right)} \left(0.01 \int_{1}^{\infty} K^{2}_{0}\left(\frac{r}{2}\right)\, dr + 0.01 \int_{1}^{\infty} K^{2}_{1}\left(\frac{r}{2}\right)\, dr\right)$$




In [35]:

integrate((theta_s/(d*besselk(1,1/d)))**2*(besselk(0,r/d)**2+besselk(1,r/d)**2), (r,1,oo))




Out[35]:

$$\frac{\theta_{s}^{2}}{d^{2} K^{2}_{1}\left(\frac{1}{d}\right)} \left(\int_{1}^{\infty} K^{2}_{0}\left(\frac{r}{d}\right)\, dr + \int_{1}^{\infty} K^{2}_{1}\left(\frac{r}{d}\right)\, dr\right)$$



## Two Cylinders E-Field



In [3]:

def nx_cart(r, x, y, theta_s, R, d):
return (-theta_s*y*besselk(1, r/d))/(r*besselk(1, R/d))

def ny_cart(r, x, y, theta_s, R, d):
return (theta_s*x*besselk(1, r/d))/(r*besselk(1, R/d))




In [ ]:

a = 2
d = 1
R_L = 1
R_R = 1
theta_s = 0.2




In [4]:

nx_L = nx_cart(sqrt((x+a)**2 + y**2), x, y, theta_2, 1, d)
ny_L = ny_cart(sqrt((x+a)**2 + y**2), x+a, y, theta_2, 1, d)
nx_R = nx_cart(sqrt((x-a)**2 + y**2), x, y, theta_1, 1, d)
ny_R = ny_cart(sqrt((x-a)**2 + y**2), x-a, y, theta_1, 1, d)




In [5]:

nx_L, nx_R




Out[5]:

$$\left ( - \frac{\theta_{2} y K_{1}\left(\frac{1}{d} \sqrt{y^{2} + \left(a + x\right)^{2}}\right)}{\sqrt{y^{2} + \left(a + x\right)^{2}} K_{1}\left(\frac{1}{d}\right)}, \quad - \frac{\theta_{1} y K_{1}\left(\frac{1}{d} \sqrt{y^{2} + \left(- a + x\right)^{2}}\right)}{\sqrt{y^{2} + \left(- a + x\right)^{2}} K_{1}\left(\frac{1}{d}\right)}\right )$$




In [ ]:

#ny_L, ny_R




In [20]:

nx_tot = nx_L + nx_R
ny_tot = ny_L + ny_R




In [31]:

(nx_tot**2 + ny_tot**2)/d**2




Out[31]:

$$\frac{1}{d^{2}} \left(\left(- \frac{\theta_{1} y K_{1}\left(\frac{1}{d} \sqrt{y^{2} + \left(- a + x\right)^{2}}\right)}{\sqrt{y^{2} + \left(- a + x\right)^{2}} K_{1}\left(\frac{1}{d}\right)} - \frac{\theta_{2} y K_{1}\left(\frac{1}{d} \sqrt{y^{2} + \left(a + x\right)^{2}}\right)}{\sqrt{y^{2} + \left(a + x\right)^{2}} K_{1}\left(\frac{1}{d}\right)}\right)^{2} + \left(\frac{\theta_{1} \left(- a + x\right) K_{1}\left(\frac{1}{d} \sqrt{y^{2} + \left(- a + x\right)^{2}}\right)}{\sqrt{y^{2} + \left(- a + x\right)^{2}} K_{1}\left(\frac{1}{d}\right)} + \frac{\theta_{2} \left(a + x\right) K_{1}\left(\frac{1}{d} \sqrt{y^{2} + \left(a + x\right)^{2}}\right)}{\sqrt{y^{2} + \left(a + x\right)^{2}} K_{1}\left(\frac{1}{d}\right)}\right)^{2}\right)$$




In [21]:

dxnx_2cyl = nx_tot.diff(x)
dynx_2cyl = nx_tot.diff(y)
dxny_2cyl = ny_tot.diff(x)
dyny_2cyl = ny_tot.diff(y)



## First term



In [22]:

simplify((dxnx_2cyl + dyny_2cyl)**2)




Out[22]:

$$0$$



## Second term



In [36]:

simplify((dxny_2cyl - dynx_2cyl)**2/2)




Out[36]:

$$\frac{1}{2 d^{2} K^{2}_{1}\left(\frac{1}{d}\right)} \left(\theta_{1} K_{0}\left(\frac{1}{d} \sqrt{a^{2} - 2 a x + x^{2} + y^{2}}\right) + \theta_{2} K_{0}\left(\frac{1}{d} \sqrt{a^{2} + 2 a x + x^{2} + y^{2}}\right)\right)^{2}$$



## n term



In [39]:

simplify((nx_tot**2 + ny_tot**2 - 1)/(2*d**2))




Out[39]:

$$\frac{1}{2 d^{2} \left(y^{2} + \left(a - x\right)^{2}\right) \left(y^{2} + \left(a + x\right)^{2}\right) K^{2}_{1}\left(\frac{1}{d}\right)} \left(y^{2} \left(\theta_{1} \sqrt{y^{2} + \left(a + x\right)^{2}} K_{1}\left(\frac{1}{d} \sqrt{y^{2} + \left(a - x\right)^{2}}\right) + \theta_{2} \sqrt{y^{2} + \left(a - x\right)^{2}} K_{1}\left(\frac{1}{d} \sqrt{y^{2} + \left(a + x\right)^{2}}\right)\right)^{2} - \left(y^{2} + \left(a - x\right)^{2}\right) \left(y^{2} + \left(a + x\right)^{2}\right) K^{2}_{1}\left(\frac{1}{d}\right) + \left(\theta_{1} \left(a - x\right) \sqrt{y^{2} + \left(a + x\right)^{2}} K_{1}\left(\frac{1}{d} \sqrt{y^{2} + \left(a - x\right)^{2}}\right) - \theta_{2} \left(a + x\right) \sqrt{y^{2} + \left(a - x\right)^{2}} K_{1}\left(\frac{1}{d} \sqrt{y^{2} + \left(a + x\right)^{2}}\right)\right)^{2}\right)$$




In [ ]:

def f_2(n1, n2, theta_s, a, R_L, R_R, d, dxnx, dynx, dxny, dyny):
return (1/2)*(dxnx**2 + dynx**2 + dxny**2 + dyny**2 + 2*dxnx*dyny - 2*dxny*dynx) + (1/2)*(n1**2 + n2**2) - (1/2)




In [ ]:

f_2cyl = f_2(nx_tot, ny_tot, theta_s, a, R_L, R_R, d, dxnx_2cyl, dynx_2cyl, dxny_2cyl, dyny_2cyl)




In [ ]:

I, e = quad(simple_f, 1, 200, args=(0.2, 10, 1, 500))




In [ ]:

I




In [ ]: