% Example2_29.m
%Program for Example 2.29, For heat losses, find q and initial T
%Sequential linear estimation for m=1.
clear all; format short g
load Example2_29_data 
% note the "data" are in dimensioness form...
X2=TX22;   % two columns:  x=0 and x=L
X3=TX23;   
T0 = 20; q0=50; % the true parameters
L = 0.05; k = 0.05; alpha = 25e-6;  h=0.1; % other parameters from the problem
%mm gives location of columns. mm=1 for x=0, mm=2 for x=L
mm=[1; 2]; n=length(X2); p=2; mu=zeros(p,1); %mu is the prior parameter estimate vector
sigma=0.5; sig=sigma*ones(n,1); Pz=1.0e+5 ;
rng(129);  % always use the same noise
for jj=1:2 %jj =1 is for x = 0; jj = 2 is for x = L
   y=T0+q0*L/k*X3(:,mm(jj)); % exact measurements from X3 case
%   y=T0+q0*X3(:,mm(jj)) + randn(n,1)*sigma; ; % noisy X3 measurements
%   y=T0+q0*X2(:,mm(jj)) + randn(n,1)*sigma; % noise X2 measurements
   X=[ones(n,1) X2(:,mm(jj))*L/k ];  % model is T0 + q0*L/k*X22
   if jj==1; 
       y1=y; Z1=X;
   else; 
       y2=y; Z2=X; 
   end
   b=mu; P0=Pz*eye(p,p); P=P0; B=mu; %P0 is V sub(beta), prior covariance matrix of parameters
   for ii=1:n;			%SEQUENTIAL
       % begin sequential steps here
        A=P*X(ii,:)';	%eq. (2.9.8a)	%SEQUENTIAL
        Delta=sig(ii)^2+X(ii,:)*A;	%eq. (2.9.8)	%SEQUENTIAL
        K=A/Delta;	%eq. (2.9.8c)	%SEQUENTIAL
        e=y(ii)-X(ii,:)*b;	%eq. (2.9.8d)	%SEQUENTIAL
        b=b+K*e;	%eq. (2.9.8e)	%SEQUENTIAL
        P=P-K*A';	%eq. (2.9.8f)	%SEQUENTIAL
      % end sequential steps 
      % results saving lines below to "end"
       if jj == 1
           B1(ii,:)=[ii b(1) b(2)];
       else
           B2(ii,:)=[ii b(1) b(2)];
       end 
   end
end
B1;t=B1(:,1);
figure(1); clf;
plot(t,B1(:,2),'k',t,B1(:,3),'k',t,B2(:,2),'r',t,B2(:,3),'r', 'linewidth', 1 )
text(40,46,'Est. q, data at x=L')
text(50,49,'Est. q, data at x=0')
text(25,22,'Est. T_0, data at x=0 and x=L')

xlabel('Time, s')
ylabel('Parameter estimates')
figure(2); clf;
residuals(:,1) = y1 - Z1*B1(n,2:3)';
residuals(:,2) = y2 - Z2*B2(n,2:3)';
plot(t,residuals(:,1),'k',t,residuals(:,2),'r', 'linewidth', 1 )
xlabel('Time, s')
ylabel('residuals y - X*b, C')
text(22,-0.1,'Data at x=L')
text(5, -0.4, 'Data at x=0')
% print the 600 dpi figure to tiff format
fig=figure(1);
print(fig,fullfile(pwd, 'Figure2_11.tif'),'-dtiff','-r600')
fig=figure(2);
print(fig,fullfile(pwd, 'Figure2_12.tif'),'-dtiff','-r600')
