% OptimalMagnet.m
%   This program takes a range of possible values for the radii of a finite
% length, nonzero thickness cylindrical shell superconducting magnet and 
% optimizes the radii and length to meet a user defined specification on 
% percent homogeneity over a specified region of uniform length.  The magnetic
% field profile is calculated on-axis.
%   The user specifies the minimum inner radius, maximum outer radius, 
% minimum thickness and maximum length.  The output of the program is the  
% set of radii and lengths which meet the given specs, if any can be found.
%   From these optimal data sets, several can be chosen and run in GKLmagnet.m 
% (available from http://web.mit.edu/cjoye/www/research/matlab/GKLmagnet.m)
% to visualize the field profile and check the result.
%  Note: if too many optimal sets result, increase 'UL' or decrease 'Tol'.
%
% run research/optimalmagnet
%
% Last edit by Colin Joye, 7/14/03
clear all;
% -------------- USER INPUT -------------------------
Bmax = 5.12;             %[Telsa] maximum magnetic field in uniform region
Tol  = 1.0;             %[%] percent homogeneity of uniform region
UL   = 5.8;               %[cm] length to be uniform to within 'Tol'

% 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.
R1min  = 18;            %[cm] minimum inner radius of coil
R2max  = 45;            %[cm] maximum outer radius of coil > R1min
Tmin   = 0.5;             %[cm] minimum thickness of coil
Lmax   = 45;            %[cm] maximum coil length > R1min
zaxis  = 40;            %[cm] length above z=0 to view the field profile

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

step = 10; % sets the amount of data to compute through.
% 4 degrees of freedom: r1, r2, L, and z
r1 = [R1min:R1min/step:(R2max-Tmin)];
r2 = [(R1min+Tmin):R2max/step:R2max];
L  = R1min:Lmax/step:Lmax;
z  = 0:zaxis/5e3:zaxis;

% define lengths of vectors
Lr1 = length(r1);
Lr2 = length(r2);
LL  = length(L);
Lz  = length(z);

% define Hz, Bz and uniformity vectors/matrices
Hz = zeros(Lr1,Lr2,LL,Lz);
Bz = Hz;
u  = zeros(Lr1,Lr2,LL); % uniformity matrix

% First loop computes Hz field on axis for each data set.
h = waitbar(0,'Computing Field Matrix');
for n1 = 1:Lr1,
  for n2 = 1:Lr2,
    if(r1(n1) <= (r2(n2)-Tmin)),
      for nL = 1:LL,
        for nz = 1:Lz,
            r1p = r1(n1) + sqrt(r1(n1)^2 + (L(nL)/2 + z(nz))^2);
            r1n = r1(n1) + sqrt(r1(n1)^2 + (L(nL)/2 - z(nz))^2);
            r2p = r2(n2) + sqrt(r2(n2)^2 + (L(nL)/2 + z(nz))^2);
            r2n = r2(n2) + sqrt(r2(n2)^2 + (L(nL)/2 - z(nz))^2);

            C0 = z(nz) .* log( r1n/r1p * r2p/r2n );
            C1 = L(nL)/2 * log( r2n * r2p / r1n / r1p );
            A0 = 1 / (r2(n2)-r1(n1)) / L(nL); % multiplicative value
            Hz(n1,n2,nL,nz) = A0 * ( C0 + C1 );    
        end
      end
    end
  end
waitbar(2*n1/3/Lr1,h);
end

% second loop computes scaled Bz field and finds optimum parameters.
for n1 = 1:Lr1,
  for n2 = 1:Lr2,
    for nL = 1:LL,
      if( any(Hz(n1,n2,nL,:)) )
        Bz(n1,n2,nL,:) = Hz(n1,n2,nL,:) / max(max(max( Hz(n1,n2,nL,1) ))) * Bmax; % normalizing for desired value of Bmax.
        bn = abs(Bz(n1,n2,nL,:) - Bmax)/Bmax;
        s = find( bn < (Tol/100) );
        if( ~isempty(s) ),
          u(n1,n2,nL) = 2*z(s(end));
        end
      end    
    end
  end
waitbar(n1/Lr1,h);
end
close(h);

% find configurations where maximum uniformity occurs to within 'Tol' percent.
%[n1 n2 nL] = ind2sub(size(Hz),find( abs(Bmax - Bz(:,:,:,1))/Bmax < Tol/100 ) );
uopt = u( ind2sub( size(u), find( u >= UL ) ) ); % find sets meeting spec.

% find where optimal values occur in input ranges.
Luopt = length(uopt);
n1 = zeros(1,1:Luopt);
n2 = n1;
nL = n1;
for k = 1:Luopt,
  [a b c] = ind2sub( size(u), find( u == uopt(k)) );
  for m=1:length(a),
    n1(k+m-1) = a(m);
    n2(k+m-1) = b(m);
    nL(k+m-1) = c(m);
  end
  k = k+length(a)-1;
end

% show results.
umax = max(max(max(u)))
if( all([n1 n2 nL]~=0) & ~isempty([n1 n2 nL]) ),
  opt_umax = uopt'
  opt_r1 = r1(n1)
  opt_r2 = r2(n2)
  opt_L  = L(nL)
else
  disp(['** No parameters found to meet requirements **']);
end

% ---------------- END OptimalMagnet.m ------------------