function B = currentLoop(x,y,z,dTh)
%B = currentLoop(x,y,z,dTh);
% calculate B-field at (x,y,z) from a unit-radius current loop centered in
% the x-y plane, w/ current 4pi/mu_o

% pre-allocate B
B = zeros(3,1); % note vector quantity

% loop over all theta
for th = 0:dTh:2*pi-dTh
    
    % r-vector parameters
    r = [
        x - cos(th)
        y - sin(th)
        z
        ];
    r3 = sqrt(sum(r.^2)).^3;
    
    % delta s
    ds = dTh * [
        -sin(th)
        cos(th)
        0 ];
    
    % cross product
    cp = cross(ds,r);
    
    % sum up B results
    B = B + cp/r3;
    
end

% theoretical result for x = y = 0
% Be = 2*pi/(1+z^2)^(3/2);

return