% solve whole domain IHCP with zeroth order Tikhonov regularization
% assume dimensionless variables
A = 15;
order = 'zero';
tlast = 2.34;   % 1.8 for part a); change to 2.34 for part b)
x = [ 1. ];   % measurement locations
% load HeatEqOut_d   % has t_d, TL_d, and qbc_d  - these are dimensional
t_d = [ -34.0 : 8.5 : 331.5 ]';
TL_d = [ 0.01 -0.04 0.02 -0.04 0.14 0.00 0.02 -0.06 -0.03 0.00 ... 
         0.10  0.28 0.42  0.65 1.00 1.28 1.74  2.07  2.50 3.02 ... 
         3.44  3.83 4.30  4.71 5.02 5.36 5.77  5.91  6.28 6.69 ... 
         6.84  7.13 7.21  7.42 7.55 7.64 7.80  7.87  7.95 8.04 ...
         7.99 8.15 8.20 8.20 ]';
qbc_d = [   0.0    0.0    0.0    0.0    0.0  204.3  701.1 1481.5 ... 
         2822.0 3900.3 4911.1 5860.0 6742.3 7618.6 8578.1 9421.8 ...
         9934.6 9884.1 9200.6 8227.0 7340.0 6872.2 6617.5 6420.9 ...
         6212.5 5922.4 5506.0 5003.5 4471.0 3963.6 3449.3 2917.7 ...
         2424.6 2025.9 1748.9 1539.0 1334.9 1075.2  591.7  189.5 ...
          133.4   36.9    0.0    0.0 ]';

% these are known from the "experiment"
% some dimensions
k = 60.5;          % W/m-K  low carbon steel
alpha = 17.7e-6;   % m^2/s   low carbon steel
q0 = 10000;          % W/m^2
L_d = 0.050;       % m

% make non-dimensional
t = t_d*alpha/L_d^2;
TL = TL_d/(q0*L_d/k);
qbc = qbc_d/q0;

% t = [ -0.24: 0.06: 1.8 ]';
dt = t(2) - t(1);
% find t = 1.8
n = find(t>=tlast,1)
% resize the arrays
tdata = t;
TLdata = TL;
qbcdata = qbc;
clear t TL qbc
t = tdata(1:n); T = TLdata(1:n); qbc = qbcdata(1:n);


Y = T;  % noise already in the data TL data
% 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); 
if ( strcmp(order,'first') )
    H = -eye(n,n);
    for in = 1:n-1
        H(in,in+1)=1;
    end
    H(1,1:3)= [-3 4 -1]/2;
    H(n,n-2:n)= [1 -4 3]/2;
    alpha_vec = [ 5e-3 1e-3 5e-4  ];
    fig_no = 91;
elseif (strcmp(order, 'second') )
    H = -2*eye(n,n);
    for in =2:n-1
        H(in,in+1)=1;
        H(in,in-1)=1;
    end
    H(1,1:4)= [2 -5 4 -1];
    H(n,n-3:n)= [-1 4 -5 2];
    alpha_vec = [ 5e-3 1e-3 5e-4  ];
    fig_no = 92;
else
    H = eye(n,n);
    alpha_vec = [ 5e-3 1e-3 5e-4  ];
    fig_no =90;
end

q_exact = qbc;
figure(fig_no); 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_T = alpha_vec(ia);
    qhat = (X'*X + alpha_T*H'*H) \ X'*Y(2:n+1);  % treat first Y(1) as initial condition
    plot( t(2:end) - dt/2, qhat, legstr(ia), 'linewidth', 1.0 )
    Results(:,ia) = qhat;
    if( alpha_T==0)
        legent(ia) = sprintf('$\\alpha$=%3.1f', alpha_T);
    else
        legent(ia) = sprintf('$\\alpha$=%7.1e', alpha_T);
    end
end
plot( t, q_exact, 'k--', 'linewidth', 1.0 )
text(-0.25,1.,sprintf('Tikhonov %s order',order) )
% Results

xlabel('time'); ylabel('Heat flux')
legend([legent 'Exact'], 'Interpreter', 'latex')
axis( [ -0.5 2.5 -0.2 1.2 ] )
figure(96); clf;
plot(t,T,'k--', t, Y, 'ro' )

switch order
    case 'zero'
        fname = 'Figure4_6c0.tif';
    case 'first'
        fname = 'Figure4_6c1.tif';
    case 'second'
        fname = 'Figure4_6c2.tif';
end

% plot with units
figure(109); clf;
fig_no = 109;  % save this figure
hold on
for ia = 1:length(alpha_vec)
    alpha_T = alpha_vec(ia);
    plot( (t(2:end) - dt/2)*L_d^2/alpha, Results(:,ia)*q0/1000, legstr(ia), 'linewidth', 1.0 )
    if( alpha_T==0)
        legent(ia) = sprintf('$\\alpha$=%3.1f', alpha_T);
    else
        legent(ia) = sprintf('$\\alpha$=%7.1e', alpha_T);
    end
end
plot( t_d, qbc_d/1000, 'k--', 'linewidth', 1.0 )
% text(-0.25,1.,sprintf('Tikhonov %s order',order) )
xlabel('time, seconds'); ylabel('Heat flux, kW/m^2')
legend([legent 'Exact'], 'Interpreter', 'latex')
axis( [ -50 350 -2 12 ] )
text(230,6.5,sprintf('Tikhonov %s order',order) )

% print a tableau
fprintf('%5.2f %5.2f %6.1f\n',[ t_d(1:n+1) Y*q0*L_d/k qbc_d(1:n+1) ]' )

% fig =figure(fig_no);
%      print(fig,fullfile(pwd, fname),'-dtiff','-r600')

