
%Clear the workspace
clear

%Load the data
ser_num = 9;
run_num = 21;
file_name = 'C:\Magdata\SN25 y-scans\bscany_dat.021';
[header, y, by] = idanal_read_meas_yscan(file_name);

%Fit the minimum
by_abs = abs(by);
index_array = find(by_abs < min(by_abs) + 0.5 * (max(by_abs) - min(by_abs)));
y_fit = y(index_array);
by_abs_fit = by_abs(index_array);
p = polyfit(y_fit, by_abs_fit, 4);
by_abs_fit_val = polyval(p, y_fit);

%Find the field minimum
[by_abs_fit_min, index_fit_min] = min(by_abs_fit_val);
y_min = y_fit(index_fit_min);

%Plot the data
figure;
hold on;
plot(y, by_abs, 'b');
plot(y_fit, by_abs_fit, 'r');
plot(y_min, by_abs_fit_min, 'r*');
xlabel('Y (m)');
ylabel('|By| (T)');
title(['SN', int2str(ser_num), ', Run ', int2str(run_num), ', Y-Scan Data']);
legend('Measurements', 'Fit');
msg_b = ['|By| min = ', num2str(by_abs_fit_min), ' T'];
msg_y = ['Y min = ', num2str(y_min * 1000), ' mm'];
text(y_min, by_abs_fit_min + .008, msg_b);
text(y_min, by_abs_fit_min + .006, msg_y);
grid on;
hold off;

%Calculate changes
dby_abs = by_abs - by_abs_fit_min;
dby_abs_by = dby_abs / by_abs_fit_min;

%Compare to cosh function
k = 2 * pi / .03;
dby_abs_by_nom = cosh(k * (y - y_min)) - 1;

%Plot the changes
figure;
hold on;
plot(y, dby_abs_by, 'b');
plot(y, dby_abs_by_nom, 'r');
xlabel('Y (m)');
ylabel('Change in |By| / |By|');
title(['SN', int2str(ser_num), ', Run ', int2str(run_num), ', Change in |By|']);
legend('Measurements', 'Cosh(k dy) - 1');
grid on;
hold off;

%Calculate fields at different y-positions from the expansion parameters----------

%Constants
q = -1.60217733e-19;    %Electron charge (C)
m = 9.1093897e-31;      %Electron mass (kg)
c = 2.99792458e8;       %Speed of light (m/s)
gamma = 1;  %Drops out of K calculation

%Put in known expansion coefficients from SN09
bn(1) = 12475.7e-4;
bn(2) = 0;
bn(3) = -410.0e-4;
bn(4) = 0;
bn(5) = -106.6e-4;
bn(6) = 0;
bn(7) = 25.5e-4;
bn(8) = 0;
bn(9) = -1.6e-4;
bn(10) = 0;

%Set the range of z-values and y-values
l_u = .030;  %undulator period (m)
z_calc = l_u / 4;
y_calc = 0 : 10.e-6 : .95e-3;

%Calculate By as a function of y
%by = Sum bn sin(nkz) cosh(nky)
%bz = Sum bn cos(nkz) sinh(nky)
by_calc = zeros(1, length(y_calc));
for nn = 1 : length(bn)
    by_calc = by_calc + bn(nn) * sin(nn * (2 * pi / l_u) * z_calc) * cosh(nn * (2 * pi / l_u) * (y_calc - y_min));
end

%Calculate changes
dby_calc = by_calc - min(by_calc);
dby_calc_by = dby_calc / min(by_calc);

%Plot the calculated changes
figure;
hold on;
plot(y, dby_abs_by, 'b');
plot(y, dby_abs_by_nom, 'r--');
plot(y_calc, dby_calc_by, 'r');
xlabel('Y (m)');
ylabel('Change in |By| / |By|');
title(['SN', int2str(ser_num), ', Run ', int2str(run_num), ', Change in |By|']);
legend('Measurements', 'Cosh(k dy) - 1', 'Calculation Using Bn', 'location', 'northwest');
grid on;
hold off;

%Done