1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133 | #!/usr/bin/env python
# -*- coding: utf-8 -*-
from numpy import *
#from matplotlib.pyplot import *
from numpy.linalg import norm
from scipy.special import ellipk, ellipe
#import time
# functions for the calculation of the magnetic field around the current loop, by function Bloop_integ
#is values calculated by integration, in the function Bloop_analytic is used analytical formula
#Autor :Tomas Odstrcil
# obecný pruodový profil může být nahrazen skupinou bodých proudových profilů v těžištích příslušných mezikruží nesoucích proud
global mu_0,a
a = 0.085 #[m]
R_0 = 0.4 #[m]
mu_0 = 4*pi*1e-7
p = 2 #píkovací faktor
def CurrentProfile(rho,zeta,Rc,zc,I, Shaf_shift, centered = 'False'):
#parabolický profil, rho,zeta jsou souřadnice, Rc,zeta - střed plazmatu (ale ne proudového profilu!),
#Shaf_shift veličina udávají shafranův posuv většinou je 0-0.2
j= zeros(shape(rho)+shape(Rc)) #předpokládám že jeden z nich je číslo
DRc = Rc-R_0
Drho = rho-R_0
ao = a-hypot(DRc,zc)
n = size(rho,0)
if ao/a <=5./n: #plasma escaped out of the chamber
DRc/= hypot(DRc,zc)*(1-5./n)*a
zc/= hypot(DRc,zc)*(1-5./n)*a
ao = 5*a/n
if centered:
r2 = (zeta)**2+(Drho)**2 #centered profile
else:
r2 = (zeta-zc)**2+(Drho-DRc)**2 #shifted profile
plazma = (r2 < ao**2)
j[plazma] = (p+1)*I/(pi*ao**2)*(1-r2[plazma]/ao**2)**p
j[plazma]*= (1+(j[plazma])/amax(j[plazma])*Shaf_shift*sin(pi*(rho[plazma]-Rc)/ao))
return j,ao
def Bloop_integ(R,z,rho,zeta): # toroidální vodič - vypočtenu přímou integrací
#R,Z poloha detektoru, rho, zeta poloha vlákna
m = 1000
dphi = pi/m #protože jsem integroval jen přes půlkruh
phi_i = (arange(0,m))*dphi
B = zeros(3)
R_Q = zeros(3)
R_Q[2] = zeta
dR_Q = zeros(3)
dR = zeros(3)
R_D = zeros(3)
R_D[0] = R
R_D[2] = z
S = sin(phi_i)
C = cos(phi_i)
N = arange(m)
Bvec= empty((m,3))
for i,c,s in zip(N,C,S):
R_Q[0] = rho*c
R_Q[1] = rho*s
dR_Q[0] = -s
dR_Q[1] = c
dR = R_D-R_Q
dB = cross(dR_Q, dR)/norm(dR)**3
Bvec[i,:] = dB
if i == 0 or i == m: #simpson integration formula
B+= dB
elif i%2 == 0:
B+= 2*dB
else:
B+= 4*dB
B*=mu_0/(4*pi)*rho*(dphi/3)
B*= 2 #protože jsme integroval jen přes půlkruh
B[1] = 0#protože jsme integroval jen přes půlkruh
Br = B[0]
Bz = B[2]
return Bz, Br
def Bloop_analytic(R,z,rho,zeta,out=None):
#R,Z poloha detektoru, rho, zeta poloha vlákna
#http://www.netdenizen.com/emagnettest/offaxis/?offaxisloop
r = abs(R)
x = z-zeta
B0 = mu_0/(2.*pi*rho)
alpha = divide(r,rho)
gamma = divide(x,r,r)
beta_sq = square(divide(x,rho,x),x)
Q = square(alpha+1.)+beta_sq
sQ = sqrt(Q)
k_sq = divide(alpha,Q)*4.
E_val = ellipe(k_sq)
K_val = ellipk(k_sq)
if out == None:
B = empty(shape(rho)+shape(R)+ (2,)) # předpokládám že jedno z nich je pole a druhé konstanta
else:
B = out
Bz = (E_val*(1.-square(alpha)-beta_sq)/(Q-alpha*4.)+K_val)*B0/sQ# B_z
Br = B0*gamma/sQ*(E_val*(Q-alpha*2.)/(Q-alpha*4.)-K_val)
B[...,1] = Bz
B[...,0] = Br
return B
|
