# Transforming a matlab code!

1. Jan 24, 2010

Hi fellows!!
I need a bit of help with transforming a matlab code that simulate magnetic field of two concentric solenoid TO a plot that simulate the magnetic fields of high voltage transmission lines (specificly two and four lines), how can it be modified? here is the code:

% 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.
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;
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:
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;

% Plot dressups
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 ---------------------------

PS: It dooes what i need but i need the magnetic fields of those high voltage transsmision lines, and i dont know to much of programing, sooo i will apreciate a loooot if you can help me out =)

2. Jan 30, 2010

3. Jan 30, 2010

### uart

Ok then. Hang around the Computer Lab while other students from your class are working on this project. After they leave, look though all the waste paper bins in the hope they've disgarded some hard copies of partially debugged code. Then you'll have some code that's a little bit closer to what you're trying to achieve to try and hack.

Last edited: Jan 31, 2010
4. Jan 31, 2010