#!/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
