
%Load simulation data
sim_data = load('C:\Magdata\Phase Shifter Err vs Gap Sim\Baseline\Gap10mm\phase_anal.mat');
z_sim = sim_data.field.z;
bx_sim = sim_data.field.bx;
by_sim = sim_data.field.by;

%Load measured data
meas_data = load('C:\Magdata\Phase Shifter Err vs Gap Sim\Block3 dy200um vs Gap\Gap10mm\phase_anal.mat');
z_meas = meas_data.field.z;
bx_meas = meas_data.field.bx;
by_meas = meas_data.field.by;

%Get the z positions of the field peaks
z_peak_meas = meas_data.field.z_peak;
by_peak_meas = meas_data.field.by_peak;
z_peak_sim = sim_data.field.z_peak;
by_peak_sim = sim_data.field.by_peak;

%Check the peaks
if length(z_peak_meas) ~= length(z_peak_sim)
    warndlg('The measurements and simulation have different number of peaks');
end
z_peak_diff = z_peak_meas - z_peak_sim;
if abs(z_peak_diff) > .002
    warndlg('The measurements and simulation have different peak positiions');
end

%Extract simulated data over the same z range as the measurements
index_array = find((z_sim > z_meas(1)) & (z_sim < z_meas(end)));
z_s = z_sim(index_array);
bx_s = bx_sim(index_array);
by_s = by_sim(index_array);

%Fit the simulated data to find its value at the measured positions
bx_sim_fit = spline(z_s, bx_s, z_meas);
by_sim_fit = spline(z_s, by_s, z_meas);

%Use data with no fit
% bx_sim_fit = bx_sim;
% by_sim_fit = by_sim;

%Use the peak fields to scale the simulation to match the measurements
for ii = 1 : length(by_peak_meas)
    scale_fact(ii) = abs(by_peak_meas(ii)) / abs(by_peak_sim(ii));
    dis(ii) = sum(abs(by_meas - scale_fact(ii) * by_sim_fit));
end
[dis_min, index_min] = min(dis);
by_sim_fit_scaled = scale_fact(index_min) * by_sim_fit;
by_peak_sim_scaled = scale_fact(index_min) * by_peak_sim;

%Scale the By simulation to match the measurements, minimize variance
% scale_fact = sum(by_meas .* by_sim_fit) / sum(by_sim_fit.^2);
% by_sim_fit_scaled = scale_fact * by_sim_fit;

%Plot the measurements and the simulation
figure;
hold on;
plot(z_sim, bx_sim, '-r', 'linewidth', 0.5);
plot(z_meas, bx_meas, '-b', 'linewidth', 0.5);
legend('Sim', 'Meas');
title('Bx vs Z');
xlabel('Z (m)');
ylabel('Bx (T)');
grid on;
zoom on;
hold off;

figure;
hold on;
plot(z_sim, by_sim, '-r', 'linewidth', 0.5);
plot(z_meas, by_meas, '-b', 'linewidth', 0.5);
plot(z_peak_meas, by_peak_meas, 'b*');
legend('Sim Scaled', 'Meas');
title('By vs Z');
xlabel('Z (m)');
ylabel('By (T)');
grid on;
zoom on;
hold off;


%Plot the difference between the measurements and the simulation
figure;
hold on;
plot(z_meas, bx_meas - bx_sim_fit, 'linewidth', 0.5);
title('Bx Meas - Bx Sim vs Z');
xlabel('Z (m)');
ylabel('Bx Meas - Bx Sim (T)');
for ii = 1 : length(z_peak_meas)
    line([z_peak_meas(ii), z_peak_meas(ii)], ylim, 'linestyle', '--');
end
xlim([z_peak_meas(1) - 0.05, z_peak_meas(end) + 0.05]);
legend('Bx Difference', 'Field Peaks');
grid on;
zoom on;
hold off;

figure;
hold on;
plot(z_meas, by_meas - by_sim_fit_scaled, 'linewidth', 0.5);
title('By Meas - By Sim Scaled vs Z');
%title('By Err - By NoErr vs Z');
xlabel('Z (m)');
%ylabel('By Err - By NoErr (T)');
ylabel('By Meas - By Sim (T)');
for ii = 1 : length(z_peak_meas)
    line([z_peak_meas(ii), z_peak_meas(ii)], ylim, 'linestyle', '--');
end
xlim([z_peak_meas(1) - 0.05, z_peak_meas(end) + 0.05]);
legend('By Difference', 'Field Peaks');
grid on;
zoom on;
hold off;

figure;
hold on;
plot(z_peak_meas, by_peak_meas - by_peak_sim_scaled, '*');
title('By Peak Meas - By Peak Sim Scaled vs Z');
%title('By Err - By NoErr vs Z');
xlabel('Z (m)');
%ylabel('By Err - By NoErr (T)');
ylabel('By Meas - By Sim (T)');
for ii = 1 : length(z_peak_meas)
    line([z_peak_meas(ii), z_peak_meas(ii)], ylim, 'linestyle', '--');
end
xlim([z_peak_meas(1) - 0.05, z_peak_meas(end) + 0.05]);
legend('By Peak Difference', 'Field Peaks');
grid on;
zoom on;
hold off;
