% Perform Beck's sequential function specification method
%  ** using the temperature based method **
%   -dimensionless data must be supplied
%   -geometry is the insulated finite plate
%   -uses function finite_T.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_T( dt, x, r, Y )
function [qhat,Tsur] = Beck_T( dt, x, r, Y )
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 = fdX12B10T0( x, time ); % get the sensitivity values at the 'x' location
denom = phi(1:r)'*phi(1:r);
K = phi(1:r)/denom;
Tsur = zeros(N-r, 1);
Mat = zeros(r, N-r);
qhat = zeros( N-r, 1 );
That = zeros( r, 1 );
% do the first one separately
Tsur(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) ) * Tsur(1:M-1);
   % now use the Beck formula to compute the next qhat
   Tsur(M) = ( Y(M:(M+r-1)) - That(1:r) )'*K;
end
%
% now we know all the Tsur components.  
% Compute the surface heat flux 
%
% use subdivisions on each time interval to improve evaluation accuracy
ndiv = 20;   % must be even?
time_0 = [ dt/ndiv : dt/ndiv : N*dt ]';
[ foo, deriv0] = fdX12B10T0( 0, time_0 ); % get the sensitivity values @ x=0
deriv0 = - deriv0;  % want dT/dx not -dT/dx
dderiv0 = diff([0;deriv0]);
qhat = zeros(N-r, 1);
for M = 1 : N-r
	qhat(M) = 0;
	MAX = M*ndiv;
   for n = 1:MAX
      T_n = interp1( [0;time(1:N-r)], [0;Tsur], time_0(n));
      qhat(M) = qhat(M) - T_n*dderiv0(MAX-n+1);
   end
end