dodal tokamak

boris.py

@@ -2,8 +2,8 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib import rc
-rc('font',**{'family':'serif','serif':['Computer Modern']})
-rc('text', usetex=True)
+#rc('font',**{'family':'serif','serif':['Computer Modern']})
+#rc('text', usetex=True)
 from scipy.special import ellipe,ellipk
 
 e = 1.602176565e-19 #osnovni naboj [C]
@@ -60,7 +60,7 @@ def B_loop(p, a):
     mu = 1e-7 # mu/4/pi
     x,y,z = p
     r = np.sqrt(x**2+y**2)
-    R0 = np.sqrt((r+a)**2+z**2))
+    R0 = np.sqrt((r+a)**2+z**2)
     m = 4*a*r/( (r+a)**2+z**2)
     K = ellipk(m)
     E = ellipe(m)
@@ -79,6 +79,48 @@ def B_bottle(p, a, h, I):
     B2 = B_loop(p+np.array([0.,0.,h]),a)*I
     return B1+B2
     
+def test_transform(p, n):
+    ''' test pravilnega transforma koordinatnih sistemov'''
+    z_obrat = np.array([[1.,0.,0.],[0.,0.,1.],[0.,-1.,0.]])
+    v0 = np.array([0.,0.,.1])
+    X = np.zeros([n,3])
+    V = np.zeros([n,3])
+    U = np.zeros([n,3])
+
+    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]])
+        print(x_obrat)
+        p_ = np.matmul(x_obrat, np.matmul(z_obrat,p)) - np.array([0.5,0.,0.])
+        v_ = np.matmul(x_obrat, np.matmul(z_obrat,v0))
+        u_ = np.matmul(z_obrat.T, np.matmul(x_obrat.T,v_))
+        X[i,:]=p_
+        V[i,:]=v_
+        U[i,:]=u_
+    
+    return X, V, U
+
+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])
+
+    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.])
+        b_ = I*B_loop(p_,a)
+        b = np.matmul(x_obrat, np.matmul(z_obrat,b_))
+        B = B+b
+    
+    print(B)
+    return B
 
 def Wk(m,V):
     return 0.5*m*np.sum(V**2,1)
@@ -133,3 +175,10 @@ if __name__ == "__main__": # testni del kode
     #ax.quiver(X[0:-1:20,0],X[0:-1:20,1],X[0:-1:20,2],V[0:-1:20,0],V[0:-1:20,1],V[0:-1:20,2], color='black',normalize=True,length=.5)
     # quiver res naredi grdo sliko
     plot3(X)
+
+    X, V, U = test_transform([0.2,0.,0.], 12)
+    fig = plt.figure()
+    ax = fig.add_subplot(projection='3d')
+    ax.quiver(X[:,0],X[:,1],X[:,2],V[:,0],V[:,1],V[:,2])
+    ax.quiver(X[:,0],X[:,1],X[:,2],U[:,0],U[:,1],U[:,2])
+    plt.show()

tokamak.py

@@ -0,0 +1,27 @@
+#!/usr/bin/python3
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import rc
+#rc('font',**{'family':'serif','serif':['Computer Modern']})
+#rc('text', usetex=True)
+from scipy.special import ellipe,ellipk
+from boris import *
+
+m = m_pr
+q = e
+i = 80 # tok v ovoju
+N = 15000 # stevilo ovojev
+I = i*N
+
+a = 1.5 # ovoj
+b = 0.8 + 0.5*a 
+
+v0 = np.array([-0.1,-0.15,0.])*c
+x0 = np.array([2.3, 0.,0.])
+E = lambda x: np.array([0.,0.,0.])
+B = lambda x: tokamak(x, a, b, I, 16)
+dt = 2e-11
+tdur = 1e-7
+print(tdur/dt)
+
+X,V = boris(x0, v0, E, B, dt, tdur, q, m)