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

load vdp_period_datafile2;

% the first order term
firstord = (3-2*log(2))*muset;


%alpha = 2.338107410459767;  % from maple
% or get it from matlab:
alpha = vdp_first_airy_root;
% second order approx (first ord + 2nd order correction)
secondord = firstord + 3*alpha*muset.^(-1/3);

% estimates of the second order and third order corrections
estsecondord = periods - firstord;
estthirdord = periods - secondord;


figure(4);  clf;  hold on;
plot(muset, firstord, 'b-');
plot(muset, secondord, 'go-');
plot(muset, periods, 'kx-');
xlabel('\mu');
ylabel('period');
legend('analytic 1st order', 'analytic 2nd order', 'computation', 4);
%print('-dpng', '-r40', 'vdp_period_vs_mu_2.png');

figure(5); clf;
plot(muset, estthirdord, 'rx-');
title('Residual third order correction');
xlabel('\mu');
%print('-dpng', '-r40', 'vdp_estimate_3rd_corr.png');

figure(6); clf; hold on;
plot(muset, estsecondord, 'rx-');
plot(muset, secondord-firstord, 'ko-');
title('second order correction');
xlabel('\mu');
legend('estimated second order correction', 'actual second order', 0);
%print('-dpng', '-r40', 'vdp_2nd_corr.png');

logmuset = log(muset);
logest = log(-estthirdord);

% throw away mu values smaller than (remember the expansions are
% valid for large mu)
startmu = 1;
% 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(7);  clf;
%H = loglog(muset, estthirdord, 'k+-');
%set(H, 'Color', .75*[1 1 1]);
%hold on;
loglog(muset(start:end), -estthirdord(start:end), 'ro-');
title(['residual error (slope = '  num2str(C(1))  ')']);
xlabel('log_{10} \mu');
legend(['data used to compute slope (\mu >= '  num2str(startmu)  ')'], 1);
grid on;
%print('-dpng', '-r40', 'vdp_estimate_3rd_corr_power.png');
