%**************************************************************************
%This function calculates the magnetic charges which are used to calculate
%the field from a permanent magnet.  It operates on a single rectangular
%permanent magnet block.
%
%Input:
%block, structure containing the rectangular permanent magnet block information
%
%Output:
%charge(i), structure containing the charge information
%   q                   %Magnetic charge, M.n dA (Am)
%   pos_block(3)        %Position of charge [x, y, z] in block coordinates (m)
%   pos(3)              %Position of charge [x, y, z] in global coordinates (m)
%   x                   %x-position of charge in global coordinates, pos(1), for vectorization (m)
%   y                   %y-position of charge in global coordinates, pos(2), for vectorization (m)
%   z                   %z-position of charge in global coordinates, pos(3), for vectorization (m)
%   block_num           %Block number of the parent block
%
%Zachary Wolf, 4/8/10

function [charge] = pm_block_rect_calc_charge(block)

%Block surface division parameters
num_div = round(block.dim / block.div_size);  %number of divisions to give approximate division size
index_array = find(num_div < 1);
num_div(index_array) = 1;
div_size = block.dim ./ num_div;  %division size is approximately the specified size

%Extract the block information
block_num = block.num;
mx = block.m(1);
my = block.m(2);
mz = block.m(3);
x_w = block.dim(1);
y_w = block.dim(2);
z_w = block.dim(3);
num_xdiv = num_div(1);
num_ydiv = num_div(2);
num_zdiv = num_div(3);
dx = div_size(1);
dy = div_size(2);
dz = div_size(3);
da_x = dy * dz;
da_y = dz * dx;
da_z = dx * dy;
rot = block.rot;
trans = block.trans;

%Initialize
ch_num = 0;
charge(1 : 2 * num_ydiv * num_zdiv + 2 * num_xdiv * num_zdiv + 2 * num_xdiv * num_ydiv) = ...
    struct('q', 0, 'pos_block', [0, 0, 0], 'pos', [0, 0, 0], 'x', 0, 'y', 0, 'z', 0, 'block_num', 0);

%X+ face
for ii = 1 : num_ydiv
    for jj = 1 : num_zdiv
        q = mx * da_x;  %Leave inside loop for possible extension to non-uniform m
        if abs(q) > 0
            ch.q = q;
            ch.pos_block = [(x_w / 2), - (y_w / 2) + (dy / 2) + dy * (ii - 1), - (z_w / 2) + (dz / 2) + dz * (jj - 1)];  
            [ch.pos] = pm_charge_calc_rot_trans_pos(ch.pos_block, rot, trans);
            ch.x = ch.pos(1);
            ch.y = ch.pos(2);
            ch.z = ch.pos(3);
            ch.block_num = block_num;
            ch_num = ch_num + 1;
            charge(ch_num) = ch;
        end
    end
end

%X- face
for ii = 1 : num_ydiv
    for jj = 1 : num_zdiv
        q = -mx * da_x;
        if abs(q) > 0
            ch.q = q;
            ch.pos_block = [- (x_w / 2),  - (y_w / 2) + (dy / 2) + dy * (ii - 1), - (z_w / 2) + (dz / 2) + dz * (jj - 1)];
            [ch.pos] = pm_charge_calc_rot_trans_pos(ch.pos_block, rot, trans);
            ch.x = ch.pos(1);
            ch.y = ch.pos(2);
            ch.z = ch.pos(3);
            ch.block_num = block_num;
            ch_num = ch_num + 1;
            charge(ch_num) = ch;
        end
    end
end

%Y+ face
for ii = 1 : num_zdiv
    for jj = 1 : num_xdiv
        q = my * da_y;
        if abs(q) > 0
            ch.q = q;
            ch.pos_block = [- (x_w / 2) + (dx / 2) + dx * (jj - 1), (y_w / 2), - (z_w / 2) + (dz / 2) + dz * (ii - 1)];
            [ch.pos] = pm_charge_calc_rot_trans_pos(ch.pos_block, rot, trans);
            ch.x = ch.pos(1);
            ch.y = ch.pos(2);
            ch.z = ch.pos(3);
            ch.block_num = block_num;
            ch_num = ch_num + 1;
            charge(ch_num) = ch;
        end
    end
end

%Y- face
for ii = 1 : num_zdiv
    for jj = 1 : num_xdiv
        q = -my * da_y;
        if abs(q) > 0
            ch.q = q;
            ch.pos_block = [- (x_w / 2) + (dx / 2) + dx * (jj - 1), - (y_w / 2), - (z_w / 2) + (dz / 2) + dz * (ii - 1)];
            [ch.pos] = pm_charge_calc_rot_trans_pos(ch.pos_block, rot, trans);
            ch.x = ch.pos(1);
            ch.y = ch.pos(2);
            ch.z = ch.pos(3);
            ch.block_num = block_num;
            ch_num = ch_num + 1;
            charge(ch_num) = ch;
        end
    end
end

%Z+ face
for ii = 1 : num_xdiv
    for jj = 1 : num_ydiv
        q = mz * da_z;
        if abs(q) > 0
            ch.q = q;
            ch.pos_block = [- (x_w / 2) + (dx / 2) + dx * (ii - 1), - (y_w / 2) + (dy / 2) + dy * (jj - 1), (z_w / 2)];
            [ch.pos] = pm_charge_calc_rot_trans_pos(ch.pos_block, rot, trans);
            ch.x = ch.pos(1);
            ch.y = ch.pos(2);
            ch.z = ch.pos(3);
            ch.block_num = block_num;
            ch_num = ch_num + 1;
            charge(ch_num) = ch;
        end
    end
end

%Z- face
for ii = 1 : num_xdiv
    for jj = 1 : num_ydiv
        q = -mz * da_z;
        if abs(q) > 0
            ch.q = q;
            ch.pos_block = [- (x_w / 2) + (dx / 2) + dx * (ii - 1), - (y_w / 2) + (dy / 2) + dy * (jj - 1), - (z_w / 2)];
            [ch.pos] = pm_charge_calc_rot_trans_pos(ch.pos_block, rot, trans);
            ch.x = ch.pos(1);
            ch.y = ch.pos(2);
            ch.z = ch.pos(3);
            ch.block_num = block_num;
            ch_num = ch_num + 1;
            charge(ch_num) = ch;
        end
    end
end

%Delete unused pre-allocated charges
charge = charge(1 : ch_num);

%Done
        