%**************************************************************************
%This function calculates the field (B) from a set of magnetic charges.
%
%Input:
%param, structure containing program parameters
%charge, structure containing the magnetic charge information
%
%Output:
%mag_field, structure containing the calculated magnetic field
%   x           %X position of the field point (m)
%   y           %Y position of the field point (m)
%   z           %Z position of the field point (m)
%   bx          %Bx magnetic field component at (x, y, z) (T)
%   by          %By magnetic field component at (x, y, z) (T)
%   bz          %Bz magnetic field component at (x, y, z) (T)
%
%Zachary Wolf, 3/31/08

function [mag_field] = ppm_calc_mag_field(param, charge)

%Global variables
global mu_0;

%Put all points where the field is to be calculated into vectors
index = 0;
for ii = 1 : length(param.x_field_pts)
    for jj = 1 : length(param.y_field_pts)
        for kk = 1 : length(param.z_field_pts)
            index = index + 1;
            mag_field.x(index) = param.x_field_pts(ii);
            mag_field.y(index) = param.y_field_pts(jj);
            mag_field.z(index) = param.z_field_pts(kk);
        end
    end
end

%Initialize the magnetic field
mag_field.bx = zeros(1, length(mag_field.x));
mag_field.by = zeros(1, length(mag_field.x));
mag_field.bz = zeros(1, length(mag_field.x));
    
%Simplify the field notation
x_b = mag_field.x;
y_b = mag_field.y;
z_b = mag_field.z;

%Loop over the charges
for ii = 1 : charge.num_charge

    %Extract the charge information
    q = charge.q(ii);
    x_q = charge.x(ii);
    y_q = charge.y(ii);
    z_q = charge.z(ii);

    %Calculate the field from the charge
    bx_add = (mu_0 / (4 * pi)) * (x_b - x_q) * q ./ sqrt((x_b - x_q).^2 + (y_b - y_q).^2 + (z_b - z_q).^2).^3;
    by_add = (mu_0 / (4 * pi)) * (y_b - y_q) * q ./ sqrt((x_b - x_q).^2 + (y_b - y_q).^2 + (z_b - z_q).^2).^3;
    bz_add = (mu_0 / (4 * pi)) * (z_b - z_q) * q ./ sqrt((x_b - x_q).^2 + (y_b - y_q).^2 + (z_b - z_q).^2).^3;
    
    %Avoid singularities
    %index_array = find(sqrt((x_b - x_q).^2 + (y_b - y_q).^2 + (z_b - z_q).^2) < 1.e-6);
    %bx_add(index_array) = 0;
    %by_add(index_array) = 0;
    %bz_add(index_array) = 0;

    %Add the field to the field from the other charges
    mag_field.bx = mag_field.bx + bx_add;
    mag_field.by = mag_field.by + by_add;
    mag_field.bz = mag_field.bz + bz_add;

%End loop over the charges
end

%Done