%**************************************************************************
%This program plots magnetic field measurements vs time.
%Zachary Wolf
%1/6/17

%Initialize
clear;

%Parameters
root_dir = 'C:\Magdata\Dieter Mag Shield\HH Calib';

%Make a list of data files
file_list = dir(root_dir);
ii = 0;
for jj = 1 : length(file_list)
    if strfind(file_list(jj).name, '.dat')
        ii = ii + 1;
        data_file_name{ii} = strcat(root_dir, '\', file_list(jj).name);
    end
end

%Loop over data files -----
for nn = 1 : length(data_file_name)

%Names for input and output files
input_file_name = data_file_name{nn};
plt_file = strrep(data_file_name{nn}, '.dat', '.pdf');

%Initialize the plot file
util_plot_init(plt_file);

%Open the data file
fid = fopen(input_file_name, 'r');
if fid == -1
    disp('magnetic_shield: unable to open input file');
end

%Look for the column headers
while 1
    in_line = fgets(fid);
    if strfind(in_line, '-----') > 0
        break;
    end
end

%Read the data
in_line = 'init';
ii = 0;
while 1
    in_line = fgets(fid);
    if ischar(in_line) ~= 1
        break;
    end
    ii = ii + 1;
    data_array = sscanf(in_line, '%f %f %f %f');
    t1(ii) = data_array(1);
    v1(ii) = data_array(2);  %Current
    t2(ii) = data_array(3);
    v2(ii) = data_array(4);  %Field
end

%Calculate the current and the field
cur = v1 * 40 / 10;  %10 V = 40 A
field = v2 * 10. / 10.;  %10 V = 10 G

%Close the file
fclose(fid);

%Fit the current vs time
f = fit(t1', cur', 'sin1');
amp_i = f.a1;
omega_i = f.b1;
phase_i = f.c1;

%Fit the field vs time
f = fit(t1', field', 'sin1');
amp_b = f.a1;
omega_b = f.b1;
phase_b = f.c1;

%Save the fit results
omega_sav(nn) = omega_i;
cur_sav(nn) = amp_i;
field_sav(nn) = amp_b;

%Make a message for the current plot
msg_i{1} = 'Fit i = i0 sin(omega t + ph)';
msg_i{2} = sprintf('  i0 = %6.3f A', amp_i);
msg_i{3} = sprintf('  omega = %7.3f 1/s', omega_i);
msg_i{4} = sprintf('  ph = %8.2f', phase_i);

%Make a message for the field plot
msg_b{1} = 'Fit b = b0 sin(omega t + ph)';
msg_b{2} = sprintf('  b0 = %6.3f G', amp_b);
msg_b{3} = sprintf('  omega = %7.3f 1/s', omega_b);
msg_b{4} = sprintf('  ph = %8.2f', phase_b);

%Plot the field outside the shield
figure;
subplot(2, 1, 1);
hold all;
plot(t1, cur, 'b-');
title('Helmholtz Coil Current vs Time');
xlabel('Time (sec)');
ylabel('Current (A)');
for jj = 1 : length(msg_i)
    text(min(xlim) + .1 * (max(xlim) - min(xlim)), max(ylim) - (jj + 1) * .1 * (max(ylim) - min(ylim)), msg_i{jj});
end
grid on;
zoom on;
hold off;

%Plot the field inside the shield
subplot(2, 1, 2);
hold all;
plot(t2, field, 'b-');
title('Helmholtz Coil Field vs Time');
xlabel('Time (sec)');
ylabel('Magnetic Field (Gauss)');
for jj = 1 : length(msg_b)
    text(min(xlim) + .1 * (max(xlim) - min(xlim)), max(ylim) - (jj + 1) * .1 * (max(ylim) - min(ylim)), msg_b{jj});
end
grid on;
zoom on;
hold off;
util_print_fig(plt_file);

%End of loop over data files -----
end

%Make summary plots
summ_plt_file = strcat(root_dir, '\', 'calib_summ.pdf');
util_plot_init(summ_plt_file);
%
figure;
hold all;
plot(omega_sav / (2*pi), field_sav ./ curr_sav, 'b*');
title('Helmholtz Coil Transfer Function');
xlabel('Frequency (1/sec)');
ylabel('Field / Current (G / A)');
grid on;
zoom on;
hold off;
%
figure;
hold all;
plot(omega_sav / (2*pi), field_sav ./ curr_sav, 'b*');
title('Helmholtz Coil Transfer Function');
xlabel('Frequency (1/sec)');
ylabel('Field / Current (G / A)');
grid on;
zoom on;
hold off;

%Done