% This file synthesizes measurement data for the Inverse Heat Conduction
% Problem
%
% u_t = u_xx,  0 < x < L
% 
% Noisy measurement of the temperature distribution at time Tmax 
% (and x values between x=0 and x=L) is created and saved to file data.mat.
%
% To avoid inverse crime, we use different discretization for simulating
% data than we use for reconstruction.
%
% This method is based on the simple explicit method found, for example, in Ames.
%
% Jennifer Mueller and Samuli Siltanen, October 2012

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

% Numerical parameters  
L    = pi;
D=1;
%Tmax = .4  % use xdim = 24, tdim =108
Tmax = 1    % use xdim = 50, tdim = 500
%Tmax = .1;   % use xdim = 20, tdim = 50
%Tmax = .2;   % use xdim = 40, tdim = 100
%Tmax = .1;   % use xdim = 48, tdim = 216
xdim = 24;
tdim = 108; 

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

% Define the actual initial temperature distribution
Z = zeros(tdim,xdim);   % allocate space for the solution as a matrix
Z(1,:) = 10*sin(2*x);   % the initial temperature
Z(:,1) = 0;             % the BC u(0,t)
Z(:,xdim) = 0;          % the BC u(L,t)

figure(1)
clf
initial = Z(1,:); 
plot(x,initial,'k')
title('Actual initial temperature distribution')


% Synthesize data by time stepping up to Tmax. Rows represent a fixed
% time, columns represent spatial location
v=2:(xdim-1);
for tin=2:tdim
	tip = tin-1;
	Z(tin,v) = Z(tip,v)+(dt/dx^2)*(Z(tip,v+1)-2.*Z(tip,v)+Z(tip,v-1));
end

%%%Matrix method
% Synthesize data by time stepping up to Tmax. 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
for tin=2:tdim
	tip = tin-1;
	ZM(v,tin) =  A*ZM(v,tip);  % If BCs are not zero, need to  add a column vector to the result
    
end
uDatM=ZM(:,tdim);


% Define data at the final time T
uDat = Z(tdim,:);
% Add noise
uDat = uDat + noiselevel*max(abs(initial))*randn(size(uDat));

%  Actual solution for u0 = 10 sin 2x is u(x,t) = 10*exp(-4t)*sin 2x
uAct = 10*exp(-4*Tmax)*sin(2*x);

% compute the relative sup norm error
sn = max(uAct)
supE=(norm(Z(tdim,:) - uAct,inf))/sn


% Show simulated data
% Plot parameters
fsize  = 10;  % This controls the fontsize
lwidth = .5;  % This is the linewidth
msize = 7;

figure(11)
hold on
plot(x,Z(1,:),'k--','linewidth',1.2)
plot(x,uDat,'k.','MarkerSize',11)
plot(x,Z(tdim,:),'k','linewidth',1)
axis([0 pi -10 10])
box off
set(gca,'ytick',[-10,-5,0,5,10],'fontsize',15)
set(gca,'xtick',[0,1.57,3.14],'fontsize',15)

%break

% figure(2)
% hold on
% plot(x,uDat,'b')
% plot(x,uAct,'--')
% plot(x,uDatM,'r')
% title('Data at time T (solid)')

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

figure(3)
mesh(t,x,Z','linewidth',1)
axis([0 1 0 pi -10 10])
colormap(cool)
xlabel('t','fontsize',18)
ylabel('x','fontsize',18)
title('Time evolution from initial data','fontsize',18)
set(gca,'xtick',[0,0.5,1],'fontsize',15)
set(gca,'ytick',[0,1.57,3.14],'fontsize',15)



% figure(4)
% hold on
% plot(x,uDat-uAct,'b')
% title('Difference between computed and actual solution')
% 
% 
%save heat_cond_data_Tp4 uDat xdim tdim Tmax L

