naredil uporaben primer magnetne boce, popravil transformacije koordinatnih sistemov in naredil izris tokamaka

boris.py

@@ -104,44 +104,23 @@ 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.]])
+    Rx = np.array([[1.,0.,0.],[0.,0.,-1.],[0.,1.,0.]])
     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.],
+        Rz = np.array([[np.cos(theta), np.sin(theta),0.],
                 [-np.sin(theta),np.cos(theta), 0.],
-                [0.,0.,0.2]])
+                [0.,0.,1.]])
-        p_ = np.matmul(x_obrat, np.matmul(z_obrat,p)) - premik 
+        #print(Rz)
+        p_ = np.matmul(Rx.T, np.matmul(Rz.T,p)) + premik 
         b_ = I*B_loop(p_,a)
-        b = np.matmul(z_obrat.T, np.matmul(x_obrat.T,b_))
+        b = np.matmul(Rz, np.matmul(Rx,b_))
         B = B+b
     
     return B
@@ -149,7 +128,7 @@ def tokamak(p, a, b, I, n):
 def Wk(m,V):
     return 0.5*m*np.sum(V**2,1)
 
-def plot3(X,enote='m',zemlja=False):
+def plot3(X,enote='m',zemlja=False,lim=[],bottle=[],tokamak=[]):
     fig = plt.figure()
     ax = fig.add_subplot(projection='3d')
     if zemlja:
@@ -165,7 +144,50 @@ def plot3(X,enote='m',zemlja=False):
     ax.set_xlabel(r'$x$ [{}]'.format(enote))
     ax.set_ylabel(r'$y$ [{}]'.format(enote))
     ax.set_zlabel(r'$z$ [{}]'.format(enote))
+    if lim:
+        ax.set_xlim(lim)
+        ax.set_ylim(lim)
+        ax.set_zlim(lim)
 
+    if bottle:
+        r = bottle[0]
+        h = bottle[1]
+        fi = np.linspace(0,2*np.pi,300)
+        a = r*np.cos(fi)
+        b = r*np.sin(fi)
+        ax.plot(a,b,+h,color='red',zorder=100)
+        ax.plot(a,b,-h,color='red',zorder=100)
+    if tokamak:
+        ax=plot_tokamak(tokamak[0],tokamak[1],tokamak[2],ax)
+
+
+    plt.show()
+
+def plot_tokamak(a,b,n,ax=None):
+    fi = np.linspace(0,2.*np.pi,300)
+    x = a*np.cos(fi)
+    y = a*np.sin(fi)
+    z = np.zeros(300)
+    X=np.zeros([300,3,n])
+
+    x = x-b
+    Rx = np.array([[1.,0.,0.],[0.,0.,-1.],[0.,1.,0.]])
+
+    if ax == None:
+        fig = plt.figure()
+        ax = fig.add_subplot(projection='3d')
+
+    for i in range(n):
+        theta = i*2.*np.pi/n
+        Rz = np.array([[np.cos(theta),-np.sin(theta),0.],
+            [np.sin(theta),np.cos(theta),0.],
+            [0.,0.,1]])
+        for j in range(300):
+            p = np.array([x[j],y[j],z[j]])
+            #X=np.zeros([300,3])
+            X[j,:,i] = np.matmul(Rz, np.matmul(Rx,p.T))
+        ax.plot(X[:,0,i],X[:,1,i],X[:,2,i],color='red',zorder=1)
+        plt.draw()
     plt.show()
 
 def rplot(X, enote='m', xy=[0,1],os='xy', fname=None):
@@ -194,15 +216,8 @@ if __name__ == "__main__": # testni del kode
 
     dt = 1e-2
 
-    X, V = boris(x0,v0,E,B,dt,20)
+    #X, V = boris(x0,v0,E,B,dt,20)
 
-    #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)
+    #plot3(X)
-    # quiver res naredi grdo sliko
+    plot_tokamak(1,2,16)
-    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()

majhenprimer.py

@@ -7,19 +7,21 @@ from matplotlib import rc
 from scipy.special import ellipe,ellipk
 from boris import *
 
-m = m_el
+m = m_pr
 q = e
-i = 5 # tok v ovoju
+i = 10 #tok v ovoju
 
-a = 0.5 # ovoj
+a = 0.1 # ovoj
 
-v0 = np.array([0.0,-0.2,0.1])
+v0 = np.array([0.0,-3,1.9])
-x0 = np.array([.1, 0.,0.])
+x0 = np.array([.05, 0.,0.02])
 E = lambda x: np.array([0.,0.,0.])
 B = lambda x: B_bottle(x, a, a/2, i)
-dt = 1e-2
+dt = 1e-5
-tdur = 10
+tdur = 1e0
 #print(tdur/dt)
 
 X,V = boris(x0, v0, E, B, dt, tdur, q, m)
+plot3(X, lim=[-1.2*a,1.2*a],bottle=[a,a/2])
+X,V = boris(x0, v0, E, B, dt, tdur/10, q, m)
 plot3(X)

tokamak.py

@@ -7,26 +7,24 @@ from matplotlib import rc
 from scipy.special import ellipe,ellipk
 from boris import *
 
-m = m_el
+m = m_pr
 q = e
-i = 80 # tok v ovoju
+i = 800 # tok v ovoju
-N = 1500 # stevilo ovojev
+N = 1000 # stevilo ovojev
 I = i*N
 
-a = 1.5 # ovoj
+a = .8 # ovoj
-b = 0.8 + 0.5*a 
+b = 0.5 + 0.5*a 
 
-v0 = np.array([-0.1,-0.15,0.])*c
+v0 = np.array([-0.1,0.15,+0.])*c/2/2
-v0 = np.array([0.8, -0.48, 0.3595])*c
+x0 = np.array([1., 0.,0.4])
-x0 = np.array([2., 0.,0.])
-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 = 1.5e-8
-dt = 1.5e-12
+dt = 1.5e-4
-tdur = 3e-7
+tdur = 3e0
-#dt = 2e-11
+dt = 2e-11
-#tdur = 1e-6
+tdur = 1e-6
 
 X,V = boris(x0, v0, E, B, dt, tdur, q, m)
-plot3(X)
+plot3(X,lim=[-b-a,b+a],tokamak=[a,b,16])