%**************************************************************************
%ppm_und_sim
%This program simulates the magnetic field in a pure permanent magnet
%undulator.  Magnet block parameters are entered.  Equivalent charges are
%calculated to replace the magnet blocks.  The field from each charge is
%calculated and added to the field from the other charges.
%
%The following structures are used.
%
%param, structure containing program parameters
%   output_file_dir     %Directory for output files
%   und_gap             %Undulator gap (m)
%   und_period          %Undulator period (m)
%   pm_div_size         %Permanent magnet segmentation size along each axis for equivalent charges (m)
%   end_pm1_dir         %Direction of the field in the first upper block, convention: up = 1, right = 2, down = 3, left = 4
%   end_pm_num_z_pos    %Number of z positions (ie. upper blocks) in the undulator ends
%   end_pm_br(n)        %Remnant field in the end blocks (T)
%   end_pm_dim(n, 3)    %Magnet block dimensions for each end block in coordinate system tied to block [w_x, w_y, w_z] (m) 
%   end_spacer_dim(n)   %Thickness of the spacer after each end block (m)
%   ctr_pm_num_z_pos    %Number of z positions (ie. upper blocks) in the undulator center
%   ctr_pm_br           %Remnant field in the center blocks (T)
%   ctr_pm_dim(3)       %Magnet block dimensions for each center block in coordinate system tied to block [w_x, w_y, w_z] (m) 
%   ctr_spacer_dim      %Thickness of the spacer after each center block (m)
%   field_pos(n, 3)     %Positions where field is calculated, [x(n), y(n), z(n)] coordinates (m)
%
%param_epu, structure containing program parameters
%   output_file_dir     %Directory for output files
%   end_pm1_dir         %Direction of the field in the first upper block, convention: up = 1, right = 2, down = 3, left = 4
%   end_pm_num_z_pos    %Number of z positions (ie. upper blocks) in the undulator ends
%   end_pm_br(n)        %Remnant field in the end blocks (T)
%   end_pm_dim(n, 3)    %Magnet block dimensions for each end block in coordinate system tied to block [w_x, w_y, w_z] (m) 
%   end_spacer_dim(n)   %Thickness of the spacer after each end block (m)
%   ctr_pm_num_z_pos    %Number of z positions (ie. upper blocks) in the undulator center
%   ctr_pm_br           %Remnant field in the center blocks (T)
%   ctr_pm_dim(3)       %Magnet block dimensions for each center block in coordinate system tied to block [w_x, w_y, w_z] (m) 
%   ctr_spacer_dim      %Thickness of the spacer after each center block (m)
%   und_gap             %Undulator gap (m)
%   und_h_gap           %Horizontal gap between upper and lower EPU quadrants (m)
%   und_period          %Undulator period (m)
%   und_q1_z_shift      %Shift along z of quadrant 1 (m)
%   und_q2_z_shift      %Shift along z of quadrant 2 (m)
%   und_q3_z_shift      %Shift along z of quadrant 3 (m)
%   und_q4_z_shift      %Shift along z of quadrant 4 (m)
%   pm_div_size         %Permanent magnet segmentation size along each axis for equivalent charges (m)
%   field_pos(n, 3)     %Positions where field is calculated, [x(n), y(n), z(n)] coordinates (m)
%
%param_varphase, structure containing program parameters
%   output_file_dir     %Directory for output files
%   end_pm1_dir         %Direction of the field in the first upper block, convention: up = 1, right = 2, down = 3, left = 4
%   end_pm_num_z_pos    %Number of z positions (ie. upper blocks) in the undulator ends
%   end_pm_br(n)        %Remnant field in the end blocks (T)
%   end_pm_dim(n, 3)    %Magnet block dimensions for each end block in coordinate system tied to block [w_x, w_y, w_z] (m) 
%   end_spacer_dim(n)   %Thickness of the spacer after each end block (m)
%   ctr_pm_num_z_pos    %Number of z positions (ie. upper blocks) in the undulator center
%   ctr_pm_br           %Remnant field in the center blocks (T)
%   ctr_pm_dim(3)       %Magnet block dimensions for each center block in coordinate system tied to block [w_x, w_y, w_z] (m) 
%   ctr_spacer_dim      %Thickness of the spacer after each center block (m)
%   und_gap             %Undulator gap (m)
%   und_period          %Undulator period (m)
%   und_up_z_shift      %Shift along z of upper jaw (m)
%   und_lo_z_shift      %Shift along z of lower jaw (m)
%   pm_div_size         %Permanent magnet segmentation size along each axis for equivalent charges (m)
%   field_pos(n, 3)     %Positions where field is calculated, [x(n), y(n), z(n)] coordinates (m)
%
%block_rect(i), structure array containing the rectangular permanent magnet block information
%   num                 %Number to identify the block
%   div_size            %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)
%
%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
%
%field, structure containing the calculated magnetic field
%   pos(nx3)            %Position [x(n), y(n), z(n)] of calculated field in global coordinates (m)
%   b(nx3)              %Magnetic field [Bx(n), By(n), Bz(n)] at position [x(n), y(n), z(n)] (T)
%
%Zachary Wolf, 4/23/10

%Clear all variables
clear;

%Set the working directory
cd('c:\Simulation\PPM Und Sim');

%Global
global mu_0;
mu_0 = 4 * pi * 10^-7;  %(Tm/A)

%Set type of undulator
%UND_TYPE = 'VARGAP';
%UND_TYPE = 'EPU';
UND_TYPE = 'VARPHASE';

%Program parameters
if strcmp(UND_TYPE, 'VARGAP')
    [param] = ppm_und_sim_param();
elseif strcmp(UND_TYPE, 'EPU')
    [param] = ppm_und_sim_param_epu();
elseif strcmp(UND_TYPE, 'VARPHASE')
    [param] = ppm_und_sim_param_varphase();
end

%Make file names
dat_file = strcat(lower(param.output_file_dir), 'ppm_und_sim.dat');
plt_file = strcat(lower(param.output_file_dir), 'ppm_und_sim.ps');
mat_file = strcat(lower(param.output_file_dir), 'ppm_und_sim.mat');

%Initialize the data file
util_sim_dat_init(dat_file);

%Write a message to the data file
msg{1} = ' ';
msg{2} = 'PPM Undulator Simulation';
msg{3} = '========================';
util_sim_dat_msg(dat_file, msg);

%Write the parameter values to the data file
if strcmp(UND_TYPE, 'VARGAP')
    ppm_und_dat_param(dat_file, param);
elseif strcmp(UND_TYPE, 'EPU')
    ppm_und_dat_param_epu(dat_file, param);
elseif strcmp(UND_TYPE, 'VARPHASE')
    ppm_und_dat_param_varphase(dat_file, param);
end

%Write the magnet configuration to the data file
if strcmp(UND_TYPE, 'VARGAP')
elseif strcmp(UND_TYPE, 'EPU')
    ppm_und_dat_block_config_epu(dat_file, param);
elseif strcmp(UND_TYPE, 'VARPHASE')
    ppm_und_dat_block_config_varphase(dat_file, param);
end

%Initialize the plot file
util_sim_plt_init(plt_file, msg, param);

%Initialize the permanent magnet blocks
disp('Initializing blocks...');
%tic;
if strcmp(UND_TYPE, 'VARGAP')
elseif strcmp(UND_TYPE, 'EPU')
    [block] = ppm_und_block_init_epu(param);
elseif strcmp(UND_TYPE, 'VARPHASE')
    [block] = ppm_und_block_init_varphase(param);
end
%dt = toc;
%disp(['dt = ', num2str(dt)]);

%Add block motions ----------

% %Move blocks at z_pos to make gap smaller or larger
% z_pos = 234;
% vert_move = -300.e-6;
% num_z_pos = length(block) / 4;
% q1_block_num = z_pos;
% q2_block_num = q1_block_num + num_z_pos;
% q3_block_num = q2_block_num + num_z_pos;
% q4_block_num = q3_block_num + num_z_pos;
% block(q1_block_num).trans(2) = block(q1_block_num).trans(2) + vert_move;  %Q1
% block(q2_block_num).trans(2) = block(q2_block_num).trans(2) + vert_move;  %Q2
% block(q3_block_num).trans(2) = block(q3_block_num).trans(2) - vert_move;  %Q3
% block(q4_block_num).trans(2) = block(q4_block_num).trans(2) - vert_move;  %Q4

%Move blocks at z_pos, Q1 and Q2 toward +x by horiz_move, Q3 and Q4 toward -x by horiz_move
% z_pos = 14;
% horiz_move = 300.e-6;
% num_z_pos = length(block) / 4;
% q1_block_num = z_pos;
% q2_block_num = q1_block_num + num_z_pos;
% q3_block_num = q2_block_num + num_z_pos;
% q4_block_num = q3_block_num + num_z_pos;
% block(q1_block_num).trans(1) = block(q1_block_num).trans(1) + horiz_move;  %Q1
% block(q2_block_num).trans(1) = block(q2_block_num).trans(1) + horiz_move;  %Q2
% block(q3_block_num).trans(1) = block(q3_block_num).trans(1) - horiz_move;  %Q3
% block(q4_block_num).trans(1) = block(q4_block_num).trans(1) - horiz_move;  %Q4

%End add block motions ----------

%Calculate charges on blocks
disp('Calculating charges...');
%tic;
charge(length(block) * 2000) = struct('q', 0, 'pos_block', [0, 0, 0], 'pos', [0, 0, 0], 'x', 0, 'y', 0, 'z', 0, 'block_num', 0);
nn = 0;
for ii = 1 : length(block)
    [charge_block] = pm_block_rect_calc_charge(block(ii));
    for jj = 1 : length(charge_block)
        nn = nn + 1;
        charge(nn) = charge_block(jj);
    end
end
charge = charge(1 : nn);  %remove extra pre-allocated charges
%dt = toc;
%disp(['dt = ', num2str(dt)]);

%Calculate fields from the charges
disp('Calculating fields...');
%tic;
[field] = pm_field_calc_b(param.field_pos, charge);
%dt = toc;
%disp(['dt = ', num2str(dt)]);

%Plot the field
%ppm_und_plot_charge_vs_yz(plt_file, block, charge);
ppm_und_plot_bv_vs_yz(plt_file, block, charge, field, 0);
ppm_und_plot_bx_vs_z(plt_file, block, charge, field, 0, 0);
ppm_und_plot_by_vs_z(plt_file, block, charge, field, 0, 0);
ppm_und_plot_bxy_vs_z(plt_file, block, charge, field, 0, 0);
ppm_und_plot_bz_vs_z(plt_file, block, charge, field, 0, 0);
%ppm_und_plot_by_vs_x(plt_file, block, charge, field, 0, block(10).trans(3));
%ppm_und_plot_by_vs_y(plt_file, block, charge, field, 0, block(10).trans(3));
%ppm_und_plot_bx_vs_x(plt_file, block, charge, field, 0, block(10).trans(3));
%ppm_und_plot_bx_vs_y(plt_file, block, charge, field, 0, block(10).trans(3));

%Write the results to a data file
%ppm_und_dat_byx_vs_z(dat_file, field, 0, 0);
ppm_und_dat_bxyz_vs_z(dat_file, field, 0, 0);

%Output the results to the output file directory
save(mat_file, 'param', 'block', 'charge', 'field');

%Done
disp('Done');
