%**************************************************************************
%This function calculates the field (B) from a set of magnetic charges.
%
%Input:
%field_pos, (nx3) array containing positions where the field is to be calculated
%charge, structure containing the magnetic charge information
%
%Output:
%field, structure containing the calculated magnetic field
%   pos(nx3)            %Position [x(n), y(n), z(n)] of calculated field in global coordinates (m)
%   b(nx3)              %Magnetic field [Bx(n), By(n), Bz(n)] at position [x(n), y(n), z(n)] (T)
%
%Zachary Wolf, 4/12/10

function [field] = pm_field_calc_b(field_pos, charge)

%Global variables
global mu_0;

%Initialize the magnetic field structure
[r, c] = size(field_pos);
field = struct('pos', zeros(r, 3), 'b', zeros(r, 3));

%Store the positions where the field is calculated
field.pos = field_pos;

%Simplify the notation
x_b = field_pos(:, 1)';
y_b = field_pos(:, 2)';
z_b = field_pos(:, 3)';

%Simplify the notation
q = [charge.q];
x_q = [charge.x];
y_q = [charge.y];
z_q = [charge.z];
% q = zeros(1, length(charge));
% x_q = zeros(1, length(charge));
% y_q = zeros(1, length(charge));
% z_q = zeros(1, length(charge));
% for ii = 1 : length(charge)
%     q(ii) = charge(ii).q;
%     x_q(ii) = charge(ii).pos(1);
%     y_q(ii) = charge(ii).pos(2);
%     z_q(ii) = charge(ii).pos(3);
% end
 
% %Loop over the charges
% for ii = 1 : length(charge)
%     
%     %Extract a charge
%     q_i = q(ii);
%     x_q_i = x_q(ii);
%     y_q_i = y_q(ii);
%     z_q_i = z_q(ii);
% 
%     %Calculate the field from the charge at all field points
%     bx_add = (mu_0 / (4 * pi)) * (x_b - x_q_i) * q_i ./ sqrt((x_b - x_q_i).^2 + (y_b - y_q_i).^2 + (z_b - z_q_i).^2).^3;
%     by_add = (mu_0 / (4 * pi)) * (y_b - y_q_i) * q_i ./ sqrt((x_b - x_q_i).^2 + (y_b - y_q_i).^2 + (z_b - z_q_i).^2).^3;
%     bz_add = (mu_0 / (4 * pi)) * (z_b - z_q_i) * q_i ./ sqrt((x_b - x_q_i).^2 + (y_b - y_q_i).^2 + (z_b - z_q_i).^2).^3;
% 
%     %Add the field to the field from the other charges
%     field.b = field.b + [bx_add', by_add', bz_add'];
% 
% %End loop over the charges
% end

%Loop over the field points
%tic;
for ii = 1 : r
    
    %Select a field point
    %if ii == 1, tic; end;
    x_b_i = x_b(ii);
    y_b_i = y_b(ii);
    z_b_i = z_b(ii);
    %if ii == 1, dt = toc; disp(['dt get x_b = ', num2str(dt)]); end;

    %Calculate the field from all the charges
    %if ii == 1, tic; end;
    r_mag = sqrt((x_b_i - x_q).^2 + (y_b_i - y_q).^2 + (z_b_i - z_q).^2);
    r_mag_cubed = r_mag.^3;
    bx = (mu_0 / (4 * pi)) * sum((x_b_i - x_q) .* q ./ r_mag_cubed);
    by = (mu_0 / (4 * pi)) * sum((y_b_i - y_q) .* q ./ r_mag_cubed);
    bz = (mu_0 / (4 * pi)) * sum((z_b_i - z_q) .* q ./ r_mag_cubed);
%     bx = sum((mu_0 / (4 * pi)) * (x_b_i - x_q) .* q ./ sqrt((x_b_i - x_q).^2 + (y_b_i - y_q).^2 + (z_b_i - z_q).^2).^3);
%     by = sum((mu_0 / (4 * pi)) * (y_b_i - y_q) .* q ./ sqrt((x_b_i - x_q).^2 + (y_b_i - y_q).^2 + (z_b_i - z_q).^2).^3);
%     bz = sum((mu_0 / (4 * pi)) * (z_b_i - z_q) .* q ./ sqrt((x_b_i - x_q).^2 + (y_b_i - y_q).^2 + (z_b_i - z_q).^2).^3);
    %if ii == 1, dt = toc; disp(['dt loop = ', num2str(dt)]); end;

    %Save the field values
    field.b(ii, :) = [bx, by, bz];

%End loop over the field points
end
%dt = toc;
%disp(['dt = ', num2str(dt), ' (loop total)']);

%Done
