% GKLmagnet.m
%   This code simulates the magnetic field on the z-axis 
% a concentric solenoidal superconducting magnet of specified
% size for use in the 140GHz Gyroklystron (GKL) experiment, MIT.
%
% run GKLmagnet
%
% Last edit by Colin Joye, 7/10/03
clear all;
opengl autoselect;
% -------------- USER INPUT -------------------------
Bmax = 5.12;            %[Telsa] maximum magnetic field in uniform region
Tol  = 0.6;            %[%] percent homogeneity of uniform region

% Magnet size.  The magnet is a finite-length finite-thickness hollow cylinder
% modeled as a single turn rectangular cross-section current density.
% The magnet is centered at z=0, extending 'Length'/2 above and below z=0.
RadIn  = 19.8;           %[cm] inner radius of coil
RadOut = 32;           %[cm] outer radius of coil
Length = 36;           %[cm] coil length
zaxis  = 40;           %[cm] length above z=0 to view the field profile

% -------------- End USER INPUT ---------------------

transparency = 0.28; % controls transparency of slices & patches. '= 1.0' is opaque.
cyl_frac = 0.74;    % azimuthal fraction of cylinder to plot. '= 1.0' for full cyl.

% Define z-axis vector for field profile.
z = [0:zaxis/2e4:zaxis];  
z = [0 1 2 5 10 15 20 25 30 31]; %Kreischer's Table II.

% Draw cylinders
L = Length;
[Xi,Yi,Zi] = cylinder(RadIn,100);
[Xo,Yo,Zo] = cylinder(RadOut,100);
s = floor( cyl_frac*size(Xi,2) )+1;
cylXi = (Zi(:,1:s)-0.5)*L;
cylYi = -Xi(:,1:s);
cylZi = -Yi(:,1:s);
cylXo = (Zo(:,1:s)-0.5)*L;
cylYo = -Xo(:,1:s);
cylZo = -Yo(:,1:s);

% Applying Biot-Savart law to a finite-length finite-thickness cylinder 
% with a rectangular cross-section of uniform current density, we get
r1 = RadIn;
r2 = RadOut;
r1p = r1 + sqrt(r1^2 + (L/2 + z).^2);
r1n = r1 + sqrt(r1^2 + (L/2 - z).^2);
r2p = r2 + sqrt(r2^2 + (L/2 + z).^2);
r2n = r2 + sqrt(r2^2 + (L/2 - z).^2);

C0 = z .* log( r1n./r1p .* r2p./r2n );
C1 = L/2 * log( r2n .* r2p ./ r1n ./ r1p );
A0 = 1 / (r2-r1) / L; % multiplicative value

Hz = A0 * ( C0 + C1 );
Bz = Hz * Bmax / max(Hz); % normalizing for desired value of Bmax.

zp = [-fliplr(z) z];
Bzp = [fliplr(Bz) Bz];

% Find uniformity zones.
u = find( abs(Bzp - Bmax)/Bmax <= (Tol/100) );
UniformLength = zp(u(end))-zp(u(1));


% ---------------- Figures -------------------
figure(1)
clf(1)
hold on;
surf(cylXi,cylYi,cylZi); % inner cylinder
surf(cylXo,cylYo,cylZo); % outer cylinder
colormap autumn;
shading flat;

% Plot dressups
a = 1.2*RadOut;
b = 1.1*zaxis;
axis equal;
axis([-b b -a a -a a]);
view([-15 15]);

% define scaling variable for making field plot look better.
scale = floor(a/max(Bzp)); % scale line and grid to view field better
if(scale<=1) % make sure you dont have scale = 0, 3, 6, 7, 9, etc.
    scale=1; 
elseif(mod(scale,2)~=0 & mod(scale,5)~=0)
    scale=scale-1;
end

% Add transparent slices.  Transparency controlled by "alpha"
% Horizontal plane slice
h = patch([-b -b b b],[-a a a -a],[0 0 0 0],[0.2 1 0.5]);
alpha(h,transparency);
set(h,'EdgeAlpha',0);
% Vertical plane slice
h = patch([-b -b b b],[0 0 0 0],[-a a a -a],[0.2 1 0.5]);
alpha(h,transparency);
set(h,'EdgeAlpha',0);
% Z=0 plane slice
%h=patch([0 0 0 0],[-a -a a a],[-a a a -a],[0.7 1 0]);
%alpha(h,.3);
%set(h,'EdgeAlpha',0);

% solidify cylinders
map = colormap;
px = cylXi(:,end)'*[1 0 0 1;0 1 1 0];
py = [cylYi(:,end)' cylYo(:,end)'];
pz = [cylZi(:,end)' cylZo(:,end)'];
h  = patch(px, py, pz,[0.7 0.8 0.9]);
set(h,'EdgeAlpha',transparency);
px = cylXi(:,1)'*[1 0 0 1;0 1 1 0];
py = [cylYi(:,1)' cylYo(:,1)'];
pz = [cylZi(:,1)' cylZo(:,1)'];
h  = patch(px, py, pz,[0.7 0.8 0.9]);
set(h,'EdgeAlpha',transparency);

% Put grid on field line patch
ytick = get(gca,'Ztick'); 
ytick = (ytick + max(ytick))/2; % remove negative ticks and interpolate.
ztick = get(gca,'Xtick');
Ly = length(ytick);
Lz = length(ztick);
y0 = zeros(1,Ly);
y1 = b*ones(1,Ly);
z0 = zeros(1,Lz);
z1 = a*ones(1,Lz);
line([1;1]*ztick,[z0;z0],[z0;z1],'Color',0.4*[1 1 1],'LineStyle',:); %vert lines
line([-y1;y1],[y0;y0],[1;1]*ytick,'Color',0.4*[1 1 1],'LineStyle',:);%horiz lines
% Add axis line and numbers on new grid.
% Vertical:
line(-[b b],[0 0],[0 a],'Color',[0 0.3 1]); % vert axis line for numbers
line(-[1;1.03]*y1,[y0;y0],[1;1]*ytick,'Color',[0 0.3 1]); % draw tick lines.
h = text(-1.03*y1',y0',ytick',num2str(ytick'/scale),'HorizontalAlignment','right');
set(h,'Color',[0 0.3 1]);
pz = length(get(h,'String'))-3;
h = text(-1.1*b-pz,0,mean(ytick),'Field [T] (blue)','HorizontalAlignment','center');
set(h,'Color',[0 0.3 1],'Rotation',90);
% Horizontal:
line(-[b zaxis],[0 0],[0 0],'Color',[0 0.3 1]); % horiz axis line for numbers
line( [zaxis b],[0 0],[0 0],'Color',[0 0.3 1]); % horiz axis line for numbers
line([1;1]*ztick,[z0;z0],-[0;0.03]*z1,'Color',[0 0.3 1]); % draw tick lines.
h = text(ztick',z0',-0.055*z1',num2str(ztick'),'HorizontalAlignment','center');
set(h,'Color',[0 0.3 1]);
h = line([-zaxis zaxis],[0,0],[0,0]); % z-axis line.
set(h,'Color','red','LineWidth',2.5);

% Add plot of field line
h = plot3(zp,zeros(1,length(Bzp)),scale*Bzp,'b');
set(h,'Color','blue','Linewidth',2);
h = text(0,0,0.1*Bmax,'Uniform Region');
set(h,'Rotation',90,'Color',[0 0.4 0.4])

% Add patch showing uniform region within 'Tol' percent.
px = [zp(u(1)) zp(u(end)) zp(u(end)) 0 zp(u(1))];
py = [0 0 0 0 0];
pz = scale*[0 0 Bzp(u(1)) Bmax Bzp(u(end))];
h = patch(px, py, pz,[1 0.2 0.5]);
alpha(h,1.5*transparency);
set(h,'EdgeAlpha',0);
%h = plot3(zp(u),zeros(1,length(u)),max(Bzp)*ones(1,length(u)),'m');
%set(h,'LineWidth',2');
h = text(0,0,0.1*Bmax,'Uniform Region');
set(h,'Rotation',90,'Color',[0 0.4 0.3])

% Title and Labels
title(['Magnet strucure w/ Field: ' num2str(Bmax) ' T max; R1 = ' num2str(r1) ...
    ' cm; R2 = ' num2str(r2) ' cm; L = ' num2str(L) ' cm; ' num2str(Tol) ...
    '% Uniform length: ' num2str(UniformLength) ' cm']);
xlabel('Z-axis [cm]');
ylabel('X-axis [cm]');
zlabel('Y-axis [cm]');

hold off;

% ------------- END GKLmagnet.m ---------------------------