%**************************************************************************
%This function calculates field integrals (int B dz) from a set of magnetic charges.
%
%Input:
%field_int_pos, (nx2) array containing (x, y) positions where the field integrals are to be calculated
%charge, structure containing the magnetic charge information
%
%Output:
%field_int, structure containing the calculated field integrals
%   pos(nx2)            %Position [x(n), y(n)] of calculated field integrals in global coordinates (m)
%   int_bx_dz(nx1)      %Field integral (int Bx dz) at position [x(n), y(n)] (Tm)
%   int_by_dz(nx1)      %Field integral (int By dz) at position [x(n), y(n)] (Tm)
%
%Zachary Wolf, 4/12/10

function [field_int] = pm_field_calc_b_int(field_int_pos, charge)

%Global variables
global mu_0;

%Initialize the field integral structure
[r, c] = size(field_int_pos);
field_int = struct('pos', zeros(r, 2), 'int_bx_dz', zeros(r, 1), 'int_by_dz', zeros(r, 1));

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

%Simplify the notation
x_int = field_int_pos(:, 1)';
y_int = field_int_pos(:, 2)';

%Simplify the notation
q = [charge.q];
x_q = [charge.x];
y_q = [charge.y];

%Loop over the field points
for ii = 1 : r
    
    %Select a field point
    x_int_i = x_int(ii);
    y_int_i = y_int(ii);

    %Calculate the field integral from all the charges
    rho_mag_sq = (x_int_i - x_q).^2 + (y_int_i - y_q).^2;
    int_bx_dz = (mu_0 / (2 * pi)) * sum((x_int_i - x_q) .* q ./ rho_mag_sq);
    int_by_dz = (mu_0 / (2 * pi)) * sum((y_int_i - y_q) .* q ./ rho_mag_sq);

    %Save the field integral values
    field_int.int_bx_dz(ii) = int_bx_dz;
    field_int.int_by_dz(ii) = int_by_dz;

%End loop over the field points
end

%Done
end
