% Inspired by Beck's SFSM with a temporary assumption of quadratic form for the unknown function
%   -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] = SFS_Linear( dt, x, r, Y )
function [qhat,Tsur] = SFS_Quad( dt, x, r, Y )
AA = 15;  % 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 = (0:dt:N*dt)';          
phi = fdX22B10T0( x, time, AA ); % get the sensitivity values

A = zeros(r,1);
B = zeros(r,1);
C = zeros(r,1);
qhat = zeros(N-r, 1);

denom = 0;
dphi = diff(phi);
 for i = 1:r
   for j = 0:(i-1)

        A(i) = A(i) + (j+1)*(j+2)*dphi(i-j)/2;
        B(i) = B(i) - j*(j+2)*dphi(i-j);
        C(i) = C(i) + j*(j+1)*dphi(i-j)/2;
   end 
   denom = denom + A(i)*A(i);
 end

% [A B C  ]
% [phi [0; dphi]]

for M = 3:N-r
    % create a vector for estimated temperature at each step 
    estimator = 0;
    for i = 1:r
        T_hat_qM0 = 0;
        for k = 1: M-1
            T_hat_qM0 = T_hat_qM0 + qhat(k)*(phi(M+i-k)- phi(M+i-k-1));
        end
        estimator = estimator + A(i)*(Y(M+i-1) - T_hat_qM0 - qhat(M-1)*B(i) - qhat(M-2)*C(i));
    end
     % Compute the next qhat
    qhat(M) = estimator/denom;        
end  

phi0 = fdX22B10T0(0, time, AA);
dphi0 = diff([phi0]);
Tsur = zeros(N-r, 1);
for M = 1 : N-r
   Tsur(M) = 0;
   for n = 1:M
      Tsur(M) = Tsur(M) + qhat(n)*dphi0(M-n+1);
   end
end

end