%**************************************************************************
%This program calculates the field in a planar undulator using the
%technique of Walker, NIM A237 (1985), 366.
%
%Zachary Wolf, 10/12/09

%Clear the workspace
clear

%Set parameter values by executing the parameter file
calc_und_field_param

%Calculated parameters
l_u = 2 * lz_m + 2 * lz_p;
k_u = 2 * pi / l_u;
z = - l_u / 2 : .000001 : l_u / 2;

%Output files
dat_file = [OUTPUT_DIR, 'und_field.dat'];
plt_file = [OUTPUT_DIR, 'und_field.ps'];
mat_file = [OUTPUT_DIR, 'und_field.mat'];

%File header information
msg{1} = ['Name = ', NAME];
msg{2} = ['Output Directory = ', OUTPUT_DIR];
msg{3} = ['Br = ', num2str(b_r), ' T'];
msg{4} = ['Magnet Block Relative Permeability = ', num2str(mu_r)];
msg{5} = ['Magnet Width In X = ', num2str(lx_m), ' m'];
msg{6} = ['Magnet Height In Y = ', num2str(ly_m), ' m'];
msg{7} = ['Magnet Length In Z = ', num2str(lz_m), ' m'];
msg{8} = ['Pole Width In X = ', num2str(lx_p), ' m'];
msg{9} = ['Pole Height In Y = ', num2str(ly_p), ' m'];
msg{10} = ['Pole Length In Z = ', num2str(lz_p), ' m'];
msg{11} = ['Undulator Gap = ', num2str(g_u), ' m'];
msg{12} = ['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 Hm from the parameters
[h_m] = field_calc_calc_hm(b_r, mu_r, lx_m, ly_m, lz_m, lx_p, ly_p, lz_p, g_u);

%Calculate the field expansion coefficients
[bn] = field_calc_calc_bn(h_m, lz_m, lz_p, g_u, num_bn, plt_file);

%Replace with known expansion coefficients from an undulator, if desired
% bn(1) = 12474.9e-4;
% bn(2) = 0;
% bn(3) = -416.3e-4;
% bn(4) = 0;
% bn(5) = -108.0e-4;
% bn(6) = 0;
% bn(7) = 25.8e-4;
% bn(8) = 0;
% bn(9) = -1.8e-4;
% bn(10) = 0;

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

%Field calculation
[by_g2] = field_calc_calc_by(bn, k_u, g_u / 2, z);
[bz_g2] = field_calc_calc_bz(bn, k_u, g_u / 2, z);
[by] = field_calc_calc_by(bn, k_u, 0, z);

%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, bn);
field_calc_plot_rn(plt_file, rn);
field_calc_plot_byz_g2(plt_file, z, by_g2, bz_g2);
field_calc_plot_by(plt_file, z, by);

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

%Save the results to the mat file
save(mat_file, 'h_m', 'bn', 'rn', 'z', 'by', 'k');

%Done
