From f94f50b282555feaabd046d694328bdbce3750db Mon Sep 17 00:00:00 2001
From: Andrej <andrej@baza>
Date: Wed, 16 Jun 2021 09:34:24 +0200
Subject: [PATCH] popravil polje zanke in singularnosti v osi

---
 __pycache__/boris.cpython-37.pyc | Bin 5789 -> 5914 bytes
 boris.py                         |  21 ++++++++++++---------
 tokamak.py                       |   8 ++++++--
 3 files changed, 18 insertions(+), 11 deletions(-)

diff --git a/__pycache__/boris.cpython-37.pyc b/__pycache__/boris.cpython-37.pyc
index 3663d1880b8e9d93f39f9c7c846b12979d5fe176..47de2a90832305051352e114c05d8c9565806d55 100644
GIT binary patch
delta 1309
zcmZWnO=w(I6u#%ad4Dr+rb#CAqlI)*O<|-nTCJLDOziZBq+L|4woPblQeQGmn#oMh
zdkMANN0hkg_JT-3aiE~CiUdj#_wL;YZ{tQ3K|#_@K^LC$W+J7|@V@i+oqNyE{k(Rk
zo@->YX+ghZ_`~8_;oID<B`<-N_Mt6h`NnmDE4~cHsvH%5YE&rjWnSz_WwKQgzL^*M
zMP!HKOT)K(OW94M2j_)vD<H1Lp910$doPFhPMA{8Hbib{CNk|ibOI<)7IC*8Blm#l
z86`86&WG(6DnrsGJ1K}tO$ymNT%{>W$CRJAvP7jvo|66S$8*br!blEsC!b6%b8^z7
z@{}<s!&5vJE@;qm)z!Op=e^t0Xwy`i#S6cEanJa1syTb_`#ZZo-I{Llu!VRuGf;&1
zJw^7w`fH|Yp-VPCOAyO$3lFiy6g2d0XZ)SUEi7;h&f-uaPb|m=VucwAjI-BjTwUvQ
zy2nTlo+c1R3er#n2cCl=sE~Bzn4E<viAVJ#=d03D;`DFk3#p14JL>wrdwBR9sewfx
zj4>DkoZNZn{tXj&m^`uZK_`ffptT;GU2Ls~vDID*+Z$^!yj#s;O^<ksa6t!N4NmEv
z_muoq=v?}lnn|mceyRtNl`LKQ^yg<2**EWohDvS2Npa<~$We}O>Y?oMd8=t8@ig&m
z)kjl_p6p3wm-&4dp`$Viw&}4(Ucxe0-9CHORUR`VNl#IhL}cAKeUTP-UQ})UU3OK!
zn5)loKgc7)*BGW4lDWJ=8aKm^IzbHA@MSip8D<!U2^Op_F1K6P@i==L3?~>)>if9^
z6KzrlyEsU@7?L>{fimqRN64QiaaRA6+suB%Oi3UB=XSRCJ%Ev9)MjTRT*g6#dT#iW
z$_vE7BMg=QAlbj!UF(D|ZR<aWCtz_WKXL@Nc#HTJLF@z@*Ek;FMK)7SZ25SJ!}AOa
z?6iZW&E=I=6EBlhb+N>eD-25<mb7$pWqql=jPJAaYU1>7#?t@rf3+ZAHvTD5f+B(7
z*9Ud_Do`<sG=UWLk8S-+{^QKo6o1Gh?C9<&raPp5!)90iG5Wf5hSVzHS^a$BG@R3$
zg);||4NWrT6vHgT9K#|*TmM{m<=~sFa_{|W-)1c-*}AdO+N7CY)~^<4?4n~kB`2?=
GV&z}DIS+vV

delta 1035
zcmZvZTTc@~6vyYxUc2oUs0-y5iil_vK){F>A_x*eh)5t&xl>XgB|u?Xh>&U2;0GX?
z_+nxr4{D;sXz;;Lpx<G?fr-Az8}i^e3oj;av%mi5f6m#NIp<~btEOm=>pBb{|K9i6
zWAAhHx5<Ap#jS?wJbMk&ln`fWeuhbd7-xG>u+o}Zwh~Xcw0EEo%4@;DBY7CJJSNR5
zCao&0k|B?1mow?0(n(dZa^|r@Sn_{!PIfw}{!e}cPovAQ)K0%2mU{82j4Wz<Fv^mK
znzpm4?M`9gb$;jwVIY!*2$8NQ1g=J^fkLhiYd)lWzOo@I(C}f)hi!2Y<v1T{5v74Q
zD2fB4M~w5HZ8a9EuJq7gnhp+>K_%3J4K{Cv7B1S<i_p7x8%9rd4Yn~wBN<I<)M+dk
zM42-YOsjz!5bX!M&L7ZyzLrtDkpwiVPmx&|Q6+8@^r|kmf*&#U!8zAz&_L20e>&Z`
zVB(smjkLMQX(*5PaS-|r^s=ef1zXx_OTBly!B#)r%AsJERLZHeocgcjZR{uPaGVBp
z*OE>^Lq_)a9xvd0B#*O#wLb8W<wMfi0l>LpqeH4Ux}*jpj!ykmAEV6=uA&gv2-gW6
zgrEsGP$k{d`=j<@PF$wwO+qK3h)_bfMd&8<5N@k0MTrX;l(}jI)9RoYqC5&QD964>
zp&nlkBHRYir?!eV-2&xu$_3GXa9s2iVstX!$Zq84#pwanx02`PCA3A08je-KQ?(R3
z4>Je5u~wKKBL6<s*tw0lwdHI+CmzsJUlf_Ncu3(#gfTL$TxxTEA+snZaFq;+I7KE2
zvk0GObH1^#np&9`Q)Eu7gm)H{>h=<_tWw^y@Ef!b@nvFF{qlsfit;Tjvno^CVc$iW
z1TmoYOYgvtIxHPX1YK>XlwLv~p`S2ASXPr|?R8xw=@KW~4wITj_~!aXW)mARp}v%L
PSuxwP<F==c%4+`ti%QnS

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)