% GKLmagnet2.m: magnetic profile for TWO ring magnets.
%   This code simulates the magnetic field on the z-axis of
% two concentric solenoidal superconducting magnets of specified
% size for use in the 140GHz Gyroklystron (GKL) experiment, MIT.
%   With two coils, it's possible to have one large region that 
% meets the tolerance spec, otherwise two fields are created which 
% may not meet spec.  If the field is not continuous (one field) 
% to within spec, the varialbe 'flag' is set to unity (otherwise 0),
% and a warning is displayed.
%
% run GKLmagnet2
%
% 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.3;             %[%] 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  = 15;            %[cm] inner radius of coil
RadOut = 18;            %[cm] outer radius of coil
Length = 6;             %[cm] coil length
Separ  = 18.0;          %[cm] center separation distance > Length
zaxis  = 30;            %[cm] length above z=0 to view the field profile

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

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

% check if length is larger than separation.
if(Separ<Length),
    disp(['Error, Separation must be greater than Length']);
    break;
end
flag = 0;

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

% Draw cylinders
sep = Separ/2;
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, 
% using an analytic (exact) result from Maple v8.0:
r1 = RadIn;
r2 = RadOut;
Bz = zeros(1,length(z));
for M=1:2, % loop calculates for both magnets.
r1p = r1 + sqrt(r1^2 + (L/2 + (z + sep)).^2);
r1n = r1 + sqrt(r1^2 + (L/2 - (z + sep)).^2);
r2p = r2 + sqrt(r2^2 + (L/2 + (z + sep)).^2);
r2n = r2 + sqrt(r2^2 + (L/2 - (z + sep)).^2);
C0 = (z + sep) .* 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 = Bz + Hz * Bmax / max(Hz); % normalizing for desired value of Bmax.
sep = -sep; % reverse sign for other magnet, reverse again on exit.
end
% Renormalize again, since addition of field messes up first normalization.
Bz = Bz*Bmax/max(Bz);

% Find uniformity zones.
zmid = find(z==0);
u = find( abs(Bz(zmid:end) - Bmax)/Bmax <= (Tol/100) );
UniformLength = 2*(z(u(end)+zmid)-z(u(1)+zmid));
if(u(1)~=1), %check if the field is still continuous.
    disp(['Warning, field is broken into two regions']);
    flag = 1;
end

% ---------------- Figures -------------------
figure(1)
clf(1)
hold on;
surf(cylXi-sep,cylYi,cylZi); % inner cylinder
surf(cylXo-sep,cylYo,cylZo); % outer cylinder
surf(cylXi+sep,cylYi,cylZi); % inner cylinder
surf(cylXo+sep,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(Bz)); % 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-sep, py, pz,[0.7 0.8 0.9]);
set(h,'EdgeAlpha',transparency);
h  = patch(px+sep, 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-sep, py, pz,[0.7 0.8 0.9]);
set(h,'EdgeAlpha',transparency);
h  = patch(px+sep, 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(z,zeros(1,length(Bz)),scale*Bz,'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.
uz0 = u(1)+zmid;
uz1 = u(end)+zmid;
px = [z(uz0) z(uz1) z(uz1) z(uz0)];
py = [0 0 0 0];
pz = scale*[0 0 Bz(u(1)+zmid) Bz(u(end)+zmid)];
h = patch(px, py, pz,[1 0.2 0.5]);
alpha(h,1.5*transparency);
set(h,'EdgeAlpha',0);
px = -[z(uz0) z(uz1) z(uz1) z(uz0)];
h = patch(px, py, pz,[1 0.2 0.5]);
alpha(h,1.5*transparency);
set(h,'EdgeAlpha',0);
%h = plot3(z(u),zeros(1,length(u)),max(Bz)*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; S = ' num2str(Separ) ...
    ' cm; ' num2str(Tol) '% Uniform length: ' num2str(UniformLength) ' cm']);
xlabel('Z-axis [cm]');
ylabel('X-axis [cm]');
zlabel('Y-axis [cm]');

hold off;

% ------------- END GKLmagnet2.m ---------------------------