%**************************************************************************
%This function calculates the magnetic charges which are used to calculate
%the field from a permanent magnet.
%
%Input:
%param, parameter structure
%block, structure containing the permanent magnet block information
%
%Output:
%charge, structure containing the charge information
%   num_charge	%Number of charges
%   q           %Magnetic charge, M.n dA (Am)
%   x           %X position of charge (m)
%   y           %Y position of charge (m)
%   z           %Z position of charge (m)
%   block_num   %Block number the charge came from
%
%Zachary Wolf, 3/31/08

function [charge] = ppm_calc_charge(param, block)

%Block surface division parameters
num_xdiv = param.num_block_xdiv;
num_ydiv = param.num_block_ydiv;
num_zdiv = param.num_block_zdiv;

%Initialize
charge_num = 0;

%Loop over blocks
for kk = 1 : block.num_block

    %Extract the block information
    mx = block.mx(kk);
    my = block.my(kk);
    mz = block.mz(kk);
    x_c = block.x_center(kk);
    y_c = block.y_center(kk);
    z_c = block.z_center(kk);
    x_w = block.x_width(kk);
    y_w = block.y_width(kk);
    z_w = block.z_width(kk);

    %Initialize
    dx = x_w / num_xdiv;
    dy = y_w / num_ydiv;
    dz = z_w / num_zdiv;
    da_x = dy * dz;
    da_y = dz * dx;
    da_z = dx * dy;

    %X+ face
    for ii = 1 : num_ydiv
        for jj = 1 : num_zdiv
            charge_num = charge_num + 1;
            charge.q(charge_num) = mx * da_x;
            charge.x(charge_num) = x_c + (x_w / 2);
            charge.y(charge_num) = y_c - (y_w / 2) + (dy / 2) + dy * (ii - 1);
            charge.z(charge_num) = z_c - (z_w / 2) + (dz / 2) + dz * (jj - 1);
            charge.block_num(charge_num) = kk;
        end
    end

    %X- face
    for ii = 1 : num_ydiv
        for jj = 1 : num_zdiv
            charge_num = charge_num + 1;
            charge.q(charge_num) = -mx * da_x;
            charge.x(charge_num) = x_c - (x_w / 2);
            charge.y(charge_num) = y_c - (y_w / 2) + (dy / 2) + dy * (ii - 1);
            charge.z(charge_num) = z_c - (z_w / 2) + (dz / 2) + dz * (jj - 1);
            charge.block_num(charge_num) = kk;
        end
    end

    %Y+ face
    for ii = 1 : num_zdiv
        for jj = 1 : num_xdiv
            charge_num = charge_num + 1;
            charge.q(charge_num) = my * da_y;
            charge.x(charge_num) = x_c - (x_w / 2) + (dx / 2) + dx * (jj - 1);
            charge.y(charge_num) = y_c + (y_w / 2);
            charge.z(charge_num) = z_c - (z_w / 2) + (dz / 2) + dz * (ii - 1);
            charge.block_num(charge_num) = kk;
        end
    end

    %Y- face
    for ii = 1 : num_zdiv
        for jj = 1 : num_xdiv
            charge_num = charge_num + 1;
            charge.q(charge_num) = -my * da_y;
            charge.x(charge_num) = x_c - (x_w / 2) + (dx / 2) + dx * (jj - 1);
            charge.y(charge_num) = y_c - (y_w / 2);
            charge.z(charge_num) = z_c - (z_w / 2) + (dz / 2) + dz * (ii - 1);
            charge.block_num(charge_num) = kk;
        end
    end

    %Z+ face
    for ii = 1 : num_xdiv
        for jj = 1 : num_ydiv
            charge_num = charge_num + 1;
            charge.q(charge_num) = mz * da_z;
            charge.x(charge_num) = x_c - (x_w / 2) + (dx / 2) + dx * (ii - 1);
            charge.y(charge_num) = y_c - (y_w / 2) + (dy / 2) + dy * (jj - 1);
            charge.z(charge_num) = z_c + (z_w / 2);
            charge.block_num(charge_num) = kk;
        end
    end

    %Z- face
    for ii = 1 : num_xdiv
        for jj = 1 : num_ydiv
            charge_num = charge_num + 1;
            charge.q(charge_num) = -mz * da_z;
            charge.x(charge_num) = x_c - (x_w / 2) + (dx / 2) + dx * (ii - 1);
            charge.y(charge_num) = y_c - (y_w / 2) + (dy / 2) + dy * (jj - 1);
            charge.z(charge_num) = z_c - (z_w / 2);
            charge.block_num(charge_num) = kk;
        end
    end

%End loop over blocks
end

%Add the number of charges to the structure
charge.num_charge = charge_num;

%Done