# Source code :: MagFieldCalc

[Return]
 ``` 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``` ```#!/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-sqrt(DRc**2+zc**2) n = size(rho,0) if ao/a <=5./n: #plasma escaped out of the chamber DRc/= sqrt(DRc**2+zc**2)*(1-5./n)*a zc/= sqrt(DRc**2+zc**2)*(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 = r/rho beta_sq = (x/rho)**2 gamma = x/r Q = (1+alpha)**2+beta_sq k_sq = 4*alpha/Q 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 B[...,1] = B0/sqrt(Q)*(E_val*(1-alpha**2-beta_sq)/(Q-4*alpha)+K_val)# B_z B[...,0] = B0*gamma/sqrt(Q)*(E_val*(Q-2*alpha)/(Q-4*alpha)-K_val)# B_R return B ```