clear all;

% Lower and upper bounds of input parameters

hgas  = [1500, 4500];    % TBC-gas heat transfer coef. (W/m^2)
Tgas  = [1200, 1400];    % Mixed gas temperature (C)
ktbc  = [0.9, 1.1];       % TBC   thermal conduct. (W/mK)
Ltbc  = [0.00025, 0.00075];  % TBC   thickness (m)
km    = [20.0, 23.0];    % Metal thermal conduct. (W/mK)
Lm    = [0.002, 0.004];  % Metal thickness (m)
hcool = [ 500, 1500];    % Coolant-metal heat transfer coef. (W/m^2)
Tcool = [ 150, 250];    % Coolant temperature (C)

% Place parameter bounds into a vector

Pbound = [hgas; ...
	  Tgas; ...
	  ktbc; ...
	  Ltbc; ...
	  km; ...
	  Lm; ...
	  hcool; ...
	  Tcool];

% Number of Monte Carlo trials
Ntrial = 1000;

for n = 1:Ntrial,
  
  % generate Ltbc values using a uniform distribution

  P(:,n) = Pbound(:,1) + (Pbound(:,2)-Pbound(:,1)).*rand(8,1);

  % Solve heat transfer problem

  [Ttbc, Tmh(n), Tmc, q] = blade1D(P(1,n), P(2,n), P(3,n), P(4,n), ...
				   P(5,n), P(6,n), P(7,n), P(8,n));

  if (n>1),
    mTmh(n-1) = mean(Tmh);
    sTmh(n-1) = std(Tmh);
  end

%  plot(Ltbc(n),Tmh(n),'*');hold on;
%  drawnow;
end

fprintf('Mean Tmh = %f\n',mTmh(Ntrial-1));
fprintf('Sdev Tmh = %f\n',sTmh(Ntrial-1));

subplot(311);
hist(Tmh,20);
xlabel('Tmh (C)');
subplot(312);
plot(mTmh);
subplot(313);
plot(sTmh);