% calculate B-field from a single current loop of radius 1, centered in
% the x-y plane, with I = 4pi/mu_o

% integration precision
dTh = pi/10;

%% first, let's compare to theoretical results along z-axis
z = -1:0.1:1;
B = zeros(3,length(z));
ze = -1:0.01:1;
Bz = 2*pi./(1+ze.^2).^(3/2); % exact result

for k = 1:length(z)
    B(:,k) = currentLoop(0,0,z(k),dTh);
end

figure(1)
clf
plot(ze,Bz,'b-',z,B(3,:),'rd')
xlabel('z')
ylabel('B_z')
title('B-field from current loop along axis')


%% next, image x-z plane through center of loop

% establish mesh in space & pre-allocate results
x = -6:1:6;
% x = -6:2:6;
% x = -2:1/3:2;
z = x;
z = min(x):diff(x(1:2)):max(x);
[X,Z] = meshgrid(x,z);
Bx = zeros(size(X));
Bz = Bx;

% loop over all space to find B-field
for ii = 1:length(z) % note orientation of x & z here!  z == rows
    for kk = 1:length(x)
        % calculate B at point
        B = currentLoop(X(ii,kk),0,Z(ii,kk),dTh);
        % assign to output variables
        Bx(ii,kk) = B(1);
        Bz(ii,kk) = B(3);
    end
end

% [ii,kk] = find( abs(X)-1 <= eps*1e6 & abs(Z) <= eps*1e6 )
% Bx(ii,kk) = 0;
% Bz(ii,kk) = 0;


figure(2)
clf
quiver(x,z,Bx,Bz)
ylabel('z')
xlabel('x')
powerpoint

return