
% % form WHO2004, butter outbreak data
    p=1.0e-3;% infection probability in the outbreak
    dose=8.2e3;
    r_d=-(1/dose)*log(1-p);% r for epi curve
%------------------------
dp=25;% number of doses
y_exact=nan(dp,1);
u=10.^linspace(1,12,dp);% doses
for j=1:dp
    IP=u(j);
    y_exact(j)=1-exp(-r_d*IP);
end

%figure(3)
  hold on;
   plot(log10(u(1:end)),(1-exp(-r_d*u(1:end))),'go--','LineWidth',2)
   % min and max responses (probabilities)
r_mean=[1.69E-06	1.64E-06	1.70E-06	1.72E-06	3.00E-06	3.16E-05	5.96E-04	6.62E-03	4.32E-02	1.72E-01	4.17E-01	7.00E-01	8.85E-01	9.66E-01	9.90E-01	9.94E-01	9.96E-01	9.97E-01	9.97E-01	9.97E-01	9.97E-01	9.97E-01	9.96E-01	9.97E-01	9.97E-01];


r_Max=[1.64E-06	1.61E-06	1.68E-06	1.68E-06	4.32E-06	3.38E-05	6.31E-04	6.36E-03	4.14E-02	1.63E-01	3.93E-01	6.58E-01	8.27E-01	9.06E-01	9.26E-01	9.30E-01	9.31E-01	9.31E-01	9.32E-01	9.31E-01	9.31E-01	9.32E-01	9.31E-01	9.31E-01	9.31E-01];%95% upper limit
r_Min=[1.53E-06	1.49E-06	1.50E-06	1.54E-06	2.08E-06	2.63E-05	5.08E-04	6.04E-03	3.92E-02	1.59E-01	3.87E-01	6.51E-01	8.22E-01	8.99E-01	9.23E-01	9.26E-01	9.29E-01	9.30E-01	9.30E-01	9.30E-01	9.29E-01	9.30E-01	9.29E-01	9.29E-01	9.30E-01];%95% lower limit

   hold on;
%----------shading-----------------------
x=1:.5:13;
length(x)
 plot(x,r_mean,'m-','LineWidth',2)
 xlabel('Inoculum (Log dose)','fontsize',16);ylabel('Probability of infection','fontsize',16)
x2 = [x, fliplr(x)];
inBetween = [r_Min, fliplr(r_Max)];
fill(x2, inBetween, 'g'); % shaded region
