\n", "The magnetic field of two parallel coils carrying currents in the same direction is calculated by evaluating the complete integrals of first and second kinds calculated numerically.\n", "\n", "$$ \n", "B_x = B_0 \\frac{1}{\\pi \\sqrt{Q}}\\bigl[E(k) \\frac{1-\\alpha^2-\\beta^2}{Q-4\\alpha} +K(k)\\bigr ], \\ \\ \\ \\ \n", "B_r = B_0 \\frac{\\gamma}{\\pi \\sqrt{Q}}\\bigl[E(k) \\frac{1+\\alpha^2+\\beta^2}{Q-4\\alpha} -K(k)\\bigr ]\n", " \\\\\n", "\\alpha = \\frac{r}{a}, \\ \\ \\ \\beta = \\frac{x}{a}, \\ \\ \\ \\gamma = \\frac{x}{r}, \\ \\ \\ Q = (1+\\alpha)^2 +\\beta^2, \\ \\ \\ k = \\sqrt{\\frac{4\\alpha}{Q}}, \\\\\n", " K(k) = \\int_0^{\\pi/2} \\, \\frac{d\\theta}{\\sqrt{1-k^2\\sin^2 \\theta}}, \\ \\ \\ \n", " E(k) = \\int_0^{\\pi/2} \\, \\sqrt{1-k^2\\sin^2 t}\\,dt.\n", " $$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "% matplotlib notebook\n", "\n", "from IPython.display import Image\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "ra = 0.1 # coil radius\n", "ra2 = ra**2\n", "nmax = 100\n", "dx = 0.01 # take dx = dy\n", "n = 300\n", "nm = int(nmax/2)\n", "nplot = 200\n", "Bx = np.zeros((nmax,nmax),float)\n", "By = np.zeros((nmax,nmax),float)\n", "xx = np.zeros(nplot,float)\n", "yy = np.zeros(nplot,float)\n", "tt = np.zeros(nplot,float)\n", "\n", "def funct(k,t):\n", " f = np.sqrt(1-k*np.sin(t)**2)\n", " return f\n", "\n", "def funth(k,th): # integrand elliptic integral\n", " y = np.sqrt(1-k*np.sin(th)**2)\n", " return 1./y \n", " \n", "def trap(a,b,n,func,x): # trapezoidal rule limits: a, b. n = points\n", " h = (b-a)/(n-1) # step, N-1 intervals\n", " sum = (func(x,a)+func(x,b))/2 # first, last mulp. by h/2\n", " for i in range(1, n -1): \n", " sum += func(x,a+i*h)\n", " return h*sum\n", "\n", "def Ellipints(k): # 1st and 2nd kind elliptic int\n", " a = 0\n", " b = np.pi/2.\n", " Ek = trap(a,b,n,funct,k)\n", " Kk = trap(a,b,n,funth,k)\n", " return Kk,Ek \n", "\n", "def Bfield(dx,x,y): \n", " alpha = y/ra\n", " beta = x/ra\n", " gamma = x/y\n", " xx = 1.-x\n", " betp = xx/ra\n", " gamp = xx/y\n", " term = (1+alpha)**2\n", " Q = term + beta**2\n", " Qp = term + betp**2\n", " Qsq = np.sqrt(Q)\n", " Qpsq = np.sqrt(Qp)\n", " k2 = 4*alpha/Q\n", " kp2 = 4*alpha/Qp\n", " Kk,Ek = Ellipints(k2)\n", " Kkp,Ekp = Ellipints(kp2)\n", " fac = Ek/(Q-4*alpha)\n", " albet = alpha*alpha+beta*beta\n", " facp = Ekp/(Qp-4*alpha)\n", " albetp = alpha*alpha+betp*betp\n", " Bxp = (facp*(1-albetp)+Kkp)/Qpsq\n", " Bx = (fac*(1-albet)+Kk)/Qsq+Bxp\n", " Byp = gamp*(facp*(1+albetp)-Kkp)/Qpsq\n", " By = gamma*(fac*(1+albet)-Kk)/Qsq-Byp \n", " return 100*Bx,100*By\n", "\n", "i = 0\n", "j = 1\n", "x = np.arange(0,1.0,0.01)\n", "y = np.arange(-0.5,0.5,0.01) \n", "X,Y = np.meshgrid(x,y) \n", "fig = plt.figure()\n", "ax = fig.add_subplot(111) \n", "Bx,By = Bfield(dx,X,Y)\n", "ax.streamplot(x,y,Bx,By)\n", "ax.set_aspect('equal')\n", "ax.set_title('Magnetic field lines')\n", "ax.set_xlabel('x')\n", "ax.set_ylabel('y')\n", "\n", "def Euler(dx):\n", " q = 1\n", " m = 10\n", " dt = 0.1\n", " v = [0.1,0.1]\n", " r = [0.6,0.8]\n", " i = 0\n", " for t in range(0,10):\n", " ii = int(r[0]/dx) \n", " jj = int(r[1]/dx)\n", " bx = (Bx[ii,jj]+ Bx[ii+1,jj])/2\n", " by = (Bx[ii,jj]+ Bx[ii,jj+1])/2\n", " B = [bx,by]\n", " F = q*(np.cross(v,B))\n", " a = F/m\n", " v = v+a*dt\n", " r = r+v*dt\n", " print(r)\n", " xx[i] = r[0]\n", " yy[i] = r[1]\n", " print(xx[i],yy[i])\n", " i = i+1\n", " \n", "Euler(dx)\n", "print(\"finished\")\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 2 }