%**************************************************************************
%This function is used to study the effect of long coil straightness on EPU
%measurements.
%
%Zachary Wolf, 4/1/13

function [] = long_coil_straight()

%Clear the workspace
clear

%Set parameter values by executing the parameter file
param = long_coil_straight_param();

%Output files
dat_file = [param.output_file_dir, 'long_coil_straight.dat'];
plt_file = [param.output_file_dir, 'long_coil_straight.ps'];
mat_file = [param.output_file_dir, 'long_coil_straight.mat'];

%Initialize the data file
util_dat_init(dat_file);

%Initialize the plot file
msg = {' ', ' ', 'Long coil before straightening', 'Effect on field integral in calculated field'};
util_plot_init(plt_file, msg);

%Read the long coil straightness measurements
[header, y_meas, z_meas] = long_coil_straight_read_coil_meas(param.coil_meas_dat_file);

%Take out the average y-position of the coil
p = polyfit(z_meas, y_meas, 1);
y_fit = polyval(p, z_meas);
y_coil = y_meas - y_fit;
z_coil = z_meas;

%Put in coil sag
%y_coil = y_coil - 200e-6;

%Take out peak in center of coil
z1 = -.055;
z2 = .115;
index_array = find((z_coil > z1) & (z_coil < z2));
y1 = y_coil(index_array(1));
y2 = y_coil(index_array(end));
y_coil(index_array) = y1 + ((y2 - y1) / (z2 - z1)) * (z_coil(index_array) - z1);

%Plot the data
long_coil_straight_plot_coil_meas(plt_file, y_meas, z_meas);
long_coil_straight_plot_fit_coil_meas(plt_file, y_coil, z_coil);

%Set undulator polarization mode
%1 = lp_vf, 2 = lp_hf, 3 = cp_rh, 4 = cp_lh
und_mode = 3;

%Set up parameters for the undulator fields
lambda_u = 0.15;
num_per = 16;
kz = 2 * pi /lambda_u;
if und_mode == 1
    %lp_vf---
    kx = 0;
    ky = kz;
    phi0 = 1.5 / ky;
elseif und_mode == 2
    %lp_hf---
    kx = kz;
    ky = 0;
    phi0 = 1.5 / kx;
elseif und_mode == 3 || und_mode == 4
    %cp
    kx = kz / sqrt(2);
    ky = kz / sqrt(2);
    phi0 = 1.5 / kx;
end
xc = 0;
yc = 0;
zc = 0;
x = -0.015;
y = y_coil;
%y = zeros(1, length(y_coil));
z = z_coil;

%Test
%x = -0.015 + .0005 * cos(kz * z);

%Get the fields at the coil measurement locations
if und_mode == 1, [bx, by, bz] = lp_vf_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end
if und_mode == 2, [bx, by, bz] = lp_hf_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end
%if und_mode == 3, [bx, by, bz] = cp_rh_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end
if und_mode == 3, [bx, by, bz] = cp_rh_get_field(num_per, xc, yc, zc, x, y, z); end
if und_mode == 4, [bx, by, bz] = cp_lh_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end

%Plot the fields
field_plot_bx_vs_z(plt_file, z, bx);
field_plot_by_vs_z(plt_file, z, by);
field_plot_bz_vs_z(plt_file, z, bz);

%Integrate the fields
i1x = int_calc_int_y_dx(z, bx);
i1y = int_calc_int_y_dx(z, by);
i1z = int_calc_int_y_dx(z, bz);

%Display the integrals
i1x(end)
i1y(end)
i1z(end)

%Loop over zc
zc_save = zc;
%zc_vals = -.08 : .001 : .08;
zc_vals = -.08 : .01 : .08;
for ii = 1 : length(zc_vals)
    zc = zc_vals(ii);
    if und_mode == 1, [bx, by, bz] = lp_vf_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end
    if und_mode == 2, [bx, by, bz] = lp_hf_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end
    %if und_mode == 3, [bx, by, bz] = cp_rh_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end
    if und_mode == 3, [bx, by, bz] = cp_rh_get_field(num_per, xc, yc, zc, x, y, z); end
    if und_mode == 4, [bx, by, bz] = cp_lh_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end
    i1x = int_calc_int_y_dx(z, bx);
    i1y = int_calc_int_y_dx(z, by);
    i1z = int_calc_int_y_dx(z, bz);
    i1x_vs_zc(ii) = i1x(end);
    i1y_vs_zc(ii) = i1y(end);
    i1z_vs_zc(ii) = i1z(end);
end
figure;
hold on;
plot(zc_vals, i1x_vs_zc, 'linewidth', 0.5);
title('I1x vs Zc');
xlabel('Zc (m)');
ylabel('I1x (Tm)');
grid on;
zoom on;
hold off;
orient landscape;
print('-append', '-dpsc', plt_file);
figure;
hold on;
plot(zc_vals, i1y_vs_zc, 'linewidth', 0.5);
title('I1y vs Zc');
xlabel('Zc (m)');
ylabel('I1y (Tm)');
grid on;
zoom on;
hold off;
orient landscape;
print('-append', '-dpsc', plt_file);
zc = zc_save;

%Loop over x
x_save = x;
%x_vals = -.05 : .001 : .05;
x_vals = -.05 : .1 : .05;
for ii = 1 : length(x_vals)
    x = x_vals(ii);
    if und_mode == 1, [bx, by, bz] = lp_vf_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end
    if und_mode == 2, [bx, by, bz] = lp_hf_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end
    %if und_mode == 3, [bx, by, bz] = cp_rh_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end
    if und_mode == 3, [bx, by, bz] = cp_rh_get_field(num_per, xc, yc, zc, x, y, z); end
    if und_mode == 4, [bx, by, bz] = cp_lh_calc_field(lambda_u, num_per, phi0, kx, ky, kz, xc, yc, zc, x, y, z); end
    i1x = int_calc_int_y_dx(z, bx);
    i1y = int_calc_int_y_dx(z, by);
    i1z = int_calc_int_y_dx(z, bz);
    i1x_vs_x(ii) = i1x(end);
    i1y_vs_x(ii) = i1y(end);
    i1z_vs_x(ii) = i1z(end);
end
figure;
hold on;
plot(x_vals, i1x_vs_x, 'linewidth', 0.5);
title('I1x vs X');
xlabel('X (m)');
ylabel('I1x (Tm)');
grid on;
zoom on;
hold off;
orient landscape;
print('-append', '-dpsc', plt_file);
figure;
hold on;
plot(x_vals, i1y_vs_x, 'linewidth', 0.5);
title('I1y vs X');
xlabel('X (m)');
ylabel('I1y (Tm)');
grid on;
zoom on;
hold off;
orient landscape;
print('-append', '-dpsc', plt_file);
x = x_save;

%Save all data
save(mat_file, 'param', 'header');

%Done