diff --git a/__pycache__/boris.cpython-37.pyc b/__pycache__/boris.cpython-37.pyc index 3663d18..47de2a9 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 7a3dd98..767fb82 100644 --- a/boris.py +++ b/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): diff --git a/tokamak.py b/tokamak.py index a920c33..d435e70 100644 --- a/tokamak.py +++ b/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)