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
|
