popravil polje zanke in singularnosti v osi

boris.py

@@ -57,16 +57,20 @@ def B_zemlja(x):
 def B_loop(p, a):
     ''' polje magnetne zanke polmera a v xy ravnini, normiran tok 1A
     glej: Sinigoj, ELMG polje, stran 200'''
-    mu = 1e-7 # mu/4/pi
+    mu = 2e-7 # mu/2/pi
     x,y,z = p
     r = np.sqrt(x**2+y**2)
-    R0 = np.sqrt((r+a)**2+z**2)
+    if r < 1e-15:
+        Bz = 2*np.pi*1e-7*a**2/(a**2+z**2)**1.5
+        return np.array([0.,0.,Bz])
     m = 4*a*r/( (r+a)**2+z**2)
     K = ellipk(m)
     E = ellipe(m)
 
-    Br = -mu * z/r/R0 * (2*K - (2-m)*E/(1-m))
-    Bz = mu/R0 * (2*K - (2-m*(r+a)/r)*E/(1-m))
+    C = 1./np.sqrt((a+r)**2+z**2)
+
+    Br = -mu * z/r*C * (-K + (a**2+r**2+z**2)*E/((a-r)**2+z**2))
+    Bz = mu*C * (K - (a**2-r**2-z**2)*E/((a-r)**2+z**2))
 
     return np.array([Br*x/r, Br*y/r, Bz])
 
@@ -106,20 +110,19 @@ def tokamak(p, a, b, I, n):
     ''' polje v modelu tokamaka, zanke polmera a
     nataknjene na toroid polmera b '''
     z_obrat = np.array([[1.,0.,0.],[0.,0.,1.],[0.,-1.,0.]])
-    v0 = np.array([0.,0.,.1])
-    B = np.zeros([1,3])
+    B = np.zeros(3)
+    premik = np.array([b,0.,0.])
 
     for i in range(n):
         theta = i*2.*np.pi/n
         x_obrat = np.array([[np.cos(theta), np.sin(theta),0.],
                 [-np.sin(theta),np.cos(theta), 0.],
                 [0.,0.,0.2]])
-        p_ = np.matmul(x_obrat, np.matmul(z_obrat,p)) - np.array([b,0.,0.])
+        p_ = np.matmul(x_obrat, np.matmul(z_obrat,p)) - premik 
         b_ = I*B_loop(p_,a)
-        b = np.matmul(x_obrat, np.matmul(z_obrat,b_))
+        b = np.matmul(z_obrat.T, np.matmul(x_obrat.T,b_))
         B = B+b
     
-    print(B)
     return B
 
 def Wk(m,V):

tokamak.py

@@ -17,11 +17,15 @@ a = 1.5 # ovoj
 b = 0.8 + 0.5*a 
 
 v0 = np.array([-0.1,-0.15,0.])*c
+v0 = np.array([0.8, -0.48, 0.3595])*c
 x0 = np.array([2.3, 0.,0.])
+x0 = np.array([0.,0.,0.8])
 E = lambda x: np.array([0.,0.,0.])
 B = lambda x: tokamak(x, a, b, I, 16)
+B = lambda x: B_bottle(x, 1., 5., I)
 dt = 2e-11
-tdur = 1e-7
-print(tdur/dt)
+tdur = 1e-6
+#print(tdur/dt)
 
 X,V = boris(x0, v0, E, B, dt, tdur, q, m)
+plot3(X)