%**************************************************************************
%mag_har_study
%This program reads in the Br vs s values from a single block and
%determines the harmonic content of the field.
%
%Zachary Wolf
%1/11/16

%Clear all variables
clear;

%Define the undulator parameters
lambda_u = .042;  %(m)
ku = 2 * pi / lambda_u;

%Input data file
input_dir = 'C:\Magdata\HD Delta Block Studies\d520\';
input_file = 'quadrantField_Br_vs_s_520-.txt';

%Make file names
plt_file = strcat(input_dir, 'mag_har_study.ps');
mat_file = strcat(input_dir, 'mag_har_study.mat');

%Initialize the plot file
util_plot_init(plt_file);

%Get the field values
file_name = strcat(input_dir, input_file);
fid = fopen(file_name, 'r');
while 1
    in_line = fgets(fid);
    if findstr('-----', in_line) > 0
        break;
    end
end
data_array = fscanf(fid, '%f %f', [2, inf]);
x = data_array(1, :)';
x = x / 1000;  %(m)
by = data_array(2, :)';
fclose(fid);

%Do a Fourier series expansion of the field
index_array = find(abs(x) < .007);
x_fit = x(index_array);
by_fit = by(index_array);
f = fit(x_fit, by_fit, 'sin8');
a1 = f.a1;
b1 = f.b1;
c1 = f.c1;
a2 = f.a2;
b2 = f.b2;
c2 = f.c2;
a3 = f.a3;
b3 = f.b3;
c3 = f.c3;
a4 = f.a4;
b4 = f.b4;
c4 = f.c4;
a5 = f.a5;
b5 = f.b5;
c5 = f.c5;
a6 = f.a6;
b6 = f.b6;
c6 = f.c6;
a7 = f.a7;
b7 = f.b7;
c7 = f.c7;
a8 = f.a8;
b8 = f.b8;
c8 = f.c8;
ks = b1;

%Put coefficients into arrays
amp_vals = [a1 * sin(c1), a2 * sin(c2), a3 * sin(c3), a4 * sin(c4), a5 * sin(c5), a6 * sin(c6), a7 * sin(c7), a8 * sin(c8)];
ks_vals = [b1, b2, b3, b4, b5, b6, b7, b8];
kr_vals = sqrt(ku^2 + ks_vals.^2);
phi0_vals = amp_vals ./ kr_vals;
grad_k_vals = phi0_vals .* ks_vals.^2;

%Metric to cancell gradient
metric = sum(grad_k_vals);

%Check the fit
xx = x_fit;
yy = 0;
for ii = 1 : length(amp_vals)
    yy = yy + amp_vals(ii) * cos(ks_vals(ii) * xx);
end
figure;
hold all;
plot(x, by, 'b-');
plot(xx, yy, 'r-');
legend('By', 'Fit By');
title('By vs X');
xlabel('X (m)');
ylabel('By (T)');
grid on;
zoom on;
hold off;
orient landscape;
util_print_fig(plt_file);

%Plot the Fourier series
ii = 1 : length(amp_vals);
figure;
hold all;
plot(ii, amp_vals, 'b*');
title(['FT By vs Harmonic Number, ks = ', num2str(ks), ' 1/m', ', Metric = ', num2str(metric), ' T/m']);
xlabel('N');
ylabel('By(n) (T)');
grid on;
zoom on;
hold off;
orient landscape;
util_print_fig(plt_file);

%Plot the phi0 values
ii = 1 : length(phi0_vals);
figure;
hold all;
plot(ii, phi0_vals, 'b*');
title(['Phi0 vs Harmonic Number, ks = ', num2str(ks), ' 1/m', ', Metric = ', num2str(metric), ' T/m']);
xlabel('N');
ylabel('Phi0(n) (Tm)');
grid on;
zoom on;
hold off;
orient landscape;
util_print_fig(plt_file);

%Plot the grad_k values
ii = 1 : length(grad_k_vals);
figure;
hold all;
plot(ii, grad_k_vals, 'b*');
title(['Grad K Values vs Harmonic Number, ks = ', num2str(ks), ' 1/m', ', Metric = ', num2str(metric), ' T/m']);
xlabel('N');
ylabel('Grad K Val(n) (T/m)');
grid on;
zoom on;
hold off;
orient landscape;
util_print_fig(plt_file);

%Plot the metric
ii = 1 : length(amp_vals);
metric_vals = cumsum(grad_k_vals);
figure;
hold all;
plot(ii, metric_vals, 'b*');
title(['Sum(amp * ks**2 / kr) up to Harmonic Number, ks = ', num2str(ks), ' 1/m', ', Metric = ', num2str(metric), ' T/m']);
xlabel('N');
ylabel('Metric (T/m)');
grid on;
zoom on;
hold off;
orient landscape;
util_print_fig(plt_file);

%Output the results to the output file directory
save(mat_file, 'f');

%Done
disp('Done');
