% Perform Beck's sequential function specification method
%   -dimensionless data must be supplied
%   -geometry is the insulated finite plate
%   -uses function fdX22B10T0.m
% dt is the dimensionless time step in the data
%  x is a scalar dimensionless measurement location
%  r is a scalar number of future times for the estimation
%  Y is the vector of measurements beginning at time=dt
% function [qhat,Tsur] = Beck( dt, x, r, Y )
function [qhat,Tsur] = Beck( dt, x, r, Y )
A = 15;  % needed for evaluation of fdX22B10T0
if( nargin == 0 )
   dt = 1;
   x = 0.5;
   r = 3;
   Y = [1:1:10]';
end

N = size(Y,1);   % the total number of observations
if( N == 1 )
   error('Y must be a column vector');
end
time = [dt:dt:N*dt]';
phi = fdX22B10T0( x, time, A ); % get the sensitivity values
denom = phi(1:r)'*phi(1:r);
K = phi(1:r)/denom;
qhat = zeros(N-r, 1);
Mat = zeros(r, N-r);
That = zeros( r, 1 );
% do the first one separately
qhat(1) = ( Y(1:r) )'*K;
for M = 2 : N-r
   % add another column each time to the 'Mat'
   Mat(:,(N-r)-(M-2)) = phi( M : M+r-1) - phi( M-1 : M+r-2 );
   % compute the influence of all the previous q's
   That = Mat(:, (N-r)-(M-2):(N-r) ) * qhat(1:M-1);
   % now use the Beck formula to compute the next qhat
   qhat(M) = ( Y(M:(M+r-1)) - That(1:r) )'*K;
end
%
% now we know all the q components.  
% Compute the surface temperature via Beck's Eq. 3.2.12
phi0 = fdX22B10T0( 0, time, A ); % get the sensitivity values @ x=0
dphi0 = diff([0;phi0]);
Tsur = zeros(N-r, 1);
Tsur(1) = qhat(1)*phi0(1);  % special case for first one
for M = 2 : N-r
   Tsur(M) = 0;
   for n = 1:M
      Tsur(M) = Tsur(M) + qhat(n)*dphi0(M-n+1);
   end
end



   


