diff --git a/__pycache__/boris.cpython-37.pyc b/__pycache__/boris.cpython-37.pyc index 0b74ed8..f2fbb3e 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 34f84df..9c1e87f 100644 --- a/boris.py +++ b/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) diff --git a/radii.pdf b/radii.pdf index c1d8a3b..7a53ad3 100644 Binary files a/radii.pdf and b/radii.pdf differ diff --git a/tokamak.py b/tokamak.py index 298d07e..9ef0065 100644 --- a/tokamak.py +++ b/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)