% Read raw data and construct Allan variance from rawstats data file
% Data are collected using rawstats filtered to remove all peers but usno_raw.
%
clear all;
fid = fopen('usno_raw', 'r');
[table, count] = fscanf(fid, '%f %f %f %f %f %f');
a = table(1:6:count); b = table(2:6:count);
t = table(1:6:count) * 86400;
t = t + table(2:6:count) - t(1);
t1 = table(3:6:count);
t2 = table(4:6:count);
t3 = table(5:6:count);
t4 = table(6:6:count);
offset = ((t2 - t1) - (t4 - t3)) / 2;
delay = (t4 - t1) - (t3 - t2);
%
% Wedge
%
clf reset; h = newplot;
set(h, 'FontSize', 12);
set(h, 'Position', [.12 .15 .85 .8]);
plot(delay * 1e3, offset * 1e3, '.', '-k')
xlabel('Delay (ms)');
ylabel('Offset (ms)');
print -dtiff -r300 usno_wedge
%
% Delay
%
clf reset; h = newplot;
set(h, 'FontSize', 12);
set(h, 'Position', [.12 .15 .85 .8]);
plot(t / 86400, delay * 1e3, '-k')
xlabel('Time (day)');
ylabel('Delay (ms)');
print -dtiff -r300 usno_delay
%
% Variance-time
%
x1 = t / 86400;
y1 = delay; i = 1; d = 16;
while length(y1) >= 10    
	z1(i) = var(y1);
	m1(i) = d;
    clf reset; h = newplot;
    set(h, 'FontSize', 12);
    set(h, 'Position', [.12 .15 .85 .8]);
    plot(x1, y1 * 1e3, '-k')
    axis([0 12 0 500]);
    lag = strcat('Time (ms) Interval=', num2str(d));
    xlabel(lag);
    ylabel('Delay (ms)');
    lag = strcat('usno_', num2str(d));
    print('-dtiff', '-r600', lag)
    x1 = (x1(2:2:length(x1)));
    y1 = (y1(2:2:length(y1)) + y1(1:2:length(y1) - 1)) / 2;
    i = i + 1; d = 2 * d;
end
e = m1 / m1(1);
e = 1 ./ e;
loglog(m1, z1, '-k', m1, e * z1(1))
%axis([1 1e5 1e-4 100]);
xlabel('Time Interval (s)');
ylabel('Variance (s^2)');
print -dtiff -r600 usno_var
return;
%
% fft
%
f = fft(expon);
r = f .* conj(f);
loglog(r(length(r) / 2:length(r)), '.')
%axis([1 1e5 1e-4 100]);
xlabel('Frequency (Hz)');
ylabel('Power');
print -dtiff -r600 usno_per
%
% fft
%
f = fft(delay);
r = f .* conj(f);
loglog(r(length(r) / 2:length(r)), '-k')
%axis([1 1e5 1e-4 100]);
xlabel('Frequency (Hz)');
ylabel('Power');
print -dtiff -r600 usno_per