function [B] = B_Ellip_10T_DNP(mode,magnet,I,r,z)
%
% Call: [B] = B_Ellip_10T_DNP(mode,magnet,r,z)
% Examples:
% 1. Plot axis magnetic field profile:
%
% >> z = linspace(0,1,101);
% >> r = 0*z + 0.0;
% >> I = [61.56 66.45 113.43 108.09 ];
% >> [B] = B_Ellip_Cryogenic_170('bz','cryogenic',I,r,z);
% >> figure(1); plot(z,B)
%
% 2. Plot field lines (r* aphi = cst)
%
% >> z = linspace(0,1,51);
% >> r = linspace(0,0.2,21);
% >> [Z,R] = meshgrid(z,r);
% >> I = [61.56 66.45 113.43 108.09 ];
% >> [Aphi] = B_Ellip_Cryogenic_170('aphi','cryogenic',I,R,Z);
% >> figure(2); contour(Z,R,R.*Aphi);
%
% Magnetic field associated with GT170 GHz magnet, with currents Igun and Icav,
% as read on power supplies (i.e. a current -(Icav+Igun) is actually flowing in the
% gun coil.
% The computation is based on the elliptic integral formulae that can be found in:
% Smythe: "Static and Dynamic Electricity", Third edition, McGraw-Hill, p.290.
%
% Version : 1.0 J.-P. Hogge Sep. 16, 2003
% 1.1 J.-P. Hogge Jan. 12, 2005 Included 165Ghz FZK magnet, other options deactivated
%
% mode : 'bz','br','dbzbz','dbzdr','dbrdz','dbrdr', 'd2bzdz2' or 'aphi'
% magnet : 'fzk_165ghz_oi' only so far
% z : a row vector (in [m]) or a 2D matrix containing the z locations where the field is to be evaluated
% r : a row vector (in [m]) or a 2D matrix containing the r locations where the field is to be evaluated
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 10T magnet for DNP gyrotron coils geometry [m]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch lower(magnet)
case 'bad' % Geometry of the 10Tesla sup.cond. solenoid with 76 mm bore diameter
Ncoil = 13;
Nturn = [1670 2117 2960 4995 12793 59 57 56 54 59 57 56 54];
rmin = [52.592 60.069 76.919 82.502 90.129 103.638 104.187 104.735 105.283 103.638 104.187 104.735 105.283].*10^(-3);
rmax = [60.069 66.051 82.502 90.129 103.638 104.187 104.735 105.283 105.832 104.190 104.735 105.283 105.832].*10^(-3);
zmin = [-145 -145 -150 -150 -150 -150 -150 -150 -150 -150 -150 -150 -150].*10^(-3);
zmax = [145 145 150 150 150 150 150 150 150 150 150 150 150].*10^(-3);
na = 8. *[10 10 10 10 10 10 10 10 10 10 10 10 10]; % Number of subcoils (radial direction
nb = 8. *[10 10 10 10 10 10 10 10 10 10 10 10 10]; % Number of subcoils in longitudinal direction
% I = 31.31 A => B~ 1.7T at flange
% I = 32.35 A => B< 3T at z = 0 => 235 mm from top flange
Current = [ I(1) I(2) I(3) I(4) I(5) I(6) I(7) I(8) I(9) I(10) I(11) I(12) I(13)]; % Current flowing in each power supply
rcen = (rmin + rmax)/2;
zcen = (zmin + zmax)/2;
dr = rmax - rmin;
dz = zmax - zmin;
coilsurf= dr .* dz;
JTot = Current .* Nturn; % Total current flowing in each coil
J = JTot./(dr.*dz); % Current density in each coil
Ic1 = JTot ; % Total current flowing in each coil [A]
case 'new' % Geometry of the 10Tesla sup.cond. solenoid with 76 mm bore diameter
Ncoil = 13;
Nturn = [1670 2117 2960 4995 12793 59 57 56 54 59 57 56 54];
rmin = [52.592 60.069 76.919 82.502 90.129 103.638 104.187 104.735 105.283 103.638 104.187 104.735 105.283].*10^(-3);
rmax = [60.069 66.051 82.502 90.129 103.638 104.187 104.735 105.283 105.832 104.190 104.735 105.283 105.832].*10^(-3);
zmin = [-144.5 -144.5 -149.5 -149.5 -149.5 -149.5 -149.5 -149.5 -149.5 112.661 113.658 114.655 115.652].*10^(-3);
zmax = [144.5 144.5 149.5 149.5 149.5 -112.661 -113.658 -114.655 -115.652 149.5 149.5 149.5 149.5].*10^(-3);
na = 8. *[10 10 10 10 10 10 10 10 10 10 10 10 10]; % Number of subcoils (radial direction
nb = 8. *[10 10 10 10 10 10 10 10 10 10 10 10 10]; % Number of subcoils in longitudinal direction
% I = 31.31 A => B~ 1.7T at flange
% I = 32.35 A => B< 3T at z = 0 => 235 mm from top flange
Current = [ I(1) I(2) I(3) I(4) I(5) I(6) I(7) I(8) I(9) I(10) I(11) I(12) I(13)]; % Current flowing in each power supply
rcen = (rmin + rmax)/2;
zcen = (zmin + zmax)/2;
dr = rmax - rmin;
dz = zmax - zmin;
coilsurf= dr .* dz;
JTot = Current .* Nturn; % Total current flowing in each coil
J = JTot./(dr.*dz); % Current density in each coil
Ic1 = JTot ; % Total current flowing in each coil [A]
otherwise
disp('Unknown magnet')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct array of current loops
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Number of coils
ncoil = length(rcen);
% Number of subcoils
ntot = sum(na.*nb);
% Preallocate space
r2 = zeros(1,ntot);
z2 = zeros(1,ntot);
Ic = zeros(ntot,1);
% isub is the index of subcoil
isub = 0;
for icoil=1:ncoil
for ib=1:nb(icoil)
for ia=1:na(icoil)
isub = isub +1;
r2(isub) = (rcen(icoil) - dr(icoil)/2) + (dr(icoil)/nb(icoil))*(ib-0.5);
z2(isub) = (zcen(icoil) - dz(icoil)/2) + (dz(icoil)/na(icoil))*(ia-0.5);
Ic(isub) = Ic1(icoil) / (na(icoil)*nb(icoil));
end
end
end
% printf('isub =',isub)
% Transform array of evaluation points in a one-dimensional array
%
r1d = reshape(r,1,prod(size(r)));
z1d = reshape(z,1,prod(size(z)));
%
% Suppress locations where Green function diverges (i.e. on the sources)
%
for i=1:length(r2),
[I] = find((r1d.*r1d -r2(i)*r2(i) == 0) & (z1d.*z1d - z2(i)^2)==0);
if ~isempty(I)
r1d(I)=NaN;
z1d(I)=NaN;
end
end
% Call greenem
b = greenem_jph(mode,r1d,z1d,r2,z2);
if(~iscell(mode))
nbmodes=length({mode});
else
nbmodes=length(mode);
end
B=zeros(size(r,1),size(r,2),nbmodes);
% b is a rectangular matrix (length(r1d) x length(r2)) which needs to be multiplied by the
% column vector Ic to get the proper contribution of each source.
for i=1:size(b,3)
temp = b(:,:,i) * Ic;
B(:,:,i) = reshape(temp,size(r,1),size(r,2));
end
B(find((B==Inf)|(B==-Inf)))= 0;
B(find(isnan(B))) = 0;
return