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')


    

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$$

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


    

In [23]:

    

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


    

    
Out[23]:





$$0$$

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)}$$

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)


    

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)$$

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)


    

In [22]:

    

simplify((dxnx_2cyl + dyny_2cyl)**2)


    

    
Out[22]:





$$0$$

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}$$

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 [ ]: