% Clear variables
clear all;

% Set variability on inputs

% Temperature & heat transfer coefficient range on external surfaces
dText = [1200, 1400];
dhext = [2000, 4000];

% Temperature & heat transfer coefficient range on internal surfaces
dTint = [ 150,  250];
dhint = [ 500, 1500];

% Metal conductivity range
dkmet = [20.0, 23.0];

% Blade length
Lblade     = 0.1;

% Load in the grid file
fname = input('Enter gridfile name: ','s');
load(fname);

% Place parameter bounds into a vector
Pbound = [dText; ...
	  dhext; ...
	  dTint; ...
	  dhint; ...
	  dkmet];

% Number of Monte Carlo trials
Ntrial = 1000;

% Zero statistical storage arrays 
meanT    = zeros(Nv,1);
sdevT    = zeros(Nv,1);

Tmin     = zeros(Ntrial,1);
meanTmin = zeros(Ntrial,1);
sdevTmin = zeros(Ntrial,1);

Tmax     = zeros(Ntrial,1);
meanTmax = zeros(Ntrial,1);
sdevTmax = zeros(Ntrial,1);


for n = 1:Ntrial,
  
  % generate parameters using a symmetric triangular distribution
  P = Pbound(:,1) + (Pbound(:,2)-Pbound(:,1)).*symtri(5);

  % Place parameters in correct non-dimensional groupings
  kmetal    = P(5);

  Tcool(1)  = P(1);
  hwall(1)  = P(2)*Lblade/kmetal;
  
  Tcool(2)  = P(3);
  hwall(2)  = P(4)*Lblade/kmetal;

  Tcool(3)  = P(3);
  hwall(3)  = P(4)*Lblade/kmetal;

  Tcool(4)  = P(3);
  hwall(4)  = P(4)*Lblade/kmetal;
  
  Tcool(5)  = P(3);
  hwall(5)  = P(4)*Lblade/kmetal;
  
  % Call finite element solver
  SolveBladeHeat;

  % Store statistical data
  meanT = meanT + Tsol;
  sdevT = sdevT + Tsol.^2;
  Tmin(n)  = min(Tsol);
  Tmax(n)  = max(Tsol);

  meanTmin(n) = mean(Tmin(1:n));
  sdevTmin(n) =  std(Tmin(1:n));
  meanTmax(n) = mean(Tmax(1:n));
  sdevTmax(n) =  std(Tmax(1:n));

end

fprintf('Mean Tmin = %f\n',meanTmin(Ntrial));
fprintf('Sdev Tmin = %f\n',sdevTmin(Ntrial));
fprintf('Mean Tmax = %f\n',meanTmax(Ntrial));
fprintf('Sdev Tmax = %f\n',sdevTmax(Ntrial));

figure;
hist(Tmin,20);
xlabel('Tmin');

figure;
hist(Tmax,20);
xlabel('Tmax');

figure;
plot(meanTmin);
xlabel('Samples');
ylabel('mean Tmin');

figure;
plot(meanTmax);
xlabel('Samples');
ylabel('mean Tmax');

figure;
plot(sdevTmin);
xlabel('Samples');
ylabel('sdev Tmin');

figure;
plot(sdevTmax);
xlabel('Samples');
ylabel('sdev Tmax');

% Calculate mean and sdev of temperature

meanT = meanT/Ntrial;
sdevT = sqrt((sdevT - Ntrial*meanT.^2)/(Ntrial-1));
figure;
bladeplot;