%**************************************************************************
%This program reads fluxgate measurements of the Earth's magnetic field.
%It calculates field integrals for corrections to long coil measurements.
%
%Zachary Wolf
%9/11/08

%Clear the workspace
clear

%Get the fluxgate measurements
load earthfield.dat;
point_num = earthfield(:, 1)';
x = earthfield(:, 2)';
y = earthfield(:, 3)';
z = earthfield(:, 4)';
bx = earthfield(:, 5)';
by = earthfield(:, 6)';
bz = earthfield(:, 7)';

%Correct a few data points
z(5) = 1.;
z(10) = 2.;
z(15) = 3.;
z(20) = 4.;

%Add an initial point at the start of the coil that the fluxgate couldn't reach
x = [0, x];
y = [0, y];
z = [0, z];
bx = [bx(1), bx];
by = [by(1), by];
bz = [bz(1), bz];

%Correct for the direction of the fluxgate
bx = -bx;
by = -by;

%Plot the fluxgate measurements
figure;
hold on;
plot(z, bx, 'b*');
plot(z, by, 'r*');
plot(z, bz, 'k*');
title('Building 6 Earth Field Measurements');
xlabel('Z Position (m)');
ylabel('Bx, By, Bz (T)');
legend('Bx', 'By', 'Bz');
grid on;
zoom on;
hold off;

%Calculate the first and second field integrals
weight = (4.8 - z);
bx_weight = bx .* weight;
by_weight = by .* weight;
I1x = trapz(z, bx);
I2x = trapz(z, bx_weight);
I1y = trapz(z, by);
I2y = trapz(z, by_weight);

%Display the field integrals
I1x
I2x
I1y
I2y

%Calculate the first and second field integrals outside the undulator
COIL_LENGTH = 4.7498;  %(m)
UND_START_POS = 1.51384;  %(m)
UND_END_POS = COIL_LENGTH - 1.54686;  %(m)
index_array = find((z > UND_START_POS) & (z < UND_END_POS));
bx_outside = bx;
bx_outside(index_array) = 0;
by_outside = by;
by_outside(index_array) = 0;
weight = (4.8 - z);
bx_outside_weight = bx_outside .* weight;
by_outside_weight = by_outside .* weight;
I1x_outside = trapz(z, bx_outside);
I2x_outside = trapz(z, bx_outside_weight);
I1y_outside = trapz(z, by_outside);
I2y_outside = trapz(z, by_outside_weight);

%Display the field integrals
I1x_outside
I2x_outside
I1y_outside
I2y_outside

%Save all results
save 'earthfield.mat';

%Done