% This file solves the Inverse Heat Conduction 
% Problem naively by finite differences
%
% u_t = u_xx,  0 < x < L
% u(x,Tmax) is given as noisy discrete data
% 
% Noisy measurement of the temperature distribution at time Tmax 
% (and x values between x=0 and x=L) is loaded
%
% To avoid inverse crime, we use different discretization for simulating
% data than we use for reconstruction.
%
% Jennifer Mueller and Samuli Siltanen, October 2012

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compare to inverting the forward matrix
%  This file commits an inverse crime

% Load data computed by InverseHeatCondDataSimulator
load heat_cond_data_Tp4

D=1;
x    = linspace(0,L,xdim);
t    = linspace(0,Tmax,tdim);
xn = x(2:xdim-1);
dt   = t(2)-t(1)
dx   = x(2)-x(1)


%%% Create forward matrix. Rows represent a fixed
% time, columns represent spatial location
ZM = zeros(xdim,tdim);   % allocate space for the solution as a matrix
ZM(:,1) = 10*sin(2*x);   % the initial temperature
ZM(1,:) = 0;             % the BC u(0,t)
ZM(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

cond(A)

% add noise (The saved data is noise-free)
n=length(uDat);
noiselevel = 0.0001  % 0.01 is 1 percent noise
uDat = uDat + noiselevel*max(abs(uDat))*randn(size(uDat));
hold on
plot(x,uDat,'r.')
title('Final temperature distribution (data). Black ideal, Red noisy.')

uDat_m=uDat(2:end-1);

num_of_its =4;
ss=tdim-num_of_its;
Back_step=inv(A)*uDat_m';
figure
plot(xn,Back_step,'.')

for ii=1:(tdim-ss)
  Back_step_new=inv(A)*Back_step;
  Back_step=Back_step_new;
end

% Plot parameters
fsize  = 10;  % This controls the fontsize
lwidth = .5;  % This is the linewidth

% Show reconstructed heat profile
figure
plot(xn,Back_step,'k-o','linewidth',lwidth)
axis([0 pi -3 3])
title('t_0=0.3888','fontsize',18)
set(gca,'ytick',[-2.5,0,2.5],'fontsize',15)
set(gca,'xtick',[0,1.57,3.14],'fontsize',15)
%set(gca,'PlotBoxAspectRatio' , [2 1 1])

box off
hold on
% Compare to the actual solution

plot_time = dt*(ss)
uAct = 10*exp(-4*plot_time)*sin(2*xn);
plot(xn,uAct,'k','linewidth',1)
%title('Backward evolution after 4 time steps (dots) and actual solution (solid)')


