% This file solves the Inverse Heat Conduction Problem with naive 
% inversion.
%
% NOTE: This file commits an inverse crime since we use the same
% discretization both for data simulation and inversion.
% 
% Jennifer Mueller, Samuli Siltanen and Esa Niemi.


% Choose relative amplitude of noise (0.01 equals 1 percent noise)
noiselevel = 0.0;

% Numerical parameters
L    = pi;
D    = 1;


Tmax = .1;   

% if Tmax = .1;   % use xdim = 20, tdim = 50
% if Tmax = .2;   % use xdim = 40, tdim = 100

if (Tmax-0.1 < 1e-8)
    xdim = 20;
    tdim = 50;
elseif (Tmax-0.2 < 1e-8)
    xdim = 40;
    tdim = 100;
end


x    = linspace(0,L,xdim);
t    = linspace(0,Tmax,tdim);
dt   = t(2)-t(1);
dx   = x(2)-x(1);
% check the stability criterion dt/(dx^2)<.5
if (dt/(dx^2)>.5)
    display('Stability criterion violated')
    display('dt/(dx^2) = ')
    dt/(dx^2)
end


% Synthesize data by time stepping up to Tmax. Rows represent a fixed
% time, columns represent spatial location
ZM1 = zeros(xdim,tdim);   % allocate space for the solution as a matrix
ZM1(:,1) = 10*sin(2*x);   % the initial temperature
ZM1(1,:) = 0;             % the BC u(0,t)
ZM1(xdim,:) = 0;          % the BC u(L,t)
s = D*(dt/dx^2);
n=length(x) - 2;
A = (1-2*s)*eye(n,n);
for ii =1:n-1
    A(ii,ii+1)=s;
end
for ii = 1:n-1
    A(ii+1,ii)=s;
end
v=2:(xdim-1);
for tin=2:tdim
	tip = tin-1;
	ZM1(v,tin) =  A*ZM1(v,tip);  % If BCs are not zero, need to  add a column vector to the result
end

% Add noise to the data
ZM1(:,end) = ZM1(:,end)+noiselevel*max(abs(ZM1(:,end)))*rand(size(ZM1(:,end)));


% Next we employ a naive inversion (and commit an inverse crime):
reco = ZM1(2:end-1,end);
for kkk=1:tdim-1
    reco = A\reco;
end
reco = [0; reco; 0];

% Plot the reconstruction
figure(1)
clf
plot(x,ZM1(:,1),'r',x,ZM1(:,end),'b',x,reco,'--g');
legend('initial temp. distr.','final temp. distr.', 'reconstruction')