dodal relativistiko. neuspesno

boris.py

@@ -27,23 +27,44 @@ def boris(x0, v0, E, B, dt, tdur, q=1, m=1):
     V = np.zeros((duration,3))
     x = x0
     v = v0
+    u = v0 / np.sqrt(1-(np.sum(v0**2)/c**2))
 
     q_prime = dt * q * 0.5 / m
     
     # https://en.wikipedia.org/wiki/Particle-in-cell
+    # https://warpx.readthedocs.io/en/latest/theory/picsar_theory.html
     for time in range(duration):
-        h = q_prime*B(x)
-        s = 2*h / (1+np.dot(h,h))
-        u = v + q_prime*E(x)
-        u_ = u + np.cross(u+np.cross(u,h),s)
-        
-        v = u_ + q_prime * E(x)
+        u_min = u + q_prime*E(x)
+        t = q_prime*B(x)/np.sqrt(1+np.sum(u_min**2)/c**2)
+        u_prime = u_min + np.cross(u_min,t)
+        u_plu = u_min + np.cross(u_prime,2*t)/(1+np.sum(t**2))
+        u = u_plu + q_prime*E(x)
+        v = u/np.sqrt(1+np.sum(u_min**2)/c**2)
         x = x + dt*v
         V[time,:] = v
         X[time,:] = x
 
+
+        #implementacija po wikipediji - nerelativisticno
+        #h = q_prime*B(x)
+        #s = 2*h / (1+np.dot(h,h))
+        #u = v + q_prime*E(x)
+        #u_ = u + np.cross(u+np.cross(u,h),s)
+        #
+        #v = u_ + q_prime * E(x)
+        #x = x + dt*v
+        #V[time,:] = v
+        #X[time,:] = x
+
     return X, V
 
+def vay(x0, v0, E, B, dt, tdur, q=1, m=1):
+    ''' doi: 19,1963/1.2837054 '''
+    print("ni implementiran")
+    return 0
+
+
+
 def B_zemlja(x):
     ''' model zemljinega magnentnega polja aproksimiran z 
     magnetnim dipolom '''
@@ -64,8 +85,8 @@ def B_loop(p, a):
         Bz = 2*np.pi*1e-7*a**2/(a**2+z**2)**1.5
         return np.array([0.,0.,Bz])
     m = np.sqrt(4*a*r/( (r+a)**2+z**2))
-    K = ellipk(m)
-    E = ellipe(m)
+    K = ellipk(m**2)
+    E = ellipe(m**2)
 
     C = 1./np.sqrt((a+r)**2+z**2)
 

tokamak.py

@@ -7,25 +7,26 @@ from matplotlib import rc
 from scipy.special import ellipe,ellipk
 from boris import *
 
-m = m_pr
+m = m_el
 q = e
 i = 80 # tok v ovoju
-N = 150000 # stevilo ovojev
+N = 1500 # stevilo ovojev
 I = i*N
 
 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*0.2
+v0 = np.array([0.8, -0.48, 0.3595])*c
 x0 = np.array([2., 0.,0.])
-#x0 = np.array([0.,0.5,-0.8])
+x0 = np.array([0.,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) + B_loop(x,4.)*100*8000
-dt = 2e-12
+B = lambda x: B_bottle(x, 1., 5., I) + B_loop(x,4.)*100*8000
+dt = 1.5e-12
 tdur = 3e-7
-#print(tdur/dt)
+#dt = 2e-11
+#tdur = 1e-6
 
 X,V = boris(x0, v0, E, B, dt, tdur, q, m)
 plot3(X)