From 9515e508592e3f5bc0461ac0c4286ec6cf0b3fa3 Mon Sep 17 00:00:00 2001 From: Andrej Date: Mon, 14 Jun 2021 14:56:32 +0200 Subject: [PATCH] dodal tokamak --- __pycache__/boris.cpython-37.pyc | Bin 3380 -> 5789 bytes boris.py | 55 +++++++++++++++++++++++++++++-- tokamak.py | 27 +++++++++++++++ 3 files changed, 79 insertions(+), 3 deletions(-) create mode 100644 tokamak.py diff --git a/__pycache__/boris.cpython-37.pyc b/__pycache__/boris.cpython-37.pyc index 29e0a13b4b15266d2658c87a594d391f6ec56a73..3663d1880b8e9d93f39f9c7c846b12979d5fe176 100644 GIT binary patch literal 5789 zcmb7ITZ|k>7453-p6QvLeXn=@iruz1k4ZN6#+E{2Og33!k~n061d=eZnefBqiI6Ry;aH$vGd;Oy+!}v2*b}rI9i97lWGYnw})-Xy;KY7zE zS&Z6v!)|6u8C^CT*``x+nz>RAb*q#Yc4ijJx?}(GaMUTVclztxZ+2&DkTAUGQs^%6e9nE& zxbJ4>DleTVo%9d+CvGxk_!Hv(w7+=ZCckQuuH^dG=BX{CblN{%J0ptX!CR~}>(6f9 z>z@!zus69dR`?NPP;>Ux9@4dw`z&{AspU>hTW(=(|Ie}AePi#`(y@1HyJK6M_iY*Q z7kkMlp4}$jFBO^MvQp$K<2M>X+wVW?)vx^gpEK(F*-Oc41Dh`c-YZ7-3_vdj(??c2!Vkg!y<0k7lG27&c12^(5*30dy*k@ZvhVcHvlmXok#>pom!WQ(wBZVH5}2`Z7jaS!>nN-J5gkfh)Y4@7BCJEHG}P=DC+2 z=oGc>;-}j#?%0a>qFeE>0Ik5ShtjLN+a@qj#f|IE!Zyc4+7L^{j54tdinrRz?)sRU zvZdFm`pOhxtn3O}qT;0Doq5IA=M}rG*rMVhR#q(5q+-_DSJ!c=hO9yWvGm?_3nL zsty0)wkPADfIsavx9TFg{OBWgR?Q~5CUfn{1ImYgyGVW)e=ZHB^cxLv>5s0A&8 zZPTl^{1)JbataKuEyH>gGy{0c?YbHx+A3}ZO`sW8VbiT1J@e~--8gJ41h-%G44K8J z$_&Idt$?KE&@d=klGs97(RCRsGg@mZ+mq7UDBAJ}Dfr@iiNDxY3?8p**QMeVDV@c# zo~m~awMfS_ll(BnrkQ*YZ56^BUBkHP^?Kd^iUQUhYkk%1nLVp#)9Pb|`J|CpBxm4( zUPiZO5;g-mg2ontleGF&ht}u5_xpJR-GuM z(e$O~dhWX0*>I(YjoS+5+*T-4|Hfh6edHt6^=Gu=vJZyc^^b^Sh2QZp5o25)a7hLv95do zI(U9w9zsE}=jE7gSd=Bc*gXh@%wo9_hV92;5+R{648R}{5FBMLbC|==vUx7gqrE~e zr;ETqpwxhP8ig2}f+7-y(S*~&VX4_oPRrWEmOH4O+0`PxCH4$$A>MuT;6g{SUvK!W zpyIph0l39?&#!|2cb7h0T+BdUMN`i+T`s*o894`ye3-;>2*sYmu6P9Z`Gd=V{jeOy zal?Neb`gD|VV+P^NcE%6oR=Y!Uxu=t=M0DG$-?w01DN zY=2g{IGf}sUxS5{ZdeQA>N~YPn=cw^P`+Fw>mK z4FbQtuU^28+@rs~fJD&QtYPyLA~m^v_Ya3>?#j&F_a=`d7#Tx zHwu7{X1LyukRL4c3xk;X!8_HrZ`1Y51T84`zS1okz>0%game(b4i zCpE&jx`??4X#EfscE_gp$)N zC75w8$r8qsERdD!zR;(a^)ynixgcn->D7_N9bj}#D?pi7Z`FLTP75FmWf+J7^DMjn z6Y6xpCNB9n^gbDI@-r}&kC7lkmY*d-i^x>xYOtEEoffS!3@dOEtym2Z$N#ZFI zizJAs6<;nNL|V4L0s-_9@uOjCi0y&bEFd<;3fy7x8R!Q@NXY?xosIznbTB^xFaZGM zhLA5n6goAsmiixmXCX~xs9TYz+wI=8OmRwYmxD9puhe_cMLHK@@>#MY<1bVw2+~D+ zHH%wMWTel?TLv;fvqk}l-vY}fT;eFSX$jNK6}55X(y&%`xE5ph0lJkHw)_bBYtZ#ofx4l)T!Fo^sHXCf zAD24~Z`p6)EffLzZ5wHu)_2pg6|IKf#HJ$VRkq&>W51M1GhAh2hKhR~c>+p#lmsQn zl2TtePr@AtO~E$3{=&D=m~IQOY=K$W2xDN{!^jFI%tukikIOQwaMIcsur4L#14RB7 z71|U)q@CEC9M4RW*`(8jMgjPZ9nt@O>}i4Twg1g66W-2jTI5aWf#W^a_HGSPI@{Ox z`=UnrehqqHc@sR2EWxBPfFWs#TSN zkgt&XRT5vP#%$!R`*5Iq4Ygt}4Xsy6{W_^R9e!6@9&+(()c6hEC^|UqZhXwGlKlEz zSVYqh2IaTxEN+*%<}|{o1-6q7Xp~bfH@%=$E)R~Mpq5?6qlB8=s#8Cm=wYi3m=1iO zCLYKzF6irb3pw|3n(6V&0>u%KP`>9f$U{2_YC5TZ(b(fgf#VYg+WdSmLR&-|(T#=3Y zN#Q^ZN0p$xaW-oE^#6~ffS}Juih)0x_!pxD3FcrG{Xc>})TXhNQhG_5l8$L!3V>?jA6@RPDPQ z{~D|{{3qy(qDv6OGzD{vTL2?|Jb%i`I}=XUnQ^8ZGe0>!o+sGiKN8sRY0Kf$s2A8( G{Qd)+NHc!` delta 1189 zcmah|&2Jl35Z`&;b=HpU_}fVusGZO>o0bUF9zc{3(vqSAJ%o?YRV8@0SSczyL!K!nfIG{^JZqheD<69mCbxU zYw`1N{<}w+oMrte%0JH(3md$mpRLbU7LRorub_=i7cl*W`UJpY%M`kp`B;I4S$_(r zeHU|n2J?OvmvIJ*C*V(G3CkyzpTi1P`J2b9*uYudDPRpPEOZslA#`k<$4hLT(M4R? zSM&oug--fy$^h4!0 z{8?SljmMlY!rPX&awZqN?}1UYXY5Wm92+On)axA0+ztn$$yihKtq?U0J!R6JV9@Iy zl3isdV~_R0_#JjlPa@6o>6+SB;Km=>FM}U{Yrp9FYSRANqSoim9mvP8IUX#XhK>yj zQfJCAoH-(afXm?7l#4%2&u&~|#aElp8MS?6fDG(YOGpVUvSi%N_K_a+KM2-Xi;~8r zY_L@o>Oj}xKMR{%$tg)#Mchw|mKi)n4Uy*sF2qzj`!BFLQVbT*b)hs)PiF1{WaICP zvuO!-X2-u2*Ow%qar?btGz#{#QKO!5LMfL??S+v^MZIDCX6fD6=Y;3zVK~+%eV(yN zb$VT+f&*&uCM^iag$*1UIJ%gWGwO%q8!yGbm9D_OxKjQj$%0nnO6B?E8{&IYeBEfW zE6j-2MYX%sLwZHFHw0c4<5U!UsQGNV&6<}X$&Ky^ydm3)^L9GJV4!JBjCT_w9vFp( z|H+Evu-1L1M0ExWTz=PhHKCcRa}#VR0PVy-Ro=gu6#pudJV-8jF_?R+9}IReSbxNv z-+5W|ARbk>SO3S69xyc!x*z{qy}r!hQ><}|rBUY-b>8kRZ_YTw$zXKI1NL6+a#;#z zvV&kO*XZ|lNySfU_2VV+yC#q{|DH%YLUr41DNMUf6|u<=Kyo8036uqD0!c|P%+|k# u_SPW8Nnfwaghi_iF68(vLr(GV!h%e2w&iBsRJKwtNl9H*grc%l-Tn)shWqgV 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)