%**************************************************************************
%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
%   num_block_xdiv  %Number of divisions in x for charge calculation
%   num_block_ydiv  %Number of divisions in y for charge calculation
%   num_block_zdiv  %Number of divisions in z for charge calculation
%   und_gap         %Undulator gap height (m)
%   period          %Undulator period (m)
%   num_ctr_period  %Number of periods in the undulator, excluding the ends
%   br              %Permanent magnet remnant field (T)
%   block_x_width   %Magnet block width in x direction (m)
%   block_y_width   %Magnet block width in y direction (m)
%   block_z_width   %Magnet block width in z direction (m)
%   block_end_mag_1_z_width 	%Width of end magnet block (m)
%   block_end_space_2_z_width	%Width of space between 1 and 2 end blocks (m)
%   block_end_mag_2_z_width 	%Width of block second from the end (m)
%   block_end_space_3_z_width	%Width of space between 2 and 3 end blocks (m)
%   block_end_mag_3_z_width     %Width of block third from the end (m)
%   x_field_pts     %X positions where field is calculated (m)
%   y_field_pts     %Y positions where field is calculated (m)
%   z_field_pts     %Z positions where field is calculated (m)
%
%block, structure containing the permanent magnet block information
%   num_block       %Number of blocks
%   mx              %X component of the magnetic moment density (A/m)
%   my              %Y component of the magnetic moment density (A/m)
%   mz              %Z component of the magnetic moment density (A/m)
%   x_center        %X position of the block center (m)
%   y_center        %Y position of the block center (m)
%   z_center        %Z position of the block center (m)
%   x_width         %Block width in the x direction (m)
%   y_width         %Block width in the y direction (m)
%   z_width         %Block width in the z direction (m)
%
%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
%
%field_pos, structure giving the positions at which to calculate the field
%   x               %X position of the field point (m)
%   y               %Y position of the field point (m)
%   z               %Z position of the field point (m)
%
%mag_field, structure containing the calculated magnetic field
%   x               %X position of the field point (m)
%   y               %Y position of the field point (m)
%   z               %Z position of the field point (m)
%   bx              %Bx magnetic field component at (x, y, z) (T)
%   by              %By magnetic field component at (x, y, z) (T)
%   bz              %Bz magnetic field component at (x, y, z) (T)
%
%Zachary Wolf, 3/31/08

%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)

%Program parameters
[param] = ppm_und_sim_get_param();

%Initialize output files
[dat_file, plt_file] = ppm_und_sim_file_init(param);

%Initialize the permanent magnet blocks
disp('Initializing blocks...');
[block] = ppm_und_sim_block_init(param);

%Add an error to block 39
%block.my(39) = block.my(39) * 1.01;
%block.y_center(39) = block.y_center(39) + 300.e-6;

%Add an error to block 41
%block.my(41) = block.my(41) * 1.01;
%block.y_center(39) = block.y_center(39) + 300.e-6;

%Rotate the magnetization in block 39
%block.mx(39) = - block.my(39) * sin(1 * pi / 180);  %Rotate 1 degree about z-axis
%block.my(39) = block.my(39) * cos(1 * pi / 180);  %Rotate 1 degree about z-axis

%Rotate horizontally the magnetization in block 41
%block.mx(41) = block.mz(41) * sin(1 * pi / 180);  %Rotate 1 degree about y-axis
%block.mz(41) = block.mz(41) * cos(1 * pi / 180);  %Rotate 1 degree about y-axis

%Rotate vertically the magnetization in block 41
%block.my(41) = block.mz(41) * sin(1 * pi / 180);  %Rotate 1 degree about y-axis
%block.mz(41) = block.mz(41) * cos(1 * pi / 180);  %Rotate 1 degree about y-axis

%Move blocks 39 and 40, get slope sensitivity, V blocks
%block.y_center(39) = block.y_center(39) + -100.e-6;
%block.y_center(40) = block.y_center(40) - -100.e-6;

%Add an error to pole 6, block 19
%block.my(19) = block.my(19) * 1.01;
%block.y_center(19) = block.y_center(19) + 300.e-6;

%Move blocks 19 and 20, get slope sensitivity, V blocks
%block.y_center(19) = block.y_center(19) + 150.e-6;
%block.y_center(20) = block.y_center(20) - 150.e-6;

%Move blocks 21 and 22, H blocks
%block.y_center(21) = block.y_center(21) + 300.e-6;
%block.y_center(22) = block.y_center(22) - 300.e-6;

%Move blocks 19, 20, 23, and 24, phase shift
%block.y_center(19) = block.y_center(19) + 600.e-6;
%block.y_center(20) = block.y_center(20) - 600.e-6;
%block.y_center(23) = block.y_center(23) + 600.e-6;
%block.y_center(24) = block.y_center(24) - 600.e-6;

%Move blocks 39, 40, 43, and 44, phase shift
%block.y_center(39) = block.y_center(39) + -200.e-6;
%block.y_center(40) = block.y_center(40) - -200.e-6;
%block.y_center(43) = block.y_center(43) + -200.e-6;
%block.y_center(44) = block.y_center(44) - -200.e-6;

%Add a 50 micron taper
%for ii = 1 : length(block.y_center)
%    block.y_center(ii) = block.y_center(ii) * (1. + (50.e-6 / 2) * (block.z_center(ii) - block.z_center(1)) / (block.z_center(end) - block.z_center(1)));
%end

%Add a 0 micron taper
%for ii = 1 : length(block.y_center)
%    block.y_center(ii) = block.y_center(ii) * (1. + (0.e-6 / 2) * (block.z_center(ii) - block.z_center(1)) / (block.z_center(end) - block.z_center(1)));
%end

%Calculate charges on blocks
disp('Calculating charges...');
[charge] = ppm_calc_charge(param, block);

%Rotate block 39
%index_array = find(charge.block_num == 39);
%r_charge.num_charge = charge.num_charge;
%r_charge.q = charge.q(index_array);
%r_charge.x = charge.x(index_array);
%r_charge.y = charge.y(index_array);
%r_charge.z = charge.z(index_array);
%r_charge.block_num = charge.block_num(index_array);
%[rotz_charge] = ppm_calc_rotz_charge(r_charge, block.x_center(39), block.y_center(39) - block.y_width(39) / 2, -0.25 * pi / 180);
%Translate
%rotz_charge.x = rotz_charge.x - block.y_width(39) * sin(2 * pi / 180);
%rotz_charge.y = rotz_charge.y - block.y_width(39) * (1 - cos(2 * pi / 180));
%End translate
%charge.q(index_array) = rotz_charge.q;
%charge.x(index_array) = rotz_charge.x;
%charge.y(index_array) = rotz_charge.y;
%charge.z(index_array) = rotz_charge.z;
%charge.block_num(index_array) = rotz_charge.block_num;

%Rotate block 40
%index_array = find(charge.block_num == 40);
%r_charge.num_charge = charge.num_charge;
%r_charge.q = charge.q(index_array);
%r_charge.x = charge.x(index_array);
%r_charge.y = charge.y(index_array);
%r_charge.z = charge.z(index_array);
%r_charge.block_num = charge.block_num(index_array);
%[rotz_charge] = ppm_calc_rotz_charge(r_charge, block.x_center(40), block.y_center(40) + block.y_width(40) / 2, -0.25 * pi / 180);
%charge.q(index_array) = rotz_charge.q;
%charge.x(index_array) = rotz_charge.x;
%charge.y(index_array) = rotz_charge.y;
%charge.z(index_array) = rotz_charge.z;
%charge.block_num(index_array) = rotz_charge.block_num;

%Calculate the field at the center of block 39
%param.x_field_pts = -.03 : .001 : .03;
%param.y_field_pts = 0;
%param.z_field_pts = block.z_center(39);

%Turn off all blocks except 39 and 40
%index_array = find((charge.block_num ~= 39) & (charge.block_num ~= 40));
%charge.q(index_array) = 0;

%Calculate fields from the charges
disp('Calculating fields...');
[mag_field] = ppm_calc_mag_field(param, charge);

%Add a field difference between the lab and the tunnel
%zm = mag_field.z;
%index_array = find((zm >= 0) & (zm <= .5333));
%mag_field.bx(index_array) = mag_field.bx(index_array) + 5.e-5;  %0.5 G
%mag_field.by(index_array) = mag_field.by(index_array) + 5.e-5;
%index_array = find((zm < 0) | (zm > .5333));
%mag_field.bx(index_array) = mag_field.bx(index_array) + 2.5e-5;  %0.25 G
%mag_field.by(index_array) = mag_field.by(index_array) + 2.5e-5;

%Plot the field
ppm_plot_block_charge_yz(plt_file, block, charge);
ppm_plot_b_on_ctr_z_line(plt_file, block, charge, mag_field);
ppm_plot_b_on_z_line(plt_file, block, charge, mag_field, 0, 0);
ppm_plot_b_on_x_slice(plt_file, block, charge, mag_field, 0);
ppm_plot_by_vs_z(plt_file, block, charge, mag_field, 0, 0);
ppm_plot_bz_vs_z(plt_file, block, charge, mag_field, 0, .005);
ppm_plot_bx_vs_z(plt_file, block, charge, mag_field, 0, 0);
ppm_plot_bx_vs_x(plt_file, block, charge, mag_field, 0, block.z_center(39));
ppm_plot_by_vs_x(plt_file, block, charge, mag_field, 0, block.z_center(39));

%Write the results to a data file
ppm_dat_byx_vs_z(dat_file, mag_field, 0, 0);

%Output the results to the output file directory
mat_file = strcat(lower(param.output_file_dir), 'ppm_und_sim.mat');
save(mat_file, 'param', 'block', 'charge', 'mag_field');

%Done
disp('Done');