% Nominal values of input parameters

hgas  = 3000;    % TBC-gas heat transfer coef. (W/(m^2 K))
Tgas  = 1500;    % Mixed gas temperature (K)
ktbc  = 1;       % TBC   thermal conduct. (W/mK)
km    = 20;      % Metal thermal conduct. (W/mK)
Lm    = 0.003;   % Metal thickness (m)
hcool = 1000;    % Coolant-metal heat transfer coef. (W/(m^2 K))
Tcool =  600;    % Coolant temperature (K)

% Ltbc min, mpp, and max values
Ltbc_min = 0.00025;
Ltbc_mpp = 0.00030;
Ltbc_max = 0.00075;

% Set termination criteria
Nmin     = 25;
eps_mean = 0.1;

% Zero out mean and standard deviation
sumTmh  = 0;
sumTmh2 = 0;

n = 0;
MetError = false;

while ( (n < Nmin) | (~MetError) ),
  
  % generate Ltbc values using a triangular distribution
  Ltbc = trirnd(Ltbc_min,Ltbc_mpp,Ltbc_max);

  % Solve heat transfer problem

  [Ttbc, Tmh, Tmc, q] = blade1D(hgas, Tgas, ...
				   ktbc, Ltbc, ...
				   km, Lm, ...
				   hcool, Tcool);

  n = n + 1;
  sumTmh  = sumTmh  + Tmh;
  sumTmh2 = sumTmh2 + Tmh^2;
  
  meanTmh  = sumTmh/n;
  meanTmh2 = sumTmh2/n;
  if (n>1),
    sdevTmh = sqrt( n/(n-1)*(meanTmh2 - meanTmh^2) );
  else
    sdevTmh = 0;
  end
  
  Tmh_conf = 3*sdevTmh/sqrt(n);
  MetError = (Tmh_conf < eps_mean);

end
fprintf('%i: Mean, Sdev, Conf Tmh = %f, %f, %f\n',n,meanTmh,sdevTmh,Tmh_conf);