% compute condition numbers for different delta-t
% using zeroth and first order regularization
% assume dimensionless variables
A = 15;
x = [ 1. ];   % measurement locations

dtvec = [ 0.03 0.06 0.1 0.12 0.15 0.2];
fprintf('With no regularization:\n')
for it = 1:length(dtvec)
    dt = dtvec(it);
    t = [ 0:dt:1.8 ];
    phi = fdX22B10T0( x, t, A );
    dphi = diff( phi );
    n = length( dphi );  % compute also for first data point
    X = zeros( n, n );
    for ic = 1:n
        X( ic:n, ic ) = dphi(1:n-ic+1);
    end
    fprintf(' %6.3f  %10.3e\n', dt, cond(X) )
end
alpha_vec = [ 1e-8 1e-7 1e-6 1e-5 ];
for it = 1:length(dtvec)
    dt = dtvec(it);
    H0_result( it, 1) = dt;
    H1_result( it, 1) = dt;
    t = [ 0:dt:1.8 ];
    phi = fdX22B10T0( x, t, A );
    dphi = diff( phi );
    n = length( dphi );  % compute also for first data point
    X = zeros( n, n );
    for ic = 1:n
        X( ic:n, ic ) = dphi(1:n-ic+1);
    end
    H0 = eye(n,n);
    H1 = -eye(n,n); 
    for i = 1:n-1
        H1(i,i+1 ) = 1;
    end
    H1(n,n)=1; H1(n,n-1)=-1;
    
    H0_result(it,2) = cond(X'*X);
    H1_result(it,2) = cond(X'*X);
    for iT = 1:length(alpha_vec)
        alpha = alpha_vec(iT);
        H0_result(it,iT+2) = cond(X'*X + alpha*H0'*H0);
        H1_result(it,iT+2) = cond(X'*X + alpha*H1'*H1);
    end
end
fprintf('With zeroth order Tikhonov:\n')
fprintf('  dt   ');
fprintf(' alpha=%5.1e', [0 alpha_vec] )
fprintf('\n')
for it = 1:length(dtvec)
    fprintf(' %5.3f', H0_result(it,1))
    fprintf('    %10.3e', H0_result(it,2:end))
    fprintf('\n');
end
fprintf('With first order Tikhonov:\n')
fprintf('  dt   ');
fprintf(' alpha=%5.1e', [0 alpha_vec] )
fprintf('\n')
for it = 1:length(dtvec)
    fprintf(' %5.3f', H1_result(it,1))
    fprintf('    %10.3e', H1_result(it,2:end))
    fprintf('\n');
end
