Particle-in-cell simulacija trajektorije nabitega delca.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

#### 223 lines 6.1 KiB Raw Permalink Blame History

 `#!/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` ``` ``` `e = 1.602176565e-19 #osnovni naboj [C]` `m_pr = 1.672621777e-27 #masa protona [kg]` `m_el = 9.10938291e-31 #masa elektrona [kg]` `c = 299792458 #hitrost svetlobe [m/s]` ``` ``` `def boris(x0, v0, E, B, dt, tdur, q=1, m=1):` ` ''' Borisov algoritem zabjega skoka za` ` integracijo diferencialne enačbe` ` x0 = [x,y,z] : vektor zacetne pozicije` ` v0 = [vx,vy,vz] : vektor zacetne hitrosti` ` E : funkcija, ki daje jakost polja E(x)` ` B : funkcija, ki daje jakost polja B(x)` ` dt: casovni korak '''` ``` ``` ` duration = int(tdur/dt)` ``` ``` ` X = np.zeros((duration,3))` ` 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):` ` 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 '''` ` B0 = 3.07e-5 # [T]` ` Rz = 6378137 # radii zemlje [m]` ``` ``` ` r = np.sqrt(np.sum(x**2))` ` k = -B0*Rz**3/r**5` ` return k*np.array([3*x[0]*x[2],3*x[1]*x[2],2*x[2]**2-x[0]**2-x[1]**2])` ``` ``` `def B_loop(p, a):` ` ''' polje magnetne zanke polmera a v xy ravnini, normiran tok 1A` ` glej: Sinigoj, ELMG polje, stran 200'''` ` mu = 2e-7 # mu/2/pi` ` x,y,z = p` ` r = np.sqrt(x**2+y**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 = np.sqrt(4*a*r/( (r+a)**2+z**2))` ` K = ellipk(m**2)` ` E = ellipe(m**2)` ``` ``` ` 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])` ``` ``` ``` ``` ``` ``` `def B_bottle(p, a, h, I):` ` ''' polje magnetne steklenice visine 2h in polmera a` ` '''` ` B1 = B_loop(p-np.array([0.,0.,h]),a)*I` ` B2 = B_loop(p+np.array([0.,0.,h]),a)*I` ` return B1+B2` ` ` ``` ``` `def tokamak(p, a, b, I, n):` ` ''' polje v modelu tokamaka, zanke polmera a` ` nataknjene na toroid polmera b '''` ` 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` ` Rz = np.array([[np.cos(theta), np.sin(theta),0.],` ` [-np.sin(theta),np.cos(theta), 0.],` ` [0.,0.,1.]])` ` #print(Rz)` ` p_ = np.matmul(Rx.T, np.matmul(Rz.T,p)) + premik ` ` b_ = I*B_loop(p_,a)` ` b = np.matmul(Rz, np.matmul(Rx,b_))` ` B = B+b` ` ` ` return B` ``` ``` `def Wk(m,V):` ` return 0.5*m*np.sum(V**2,1)` ``` ``` `def plot3(X,enote='m',zemlja=False,lim=[],bottle=[],tokamak=[]):` ` fig = plt.figure()` ` ax = fig.add_subplot(projection='3d')` ` if zemlja:` ` u = np.linspace(0, 2 * np.pi, 100)` ` v = np.linspace(0, np.pi, 100)` ` x = np.outer(np.cos(u), np.sin(v))` ` y = np.outer(np.sin(u), np.sin(v))` ` z = np.outer(np.ones(np.size(u)), np.cos(v))` ` ax.plot_surface(x,y,z,zorder=-4)` ` #ax.set_box_aspect(np.ptp(X,axis=0))` ``` ``` ` ax.plot(X[:,0],X[:,1],X[:,2],zorder=4)` ` 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):` ` fig = plt.figure()` ` ax = fig.subplots()` ` ax.plot(X[:,xy[0]],X[:,xy[1]], 'k')` ` ax.grid()` ` dx = X[-1,xy[0]]-X[-2,xy[0]]` ` dy = X[-1,xy[1]]-X[-2,xy[1]]` ` #ax.arrow(X[-1,xy[0]],X[-1,xy[1]],dx,dy,head_width=0.1,head_length=0.2,` ` # overhang=0.3, zorder=10,color='k')` ` ax.set_xlabel(r'\${}\$ [{}]'.format(os[0],enote))` ` ax.set_ylabel(r'\${}\$ [{}]'.format(os[1],enote))` ` if fname:` ` plt.savefig(fname,bbox_inches='tight')` ` plt.show()` ``` ``` ``` ``` ``` ``` `if __name__ == "__main__": # testni del kode` ` E = lambda x: np.array([0.0,0.2,0.3])` ` B = lambda x: np.array([0.,0.0,1.*x[0]])` ``` ``` ` x0 = np.array([-1.,0.,0.])` ` v0 = np.array([0.,1.,0.])` ``` ``` ` dt = 1e-2` ``` ``` ` #X, V = boris(x0,v0,E,B,dt,20)` ``` ``` ` #plot3(X)` ` plot_tokamak(1,2,16)` ``` ```