%**************************************************************************
%ppm_und_block_init_planar
%This file initializes the structure containing the permanent magnet
%block information.
%
%Input:
%param, parameter structure
%
%Output:
%block(i), structure array containing the rectangular permanent magnet block information
%   num                 %Number to identify the block
%   div_size            %Magnet block segmentation size along each axis for equivalent charges (m)
%   br                  %Permanent magnet remnant field (T)
%   m(3)                %Magnetic moment per unit volume in coordinate system tied to block [m_x, m_y, m_z] (A/m)
%   dim(3)              %Magnet block dimension in each direction in coordinate system tied to block [w_x, w_y, w_z] (m)
%   rot(3x3)            %Rotation matrix, rotate the block about its center, direction cosines R_ij = e_rot_i . e_orig_j
%   trans(3)            %Translate the rotated block in x, y, and z in global coordinates parallel to the original block coordinates [t_x, t_y, t_z] (m)
%
%Zachary Wolf, 4/22/10

function [block] = ppm_und_block_init_planar(param)

%Global
global mu_0;

%Undulator parameters
und_gap = param.und_gap;
und_period = param.und_period;

%General pm block parameters
div_size = param.pm_div_size;

%End pm block parameters
end_pm_num_z_pos = param.end_pm_num_z_pos;
end_pm1_dir = param.end_pm1_dir;
end_pm_br = param.end_pm_br;
end_pm_dim = param.end_pm_dim;
end_spacer_dim = param.end_spacer_dim;

%Center pm block parameters
ctr_pm_num_z_pos = param.ctr_pm_num_z_pos;
ctr_pm_br = param.ctr_pm_br;
ctr_pm_dim = param.ctr_pm_dim;
ctr_spacer_dim = param.ctr_spacer_dim;

%Determine the number of block z positions along the undulator
num_z_pos = end_pm_num_z_pos + ctr_pm_num_z_pos + end_pm_num_z_pos;

%Determine the block x positions
x_pos_up = zeros(1, num_z_pos);
x_pos_lo = zeros(1, num_z_pos);

%Determine the block y positions
y_pos_up = zeros(1, num_z_pos);
y_pos_lo = zeros(1, num_z_pos);
nn = 0;
for ii = 1 : end_pm_num_z_pos
    nn = nn + 1;
    y_pos_up(nn) =  (und_gap / 2) + (end_pm_dim(ii, 2) / 2);
    y_pos_lo(nn) = -(und_gap / 2) - (end_pm_dim(ii, 2) / 2);
end
for ii = 1 : ctr_pm_num_z_pos
    nn = nn + 1;
    y_pos_up(nn) =  (und_gap / 2) + (ctr_pm_dim(2) / 2);
    y_pos_lo(nn) = -(und_gap / 2) - (ctr_pm_dim(2) / 2);
end
for ii = end_pm_num_z_pos : -1 : 1
    nn = nn + 1;
    y_pos_up(nn) =  (und_gap / 2) + (end_pm_dim(ii, 2) / 2);
    y_pos_lo(nn) = -(und_gap / 2) - (end_pm_dim(ii, 2) / 2);
end

%Determine the block z positions
z_pos = zeros(1, num_z_pos);
nn = 1;
z_pos(nn) = end_pm_dim(1, 3) / 2;
z_ref = 0 + end_pm_dim(1, 3) + end_spacer_dim(1);
for ii = 2 : end_pm_num_z_pos
    nn = nn + 1;
    z_pos(nn) = z_ref + end_pm_dim(ii, 3) / 2;
    z_ref = z_ref + end_pm_dim(ii, 3) + end_spacer_dim(ii);
end
for ii = 1 : ctr_pm_num_z_pos
    nn = nn + 1;
    z_pos(nn) = z_ref + ctr_pm_dim(3) / 2;
    z_ref = z_ref + ctr_pm_dim(3) + ctr_spacer_dim;
end
z_ref = z_ref - ctr_spacer_dim;  %Use end spacer instead, now ref goes to block edge
for ii = end_pm_num_z_pos : -1 : 1
    nn = nn + 1;
    z_pos(nn) = z_ref + end_spacer_dim(ii) + end_pm_dim(ii, 3) / 2;
    z_ref = z_ref + end_spacer_dim(ii) + end_pm_dim(ii, 3);
end
z_pos_up = z_pos;
z_pos_lo = z_pos;

%Determine the block remnant fields
br = zeros(1, num_z_pos);
nn = 0;
for ii = 1 : end_pm_num_z_pos
    nn = nn + 1;
    br(nn) = end_pm_br(ii);
end
for ii = 1 : ctr_pm_num_z_pos
    nn = nn + 1;
    br(nn) = ctr_pm_br;
end
for ii = end_pm_num_z_pos : -1 : 1
    nn = nn + 1;
    br(nn) = end_pm_br(ii);
end
br_up = br;
br_lo = br;

%Set up magnetization constants
%Convention: up = 1, right = 2, down = 3, left = 4
m_vec(1, :) = [0, 1, 0];   %up
m_vec(2, :) = [0, 0, 1];   %right
m_vec(3, :) = [0, -1, 0];  %down
m_vec(4, :) = [0, 0, -1];  %left

%Determine the block magnetization directions
%Convention: up = 1, right = 2, down = 3, left = 4
%Once the first block is specified, the rest are determined
m_dir_up = zeros(1, num_z_pos);
m_dir_lo = zeros(1, num_z_pos);
m_dir_up(1) = end_pm1_dir;
for ii = 2 : num_z_pos
    m_dir_up(ii) = m_dir_up(ii - 1) + 1;
    if m_dir_up(ii) == 5, m_dir_up(ii) = 1; end;
end
if m_dir_up(1) == 1, m_dir_lo(1) = 1; end;
if m_dir_up(1) == 2, m_dir_lo(1) = 4; end;
if m_dir_up(1) == 3, m_dir_lo(1) = 3; end;
if m_dir_up(1) == 4, m_dir_lo(1) = 2; end;
for ii = 2 : num_z_pos
    m_dir_lo(ii) = m_dir_lo(ii - 1) - 1;
    if m_dir_lo(ii) == 0, m_dir_lo(ii) = 4; end;
end

%Determine the block magnetization vectors
m_up = zeros(num_z_pos, 3);
m_lo = zeros(num_z_pos, 3);
for ii = 1 : num_z_pos
    m_up(ii, :) = (br_up(ii) / mu_0) * m_vec(m_dir_up(ii), :);
    m_lo(ii, :) = (br_lo(ii) / mu_0) * m_vec(m_dir_lo(ii), :);
end

%Store the block dimensions
dim = zeros(num_z_pos, 3);
nn = 0;
for ii = 1 : end_pm_num_z_pos
    nn = nn + 1;
    dim(nn, :) = end_pm_dim(ii, :);
end
for ii = 1 : ctr_pm_num_z_pos
    nn = nn + 1;
    dim(nn, :) = ctr_pm_dim;
end
for ii = end_pm_num_z_pos : -1 : 1
    nn = nn + 1;
    dim(nn, :) = end_pm_dim(ii, :);
end

%Fill the block structure
block(2 * num_z_pos) = 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 : num_z_pos  %Upper blocks
    nn = nn + 1;
    block(nn).num = nn;
    block(nn).div_size = div_size;
    block(nn).br = br(ii);
    block(nn).m = m_up(ii, :);
    block(nn).dim = dim(ii, :);
    block(nn).rot = [1, 0, 0; 0, 1, 0; 0, 0, 1];
    block(nn).trans = [x_pos_up(ii), y_pos_up(ii), z_pos_up(ii)];
end
for ii = 1 : num_z_pos  %Lower blocks
    nn = nn + 1;
    block(nn).num = nn;
    block(nn).div_size = div_size;
    block(nn).br = br(ii);
    block(nn).m = m_lo(ii, :);
    block(nn).dim = dim(ii, :);
    block(nn).rot = [1, 0, 0; 0, 1, 0; 0, 0, 1];
    block(nn).trans = [x_pos_lo(ii), y_pos_lo(ii), z_pos_lo(ii)];
end

%Done
