function [Bx,By,Bz] = biotsavart_helmholtzpair(r);
%----------------------------------------------------------------
% This function calculates the magnetic field due to a helmholtz
% pair.
%
% Inputs:
% r - array of points at which to calculate the magnetic field
% - three rows, Nr columns
% - each column contains coords for single point
% - first row contains x coordinates, second y, third z
% - units of meters
%
% Outputs:
% Bx - x components of magnetic field at locations specified
% by "r" input array
% By - y components of magnetic field
% Bz - z components of magnetic field
%
% Notes:
% - the specific parameters defining the helmholtz pair have
% been hard-coded into this function.
% - the radius of the coils is 15cm
% - there are 120 windings per section of the pair
% - there is assumed to be 3A of current flowing
% - the loops are discretized into 48 elements each
%
% Date: October 25, 2007
% Originator: Blaine Chronik
%----------------------------------------------------------------
% define helmholtz pair parameters:
%
Nloop = 120; % number of windings for each coil
Nseg = 40; % number of segments to discretize coil into
I = 3; % 3 A current flowing through entire structure
a = 0.15; % radius of coils in m
% create discrete approximation of single loop:
%
delta_theta = 2*pi./Nseg;
theta = (0:delta_theta:(2*pi - delta_theta));
% Note: we drop the last theta value, because 2pi is same
% as angle of zero.
% now convert array of angles into array of x,y positions for
% current elements in first loop:
%
xloop1 = a.*cos(theta);
yloop1 = a.*sin(theta);
zloop1 = (a/2).*ones(1,Nseg);
% Note: the z coords of the first loop are all the same:
% assumed to be half the radius in the positive z direction
% now, create x,y,z coords for elements in second loop (the
% negative-z loop
%
xloop2 = a.*cos(theta);
yloop2 = a.*sin(theta);
zloop2 = (-a/2).*ones(1,Nseg);
% Note: the z coords of the second loop are all the same:
% assumed to be half the radius in the negative z direction
% now we create the arrays containing the sizes and directions
% of each element:
%
dl = 2*pi*a./Nseg; % all segments are same length, given by
% circumference of circle, divided by number
% of segments in loop
dxloop1 = -dl.*sin(theta);
dyloop1 = dl.*cos(theta);
dzloop1 = zeros(1,Nseg);
dxloop2 = -dl.*sin(theta);
dyloop2 = dl.*cos(theta);
dzloop2 = zeros(1,Nseg);
% Now, we "assemble" the 6 vectors into a single array holding
% all of them together: the "elements array":
%
xcoil = [xloop1 xloop2];
ycoil = [yloop1 yloop2];
zcoil = [zloop1 zloop2];
dxcoil = [dxloop1 dxloop2];
dycoil = [dyloop1 dyloop2];
dzcoil = [dzloop1 dzloop2];
% Note: the above operation concatenates the arrays
% to create a single row array
E = [xcoil; ycoil; zcoil; dxcoil; dycoil; dzcoil];
% Note: the above operation combines the 6 1xN arrays
% into a single 6xN array
%====== now enter Biot-Savart calculation ========
% Now that we have an array "E" that contains all the information
% about the wire pattern, we can calculate the magnetic field
% at the necessary locations...
%
mu = 4*pi*1E-7; % define mu-not
const = mu*I*Nloop/(4*pi); % define constant coefficient
Bx = zeros(1,length(r));
By = zeros(1,length(r));
Bz = zeros(1,length(r));
% the idea here, is that we'll loop through the individual
% positions, and for each position, we loop through the entire
% elements array, adding up the contributions to the magnetic
% field at the individual location.
for i = 1:length(r), % loop over all positions
for j = 1:length(E), % loop over all elements of the wire
% calculate the components of the (r - r') vector:
v =[(0-xcoil) (0-ycoil) (z-zcoil)];
% calculate the denominator of the biotsavart equation:
d=r.^3;
% calculate the cross-produce section of the equation:
cross(dl,r)
% add the calculated field component to the running
% total for this position:
end
end