% solve whole domain IHCP with first order Tikhonov regularization
% assume dimensionless variables
A = 15;
x = [ 1. ];   % measurement locations
t = [ -0.24: 0.06: 1.8 ]';
dt = t(2) - t(1);
rng(129);
noise = randn( size(t) );
sigma_pct = 0.005;
% T = T_1pulse(t, 0, 1.2, x )';
% T = round(T,6);      % round to 6 decimal places
T = [ 0.000000 0.000000 0.000000 0.000000 0.000000 0.000007 ...
      0.000374 0.002171 0.006323 0.013381 0.023656 0.037319 ...
      0.054465 0.075145 0.099389 0.127201 0.157880 0.189293 ...
      0.219594 0.247679 0.272930 0.295005 0.313715 0.328954 ...
      0.340665 0.348823 0.353763 0.356545 0.358089 0.358943 ...
      0.359415 0.359677 0.359821 0.359901 0.359945 ]';
sigma = sigma_pct*max(T);
Y = T + sigma*noise;
% can't have negative times
tzero = t - t(1);
phi = fdX22B10T0( x, tzero, A );
dphi = diff( phi );
n = length( dphi );  % exclude first data point (I.C.)
X = zeros( n, n );
for ic = 1:n
    X( ic:n, ic ) = dphi(1:n-ic+1);
end
cond(X); 
H = -eye(n,n);
for in = 1:n-1
    H(in,in+1)=1;
end
H(n,n)=1; H(n,n-1)=-1;
alpha_vec = [ 0 1e-6 1e-2 1e-0 ];
q_exact = 0.6*q_1pulse( t(2:n+1)-dt/2, 0, 1.2);
figure(99); clf;
hold on
legstr = [ "ro", "b+", "md", "ks", "gv" ];
Results = zeros( n, length(alpha_vec) );
legent = string(length(alpha_vec))
for ia = 1:length(alpha_vec)
    alpha = alpha_vec(ia);
    qhat = (X'*X + alpha*H'*H) \ X'*Y(2:n+1);  % treat first Y(1) as initial condition
    q_RMS = sqrt((qhat - q_exact)'*(qhat-q_exact)/n);
    fprintf(' %10.3e   %10.3f\n', alpha, q_RMS);
    plot( t(2:end) - dt/2, qhat, legstr(ia), 'linewidth', 1.0 )
    Results(:,ia) = qhat;
    if( alpha==0)
        legent(ia) = sprintf('$\\alpha$=%3.1f', alpha);
    else
        legent(ia) = sprintf('$\\alpha$=%7.1e', alpha);
    end
end
plot( t, 0.6*q_1pulse( t, 0, 1.2), 'k--', 'linewidth', 1.0 )
Results

xlabel('time'); ylabel('Heat flux')
legend([legent 'Exact'], 'Interpreter', 'latex')
axis( [ -0.5 2.0 -.2 0.8 ] )

    fig=figure(99);
    print(fig,fullfile(pwd, 'Figure4_5.tif'),'-dtiff','-r600')

