% process the data from vdp_period
%
% 2003-11-07 cbm
% no rights reserved: public domain code

load vdp_period_datafile;

% the first order term
firstord = (3-2*log(2))*muset;
% an estimate of what the second order correction should be
estsecondord = periods - firstord;


figure(4);  clf;  hold on;
plot(muset, periods, 'kx-');
plot(muset, firstord, 'bo-');
xlabel('\mu');
ylabel('period');
legend('computation', 'analytic 1st order', 4);
%print('-dpng', '-r50', 'vdp_period_vs_mu.png');

figure(5); clf;
plot(muset, estsecondord, 'rx-');
title('estimated second order correction');
xlabel('\mu');
%print('-dpng', '-r50', 'vdp_estimate_2nd_corr.png');

logmuset = log(muset);
logest = log(estsecondord);

% throw away mu values smaller than (remember the expansions are
% valid for large mu)
startmu = 10;
% trick to compute the index of first value to keep
start = sum(muset < startmu)+1;

disp('polyfit results for all data');
C = polyfit(logmuset, logest, 1)
disp(['polyfit results for all data with \mu >= '  num2str(startmu)]);
C = polyfit(logmuset(start:end), logest(start:end), 1)

FH = figure(6);  clf;
H = loglog(muset, estsecondord, 'k+-');
set(H, 'Color', .75*[1 1 1]);
hold on;
loglog(muset(start:end), estsecondord(start:end), 'ro-');
title(['computed second order correction (slope = '  num2str(C(1))  ')']);
xlabel('\mu');
legend('all data', ['data used to compute slope (\mu >= '  num2str(startmu)  ')'], 1);
grid on;
%print('-dpng', '-r50', 'vdp_estimate_2nd_corr_power.png');