%**************************************************************************
%This program calculates the field in a pure permanent magnet planar
%undulator using the technique of Halbach, NIM 187 (1981), 109.
%
%Zachary Wolf, 4/26/11

%Clear the workspace
clear

%Set parameter values by executing the parameter file
calc_ppm_und_field_param

%Calculated parameters
l_u = 4 * lz_m;
k_u = 2 * pi / l_u;
z = - l_u / 2 : .000001 : l_u / 2;
y = 0;
z_c = z + 1i * y;

%Output files
dat_file = [OUTPUT_DIR, 'ppm_und_field.dat'];
plt_file = [OUTPUT_DIR, 'ppm_und_field.ps'];

%File header information
msg{1} = ['Name = ', NAME];
msg{2} = ['Output Directory = ', OUTPUT_DIR];
msg{3} = ['Br = ', num2str(b_r), ' T'];
msg{4} = ['Magnet Width In X = ', num2str(lx_m), ' m'];
msg{5} = ['Magnet Height In Y = ', num2str(ly_m), ' m'];
msg{6} = ['Magnet Length In Z = ', num2str(lz_m), ' m'];
msg{7} = ['Undulator Gap = ', num2str(g_u), ' m'];
msg{8} = ['Number Of Terms In Field Expansion = ', num2str(num_bn)];

%Initialize the data file
util_dat_init(dat_file);
util_dat_msg(dat_file, cellstr(' '));
util_dat_msg(dat_file, cellstr('Parameters:'));
util_dat_msg(dat_file, msg);

%Initialize the plots
util_plot_init(plt_file, msg);

%Calculate the field expansion terms, eq. 5 of Halbach paper
bn = zeros(1, num_bn);
b_r_c = b_r;
k = 2 * pi / l_u;
h = g_u / 2;
l = ly_m;
epsilon = 1;
m = 4;
for ii = 0 : num_bn - 1
    n = 1 + ii * m;
    if n > num_bn, break; end;
    bn(n) = 1i * 2 * b_r_c * exp(-n * k * h) * (1 - exp(-n * k * l)) * sin(n * epsilon * pi / m) / (n * pi / m);
end

%Calculate ratio of expansion coefficients
rn = abs(bn) / abs(bn(1));

%Field calculation
b_c_star = 0;
for ii = 1 : num_bn
    b_c_star = b_c_star + bn(ii) * cos(ii * k * z_c);
end
bz = real(b_c_star);
by = - imag(b_c_star);

%Calculate the peak field
by_peak = max(by);

%Calculate K
[k] = field_calc_calc_k(z, by);
          
%Plot the results of the field calculation
field_calc_plot_bn(plt_file, abs(bn));
field_calc_plot_rn(plt_file, rn);
field_calc_plot_by(plt_file, z, by);

%Write the results of the field calculation to the data file
field_calc_dat_bn(dat_file, abs(bn), abs(rn));
field_calc_dat_by_peak(dat_file, by_peak);
field_calc_dat_k(dat_file, k);

%Done
