kblade_vals = [ 28.5,    30,  31.5 ];
Tgas_vals   = [ 1400,  1500,  1600 ];
hgasLE_vals = [15500, 16000, 18000 ];
hgasTE_vals = [ 3750,  4000,  4250 ];
Tcool_vals  = [  550,   600,   650 ];
hcool_vals  = [ 1400,  1500,  1600 ];

Tdes = 1430;
Tlim = 1500;

Tgraddes = 70000;
Tgradlim = 80000;

chord = 0.04;
[tri2nod, xy, bedge] = loadblade('hpblade_coarse', chord);

conf_lev = 3;
err_lev  = 0.01;

Nmin = 300;
Nmax = 15000;
Tmax     = zeros(Nmax,1);
Tgradmax = zeros(Nmax,1);
D        = zeros(Nmax,1);
params   = zeros(Nmax,6);

MetError = false;
Nsize    = 0;
Psum     = 0;

while ( ((Nsize < Nmin) | (~MetError)) & (Nsize < Nmax) ),

  Nsize = Nsize + 1;
  
  params(Nsize,1) = trirnd(kblade_vals(1), kblade_vals(2), kblade_vals(3));
  params(Nsize,2) = trirnd(Tgas_vals(1), Tgas_vals(2), Tgas_vals(3));
  params(Nsize,3) = trirnd(hgasLE_vals(1), hgasLE_vals(2), hgasLE_vals(3));
  params(Nsize,4) = trirnd(hgasTE_vals(1), hgasTE_vals(2), hgasTE_vals(3));
  params(Nsize,5) = trirnd(Tcool_vals(1), Tcool_vals(2), Tcool_vals(3));
  params(Nsize,6) = trirnd(hcool_vals(1), hcool_vals(2), hcool_vals(3));
  
  [T, Tgrad] = calcblade( params(Nsize,:), chord, tri2nod, xy, bedge);
  
  Tmax(Nsize)     = max(T);
  Tgradmax(Nsize) = max(Tgrad);
  
  D(Nsize) = (Tmax(Nsize)-Tdes)/(Tlim-Tdes) + ...
         (Tgradmax(Nsize)-Tgraddes)/(Tgradlim-Tgraddes);
  
  Psum   = Psum + (D(Nsize)>1);
  Pfail  = Psum/Nsize;
  sPfail = sqrt(Pfail*(1-Pfail)/Nsize);
  
  MetError = (conf_lev*sPfail < err_lev);
  
  if (rem(Nsize,100) == 0),
    fprintf('Nsize = %i: P = %f +/- %f at 99%% confidence\n',Nsize,Pfail,conf_lev*sPfail);
  end
  
end

fprintf('\n\n Monte Carlo simulation terminated.\n\n');
fprintf('Nsize = %i: P = %f +/- %f at 99%% confidence\n',Nsize,Pfail,conf_lev*sPfail);

% Plot histograms and CDf's of Tmax, Tgradmax, and D
figure(1);
subplot(211);
hist(Tmax(1:Nsize),40);
xlabel('T_{max}');
subplot(212);
cdfplot(Tmax(1:Nsize));
xlabel('T_{max}');
hold on; plot([Tlim,Tlim],[0,1],'r'); hold off;

figure(2);
subplot(211);
hist(Tgradmax(1:Nsize),40);
xlabel('|\nabla T|_{max}');
subplot(212);
cdfplot(Tgradmax(1:Nsize));
xlabel('|\nabla T|_{max}');
hold on; plot([Tgradlim,Tgradlim],[0,1],'r'); hold off;

figure(3);
subplot(211);
hist(D(1:Nsize),40);
xlabel('Damage');
subplot(212);
cdfplot(D(1:Nsize));
xlabel('Damage');
hold on; plot([1,1],[0,1],'r'); hold off;