%TAYLORDATAFIT: MATLAB script M-file that analyzes the data stored
%in taylordata.m
%
%define and fit the data
taylordata;
%
%VERIFICATION OF THE 2/5 POWER LAW
plot(log(t),log(R),'o')
title('Data Plot: ln R Against ln t','Interpreter','latex','FontSize',16)
xlabel('$\ln t$','Interpreter','latex','FontSize',14)
ylabel('$\ln R$','Interpreter','latex','FontSize',14)
pause
p=polyfit(log(t),log(R),1);
plot(log(t),log(R),'o',log(t),p(1)*log(t)+p(2))
title('Data and Fit: ln R Against ln t','Interpreter','latex','FontSize',16)
xlabel('$\ln t$','Interpreter','latex','FontSize',14)
ylabel('$\ln R$','Interpreter','latex','FontSize',14)
legend('Data','Fit','location','Northwest','Interpreter','latex')
fprintf('Slope = %5.4f; Should be .4000.\n',p(1));
pause
%
%METHOD 1 FOR COMPUTING C: DIRECT FIT
c = t'.^(2/5)\R' %i.e., F = t'.^(2/5)
%
%Compute E
E = (c*rho^(1/5)/phi0)^5 %in Joules (kg m^2/s^2)
pause
%One metric ton (1,000 kg or 2,204.62 lbs) of TNT corresponds, by definition, to 
%an energy of 4.184e9 Joules, but the definition was different in 1950,
%and Taylor used 4.25e9.
Etons = E/4.25e9;
fprintf('This corresponds with %7.2f metric tons of TNT\n',Etons);
pause
%
%METHOD 2 FOR COMPUTING C: TAYLOR'S APPROACH
%NOTE: Taylor used cm, and log10
Y = (5/2)*log10(100*R)-log10(t);
F = ones(size(R))';
ctemp = F\Y'
fprintf('Taylor reported 11.915; see p. 176. Probably round-off error (on his part).\n')
pause
our_taylorc = 10^((2/5)*ctemp)
our_taylorc_our_units = our_taylorc/100  %Return from cm to m.
ourEtaylor = (our_taylorc_our_units*rho^(1/5)/phi0)^5 %in Joules (kg m^2/s^2)
pause
ourEtaylortons = ourEtaylor/4.25e9;
fprintf('This corresponds with %7.2f metric tons of TNT\n',ourEtaylortons);
pause
fprintf('Here is the calculation Taylor actually did.\n')
taylorc = 10^((2/5)*11.915)
fprintf('This corresponds with c^5 = %7.2e \n',taylorc^5)
fprintf('The value he computed with was c^5 = 6.67 x 10^23 cm^5 s^-2. More round-off error? \n')
pause
%There seems to be a typesetting error on p. 177 of Taylor's second paper.
%The value of c^5 is reported as c^5 = 6.67 x 10^2, suggesting the 3 was
%dropped off by the typesetter. In our units, c^5 = 6.67 x 10^13.
Etaylor = (6.67e13*rho/phi0^5)/4.25e9
fprintf('This corresponds with %7.2f metric tons of TNT\n',Etaylor);
fprintf('Taylor reported 16,800 metric tons. In most references the energy is\n')
fprintf('now reported as 20,000 metric tons.\n')


