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:
       
        
        
        % calculate the denominator of the biotsavart equation:
        
        
        
        % calculate the cross-produce section of the equation:
        
        
        
        % add the calculated field component to the running
        % total for this position:
        
        
        
        
    end
    
end