/* "Experiments to Determine the force of Fired Gunpowder"
by Benjamin Count of Rumford. Philosophical Transactions of the
Royal Society of London. Vol. 87 (1797) pp222-292.

Theory: X=charge of gunpowder   Y=elasticity of generated fluid

  Y = X**(1+0.0004X)  = (B0*X)**(B1+B2*X)
  B0=1=B1  B2=0.0004 (Count's Theory)

*/

options nodate nonumber ps=54 ls=76;
goptions reset=all dev=pdf gsfname=output gsfmode=replace;

data gunpwder;
input charge elasticity;
elasticity=elasticity/1.841;    /* Converts atmospheres back to =parts  */
lnchar=log(charge);
cards;
39 77.86
78 182.30
117 228.2
156 382.4
195 561.2
234 685.6
273 811.7
312 1164.8
351 1551.3
390 1884.3
429 2219.0
468 2573.7
507 3283.3
546 4008.0
585 4722.5
624 7090.0
702 10977.0
;
run;

proc nlin;
parms b0=1 b1=1 b2=0.0004;
model elasticity = b0*(charge**(b1+(b2*charge)));
der.b0=(charge**(b1+(b2*charge)));
der.b1=lnchar*b0*(charge**(b1+(b2*charge)));
der.b2=charge*lnchar*b0*(charge**(b1+(b2*charge)));
output out=gunout p=yhat r=e;
run;

data gunout; set gunout;
theory=charge**(1+(0.0004*charge));

proc sort;
by charge;
run;

filename output 'gunpowder.pdf';
symbol1 c=black v=:;
symbol2 c=black l=1 i=spline;
symbol3 c=black l=5 i=spline;

proc gplot;
plot elasticity*charge=1 yhat*charge=2 theory*charge=3 / overlay frame;
title 'Gunpowder Experiment';
run;

filename output 'gunpowder1.pdf';
proc gplot;
plot e*yhat=1 / frame vref=0;
run;

quit;