% Example2_15.m Weighted least squares steady state thermal cond.
clear all; format short g
%Initial data
n=11; m=2; x=(0:0.005:.05)'; fsave=[400; 200]; T0=[80; 40]; sig=[0.10; 20]; sig2=sig.*sig; kr=1/0.5;
rng(129);  % let's always use the same sequence of random numbers
Ntrials = 10;  % number of MonteCarlo trials
results = zeros(Ntrials, 4 ); % save them
figure(1); clf; figure(2); clf;
fprintf('Case Exp1   Exp2  Avg1,2  K_WLS\n')
for iM = 1:Ntrials
    for ii=1:m
       ep=randn(n,1); y(:,ii)=T0(ii)-fsave(ii)*kr*x+ep*sig(1);
    end
    f=fsave+sig(2)*randn(m,1);ff=f;
    %Calculate the average conductivity treating experiments separately.
    X=[ones(n,1) -f(1)*x]; b(:,1)=X\y(:,1); X=[ones(n,1) -f(2)*x]; b(:,2)=X\y(:,2);bavg=(b(:,1)+b(:,2))/2;
    %b(2,1) is K for first experiment. b(2,2) is K for 2nd experiment.
    %end calculation of average conductivity, Kavg=bavg(2).
    %Start calculation of K assuming heat flux values
    yb1=sum(y(:,1))/n;yb2=sum(y(:,2))/n;xb=sum(x)/n; %ik is an iteration index; 2 iterations below
    bw = [ y(1,1) f(1)  y(1,2)  f(2) ]';  % just initialize these to the measured values
    for ik=1:2
        Xk(:,3) = [ -f(1)*x' -f(2)*x' ]';
        Xk(:,1) = [ ones(1,n) zeros(1,n)]';
        Xk(:,2) = [ zeros(1,n) ones(1,n)]';
        yk = [ y(:,1); y(:,2) ];
        Khat_foo = Xk \ yk;  % include the T0 vector we don't need
        Khat(ik) = Khat_foo(3);
        %End calculation of K assuming heat flux values
        %Start calculation of T0(1), q(1),T0(2) and q(2).
        XX(:,1)=[0 0 ones(1,n) zeros(1,n)]'; 
        XX(:,3)=[0 0 zeros(1,n) ones(1,n) ]';
        XX(:,2)=[1 0 -Khat(ik)*x' zeros(1,n)]'; 
        XX(:,4)=[0 1 zeros(1,n) -Khat(ik)*x']';
        IT(ik)=ik; a=[sig2(2), sig2(2), sig2(1)*ones(1,2*n)]; 
        Phi=diag(a); 
        ya=[ff' y(:,1)' y(:,2)']';
        bw=lscov(XX,ya,Phi); %End one iteration for estimating T0(1), q(1),T0(2) and q(2).
        f(1)=bw(2);f(2)=bw(4);
     end
    if( Ntrials <= 500 )
        figure(1); plot(iM, f(1), 'bo', iM, f(2), 'rd', 'linewidth', 1.0 ); hold on;
        figure(2); plot(iM, Khat(ik), 'bo', 'linewidth', 1.0 ); hold on;
    end
     results( iM, : ) = [ b(2,1) b(2,2) bavg(2) Khat(1) ];
     if( Ntrials <= 10 )
        fprintf(' %2d', iM );
        fprintf(' %6.4f', results(iM, :) );
        fprintf('\n')
     end
end
figure(1); xlabel('Monte Carlo Trial Number'); ylabel('Flux, Btu/hr/ft^2')
figure(2); xlabel('Monte Carlo Trial Number'); ylabel('Reciprocal Thermal Conductivity, ft-hr-F/Btu')
axis( [ 1 Ntrials 1.5 2.5 ] )
fprintf('Averages for %3d trials:\n', Ntrials);
fprintf(' %6.4f', mean(results) );
fprintf('\n')
fprintf('Standard deviations for %3d trials:\n', Ntrials);
fprintf(' %6.4f', std(results) );
fprintf('\n')
