diff --git a/__pycache__/boris.cpython-37.pyc b/__pycache__/boris.cpython-37.pyc index 29e0a13..3663d18 100644 Binary files a/__pycache__/boris.cpython-37.pyc and b/__pycache__/boris.cpython-37.pyc differ diff --git a/boris.py b/boris.py index a9b5418..7a3dd98 100644 --- a/boris.py +++ b/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() diff --git a/tokamak.py b/tokamak.py new file mode 100644 index 0000000..a920c33 --- /dev/null +++ b/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)