%  Energy Model: Curve Fitting
%  ---------------------------
%
%  World oil production data measured in millions of barrels/day
%
t = 1900:2006;
po = 0:0.1:0.9;
po = [ po , 1:0.1:1.9 ];
po = [ po , 2:0.3:4.7 ];
po = [ po , 5 , 5.1, 5.2, 5.2, 5.3, 5.3, 5.7, 6, 6, 6 ];
po = [ po , 6 , 6, 6, 6.4, 6.9, 7.4, 8, 8.7, 9.5, 10.4 ];
po = [ po , 11.2, 12.3, 13.3, 13.8, 14.1, 16.7, 18.8, 19.2, 20, 20.8 ];
po = [ po, 21.9, 23.3, 26.6, 27.9, 29.6, 31.803,34.568,37.118,40.436,43.633 ];
po = [ po , 48.061,50.844,53.666,58.464,58.617,55.825,60.411,62.712,63.330,66.048 ];
po = [ po , 62.946,59.532,57.296,56.595,57.683,57.469,60.463,60.784,63.157,64.049 ];
po = [ po , 65.470,65.287,65.795,66.051,67.122,68.125,69.931,72.251,73.626,72.439 ];
po = [ po , 75.033,74.932,74.496,77.056,80.244,81.250,81.663 ];
%
%  Convert data to barrels per year
%
po = 365*1e6*po;
%
%  Compute total produced oil measured in barrels
%
pp = 0;
for i=2:107,
  pp = [ pp , pp(i-1) + po(i-1) ];
end
%
%  Identify parameters of logistic model
%
x = pp(86:107)';
VDM = [ x , x.*x ];
coef = VDM \ po(86:107)';
%
%  Calculate maximum value
%
a = coef(1);
b = coef(2);
ppmax = -a/b
%
%  Simulate the logistic model
%
t2 = t(30:107);
po2 = po(30:107);
t3 = 1985:2100;
pp3 = pp(86);
po3 = zeros(1,116);
po3(1) = po(86);
for i = 2:116,
  po3(i) = a*pp3(i-1) + b*pp3(i-1)*pp3(i-1);
  pp3 = [ pp3 , pp3(i-1) + po3(i-1) ];
end
%
%  Simulate constant oil production until the reserves are used up
%
t4 = 2007:2100;
pp4 = pp(107);
po4 = zeros(1,94);
po4(1) = po(107);
for i=2:94,
  if pp4(i-1) < ppmax,
    po4(i) = po(107);
   else
    po4(i) = 0;
  end;
  pp4 = [ pp4 , pp4(i-1) + po4(i-1) ];
end
%
%  Simulate exponential oil production based on the average exponential growth
%  measured during the last 10 years, until the reserves are used up
%
pd = 0;
for i=98:107,
  pd = pd + po(i)/po(i-1);
end
pd = pd/10
t5 = t4;
pp5 = pp(107);
po5 = zeros(1,94);
po5(1) = po(107);
for i=2:94,
  if pp5(i-1) < ppmax,
    po5(i) = pd*po5(i-1);
   else
    po5(i) = 0;
  end;
  pp5 = [ pp5 , pp5(i-1) + po5(i-1) ];
end
%
%  Plot data
%
%
subplot(2,1,1)
plot(t2,po2,'k-',t3,po3,'b-',t4,po4,'r-',t5,po5,'g-')
grid on
title('World Oil Production')
xlabel('time (years)')
ylabel('produced oil (barrels/year)')
text(1961,1.5e10,'(a)')
text(2025,2.2e10,'(b)')
text(2054,2.5e10,'(c)')
text(2043,4.5e10,'(d)')
print -djpeg TOD_08a.jpg
pause
%
%  Let us now curfe-fit a logistic population model
%
tt = 1950:10:2050;
p = [ 2555360972 , 3039585530 , 3707475887 , 4452584592 , 5281653820 , 6079603571 ];
p = [ p , 6815892190 , 7541773753 , 8175075482 , 8722646253 , 9190252532 ];
pdot = zeros(1,9);
for i=1:9,
  pdot(i) = ( p(i+2) - p(i) ) / 20;
end
px = p(2:10)';
VDM = [ px , px.*px ];
coef = VDM \ pdot';
a = coef(1);
b = coef(2);
%
%  Plot data
%
t6 = 1950:2100;
p6 = p(1);
for i=2:151,
  pdt = a*p6(i-1) + b*p6(i-1)*p6(i-1);
  p6 = [ p6 , p6(i-1) + pdt ];
end
clf
subplot(2,1,1)
plot(t6,p6,'k-')
grid on
title('World Population')
xlabel('time (years)')
ylabel('population')
print -djpeg TOD_08b.jpg
pause
%
%  Now compute pere capita oil production
%
conv = 42*3.6/365;
t2x = t2(21:78);
po2pc = conv*po2(21:78)./p6(1:58);
t3x = t3;
po3pc = conv*po3./p6(36:151);
t4x = t4;
po4pc = conv*po4./p6(58:151);
t5x = t5;
po5pc = conv*po5./p6(58:151);
clf
subplot(2,1,1)
plot(t2x,po2pc,'k-',t3x,po3pc,'b-',t4x,po4pc,'r-',t5x,po5pc,'g-')
grid on
title('World Per Capita Oil Production')
xlabel('time (years)')
ylabel('oil (liters/person/day)')
text(1963,1.7,'(a)')
text(2023,1.2,'(b)')
text(2054,0.8,'(c)')
text(2043,1.8,'(d)')
print -djpeg TOD_08c.jpg
return


