%**************************************************************************
%This file initializes the simulation program.
%
%Input:
%param, parameter structure
%
%Output:
%block, structure containing permanent magnet block information
%
%Zachary Wolf, 2/4/11

function [block] = virtual_shim_block_init(param)

%Global
global mu_0;

%Parameters
pm_br = param.pm_br;
dim = param.pm_dim;
x_nom_pos = param.pm_x_nom_pos;
y_nom_pos = param.pm_y_nom_pos;
z_nom_pos = param.pm_z_nom_pos;
m_dir = param.pm_m_dir;
div_size = param.pm_div_size;
x_shift = param.shim_x_shift;
y_shift = param.shim_y_shift;
z_shift = param.shim_z_shift;

%Determine the number of blocks
num_actual_blocks = length(x_nom_pos);
num_blocks = 2 * num_actual_blocks;

%Determine the x positions of the shifted blocks
x_pos = x_nom_pos + x_shift;
%Add the x positions of the un-shifted blocks
x_pos = [x_pos, x_nom_pos];

%Determine the y positions of the shifted blocks
y_pos = y_nom_pos + y_shift;
%Add the y positions of the un-shifted blocks
y_pos = [y_pos, y_nom_pos];

%Determine the z positions of the shifted blocks
z_pos = z_nom_pos + z_shift;
%Add the z positions of the un-shifted blocks
z_pos = [z_pos, z_nom_pos];

%Determine the block remnant fields
br = pm_br * ones(1, num_blocks);

%Set up magnetization constants
%Convention: +x = 1, -x = 2, +y = 3, -y = 4, +z = 5, -z = 6
m_vec(1, :) = [1, 0, 0];   %+x
m_vec(2, :) = [-1, 0, 0];  %-x
m_vec(3, :) = [0, 1, 0];   %+y
m_vec(4, :) = [0, -1, 0];  %-y
m_vec(5, :) = [0, 0, 1];   %+z
m_vec(6, :) = [0, 0, -1];  %-z

%Reverse the magnetization in the non-shifted blocks
m_rev_dir = zeros(1, length(m_dir));
for ii = 1 : length(m_dir)
    if m_dir(ii) == 1, m_rev_dir(ii) = 2; end
    if m_dir(ii) == 2, m_rev_dir(ii) = 1; end
    if m_dir(ii) == 3, m_rev_dir(ii) = 4; end
    if m_dir(ii) == 4, m_rev_dir(ii) = 3; end
    if m_dir(ii) == 5, m_rev_dir(ii) = 6; end
    if m_dir(ii) == 6, m_rev_dir(ii) = 5; end
end

%Determine the block magnetization vectors
m = zeros(num_blocks, 3);
nn = 0;
for ii = 1 : num_actual_blocks
    nn = nn + 1;
    m(nn, :) = (pm_br / mu_0) * m_vec(m_dir(ii), :);
end
for ii = 1 : num_actual_blocks
    nn = nn + 1;
    m(nn, :) = (pm_br / mu_0) * m_vec(m_rev_dir(ii), :);
end

%Fill the block structure
block(2 * num_actual_blocks) = struct('num', 0, 'div_size', 0, 'br', 0, 'm', [0, 0, 0,], 'dim', [0, 0, 0], 'rot', [0, 0, 0; 0, 0, 0; 0, 0, 0], 'trans', [0, 0, 0]);
nn = 0;
for ii = 1 : 2 * num_actual_blocks 
    nn = nn + 1;
    block(nn).num = nn;
    block(nn).div_size = div_size;
    block(nn).br = br(ii);
    block(nn).m = m(ii, :);
    block(nn).dim = dim;
    block(nn).rot = [1, 0, 0; 0, 1, 0; 0, 0, 1];
    block(nn).trans = [x_pos(ii), y_pos(ii), z_pos(ii)];
end

%Done
end
