% GKLmagnet3.m: magnetic profile for TWO ring magnets plus a THIRD 
% auxiliary gun coil magnet.
%   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.
%   The effect of a third weaker coil at a prescribed distance from 
% the magnetic center is also included.  This magnet is useful for 
% reducing the length required for compression, but it can distort 
% the uniformity profile.  Use a negative magnetic field strength for 
% the subtraction mode, positive for add-mode.
%   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.
%   A pink star marks the location of the electron gun cathode, which 
% is located at the point where the magnetic field falls to comp*Bmax, 
% the magnetic compression point.
%
% NOTE: Requires 'find_intersection.m' to find gun position.  This file
%   is available from http://web.mit.edu/cjoye/www/research/research.html
%
% run GKLmagnet3
%
% Last edit by Colin Joye, 8/12/03
clear all; 
opengl autoselect;
% -------------- USER INPUT -------------------------
Bmax = 6.2;             %[Telsa] maximum magnetic field in uniform region
Tol  = 1.0;             %[%] percent homogeneity of uniform region = +/-(TOL/2)

% Main magnet sizes.  The magnets are finite-length finite-thickness hollow 
% cylinders modeled as a single turn rectangular cross-section with uniform 
% current density. The magnets are centered around z=0, extending 'Length'/2 
% above and below z=0.
RadIn  = 16;            %[cm] inner radius of main coils
RadOut = 16.2;          %[cm] outer radius of main coils
Length = 2.5;           %[cm] main coil lengths
Separ  = 16.9;          %[cm] center separation distance > Length
zaxis  = 35;            %[cm] length above z=0 to view the field profile

% -- Auxiliary magnet -- 
ABmax   = -0.75;         %[T] max aux field (use -x.x for subtract mode)
ARadIn  = 18;            %[cm] inner radius of auxiliary coil
ARadOut = 20;            %[cm] outer radius of auxiliary coil
ALength = 2.0;           %[cm] auxiliary coil length
ALoc    = -35.0;         %[cm] location of center of aux coil from z=0
comp    = 1/25;          %[] magnetic compression ratio

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

% check if length is larger than separation.
if(Separ<Length),
    disp(['Error, Separation must be greater than Length']);
    break;
end
% make sure aux magnet doesn't touch main magnet.
if( (2*ALoc+ALength+Length+Separ) >= 0  ),
    disp(['Error, ALoc must be larger to prevent magnets from touching']);
    break;
end
flag = 0;

% Define z-axis vector for field profile.
z = zaxis*[(ALoc/zaxis-1/2):1/2e4:0 0:1/2e4:1] ;  
%k = [0 1 2 5 10 15 20 25 30 35 40 45 53]; %like Kreischer's Table II.
%z = [-fliplr(k) k(2:end)];

% 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:
sep = Separ/2;
L   = Length;
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

% now run calculation again for auxiliary magnet
Ar1 = ARadIn;
Ar2 = ARadOut;
AL  = ALength;
r1p = Ar1 + sqrt(Ar1^2 + (AL/2 + (z - ALoc)).^2);
r1n = Ar1 + sqrt(Ar1^2 + (AL/2 - (z - ALoc)).^2);
r2p = Ar2 + sqrt(Ar2^2 + (AL/2 + (z - ALoc)).^2);
r2n = Ar2 + sqrt(Ar2^2 + (AL/2 - (z - ALoc)).^2);
C0  = (z - ALoc) .* log( r1n./r1p .* r2p./r2n );
C1  = AL/2 * log( r2n .* r2p ./ r1n ./ r1p );
A0  = 1 / (Ar2 - Ar1) / AL; % multiplicative value
Hz  = A0 * (C0 + C1);
Bz  = Bz + Hz * ABmax / max(Hz); % normalizing for desired value of Bmax.

% Renormalize again, since addition of fields messes up first normalization.
Bz  = Bz*Bmax/max(Bz);
LBz = length(Bz);

% Find uniformity zones.
u     = find( abs(Bz - Bmax)/Bmax <= (Tol/100) );
split = find( diff(u)>1 );
if( ~isempty(split) ), %check if the field is still continuous.
    disp(['Warning, field is broken into two regions or is distorted severely']);
    UniformLength = 0;
    flag = 1;
else
    UniformLength = z(u(end)) - z(u(1));
    split = 1;
end

% ---------------- Figures -------------------
% Draw cylinders
[Xi,Yi,Zi] = cylinder(RadIn,100);
[Xo,Yo,Zo] = cylinder(RadOut,100);
s     = ceil( cyl_frac*size(Xi,2) );
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);

% Draw aux magnet cylinder
[AXi,AYi,AZi] = cylinder(ARadIn,100);
[AXo,AYo,AZo] = cylinder(ARadOut,100);
s      = ceil( cyl_frac*size(AXi,2) );
AcylXi = (AZi(:,1:s)-0.5)*AL;
AcylYi = -AXi(:,1:s);
AcylZi = -AYi(:,1:s);
AcylXo = (AZo(:,1:s)-0.5)*AL;
AcylYo = -AXo(:,1:s);
AcylZo = -AYo(:,1:s);

figure(1)
clf(1)
hold on;
surf(cylXi-sep,cylYi,cylZi);     % inner cylinder 1 
surf(cylXo-sep,cylYo,cylZo);     % outer cylinder 1
surf(cylXi+sep,cylYi,cylZi);     % inner cylinder 2
surf(cylXo+sep,cylYo,cylZo);     % outer cylinder 2
surf(AcylXi+ALoc,AcylYi,AcylZi); % inner cylinder aux
surf(AcylXo+ALoc,AcylYo,AcylZo); % outer cylinder aux
colormap autumn;
shading flat;

% Plot dressups
a = 1.2*max(RadOut,ARadOut);
b = 1.1*zaxis;
c = b/2-ALoc;
axis equal;
axis([-c 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([-c -c 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([-c -c 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.2 1 0.5]);
%alpha(h,.3);
%set(h,'EdgeAlpha',0);

% solidify cylinders
% y-axis patch for coils 1 & 2
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);
% x-axis patch for coils 1 & 2
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);
% y-axis patch for aux
px = AcylXi(:,end)'*[1 0 0 1;0 1 1 0];
py = [AcylYi(:,end)' AcylYo(:,end)'];
pz = [AcylZi(:,end)' AcylZo(:,end)'];
h  = patch(px+ALoc, py, pz,[0.7 0.8 0.9]);
set(h,'EdgeAlpha',transparency);
% x-axis patch for aux
px = AcylXi(:,1)'*[1 0 0 1;0 1 1 0];
py = [AcylYi(:,1)' AcylYo(:,1)'];
pz = [AcylZi(:,1)' AcylZo(:,1)'];
h  = patch(px+ALoc, py, pz,[0.7 0.8 0.9]);
set(h,'EdgeAlpha',transparency);

% Put grid on field line patch
ytick = get(gca,'Ztick'); 
ytick = (ytick + ytick(end))/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);
y2 = c*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([-y2;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(-[c c],[0 0],[0 a],'Color',[0 0.3 1]); % vert axis line for numbers
line(-[1;1.03]*y2,[y0;y0],[1;1]*ytick,'Color',[0 0.3 1]); % draw tick lines.
h = text(-1.03*y2',y0',ytick',num2str(ytick'/scale),'HorizontalAlignment','right');
set(h,'Color',[0 0.3 1]);
pz = length(get(h,'String'))-3;
h = text(-1.05*c-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/2+ALoc 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,LBz),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.
u0 = u(1);
u1 = u(end);
us0= u(split);
us1= u(split+1);
px = [z(u0) z(us0) z(us0) z(u0)];
py = [0 0 0 0];
pz = scale*[0 0 Bz(u1) Bz(u0)];
h = patch(px, py, pz,[1 0.2 0.5]);
alpha(h,1.5*transparency);
set(h,'EdgeAlpha',0);
px = [z(u1) z(us1) z(us1) z(u1)];
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])

% Add point showing 'comp' compression point
[zcomp,ycomp] = find_intersection(z,Bz,Bmax*comp*ones(1,LBz));
gp = [];
if ~isempty(zcomp),
  plot3(zcomp(1),0,0,'m*','MarkerSize',14);
  gp = ['; Gun Position = ' num2str(zcomp,3) ' cm'];
end

% List parameters on plot
param12 = ['R1 = ' num2str(r1) ' cm; R2 = ' num2str(r2) ' cm; L = ' num2str(L) ...
           ' cm; S = ' num2str(Separ) ' cm; '  ];
paramA  = ['BAux = ' num2str(ABmax) ' T; AR1 = ' num2str(Ar1) ' cm; AR2 = ' ...
           num2str(Ar2) ' cm; AL = ' num2str(AL) ' cm; ALoc = ' num2str(ALoc) ' cm; '];
text(-b/2,-a/2,-RadOut,param12);
text(-b/2,-a/2,-RadOut-4,paramA);

% Title and Labels
title(['Magnet strucure field profile w/ auxiliary magnet : ' num2str(Bmax) ... 
       ' T max; ' num2str(Tol) '% Uniform length: ' num2str(UniformLength,3) ' cm' gp]);
xlabel('Z-axis [cm]');
ylabel('X-axis [cm]');
zlabel('Y-axis [cm]');

hold off;
% ------------- END GKLmagnet2.m ---------------------------