% eda07_01
% example of a filter problem
% using impulse repsonse of a heat-producing layer
% with solution by damped least squares
% Note: symbol q for theta, temperture

% Note: eda07_01 uses figure 1,2
%       eda07_02 uses figure 3,4
%   and eda07_03 uses figure 5,6
% to make a comparison easy, so don't
% delete the figures until after running all
% three of these scripts

% time setup for g (length M)
M=1024;
N=M;
Dtdays = 0.5;
tdays = Dtdays*[0:N-1]';
Dtseconds = Dtdays*3600*24;
tseconds = Dtseconds*tdays;

% materal constants
rho = 1500; % kg/m3
cp = 800; % J / kg-K
kappa = 1.25e-6; % m2/s
z = 1; % m

% impulse response
g = zeros(M,1);
g(1)=0.0; % manually set first value
g(2:M) = (1/(rho*cp)) * (Dtseconds/sqrt(2*pi)) .* ...
    (1./sqrt(2*kappa*tseconds(2:M))) .* ...
    exp( -0.5*(z^2)./ (2*kappa*tseconds(2:M)));

% true heat production
htrue = zeros(N,1);
t1 = Dtdays*N/4;
t2 = Dtdays*7*N/16;
sd = Dtdays*N/256;
Ah=10;
htrue = Ah*exp(-((tdays(1:N)-t1).^2) / (2*sd^2) ) + 0.5*Ah*exp(-((tdays(1:N)-t2).^2) / (2*sd^2) );

% predict true temperature
tmp = conv( g, htrue);
qtrue=tmp(1:N);

% plot impulse response
figure(1);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(M), 0, 1.1*max(g)] );
plot(tdays(1:M),g,'k-','LineWidth',2);
xlabel('time, days');
ylabel('g(t)');

% plot true heat production
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:N),htrue,'k-','LineWidth',2);
xlabel('time, days');
ylabel('htrue(t)');

% plot true temperature
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), 0, 1.1*max(qtrue)] );
plot(tdays(1:N),qtrue,'k-','LineWidth',2);
xlabel('time, days');
ylabel('qtrue1(t)');

% matrix version of impulse response
G = toeplitz([g', zeros(N-M,1)']', [g(1), zeros(1,M-1)] );
qtrue2 = G*htrue;

% create simulated observations, qobs = qtrue+noise
sigmad=0.001*max(abs(qtrue));
qobs=qtrue+random('Normal',0,sigmad,N,1);

% plot simulated observations
figure(2);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), 0, 1.1*max(qobs)] );
plot(tdays(1:N),qobs,'k-','LineWidth',2);
xlabel('time, days');
ylabel('qobs(t)');

% solve for h given qobs and g using least squared
GTG=G'*G;
GTGmax = max(max(abs(GTG)));
e2=1e-4*GTGmax;
hest1=(GTG+e2*eye(M,M))\(G'*qobs);

% plot true heat production
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:N),htrue,'k-','LineWidth',2);
xlabel('time, days');
ylabel('htrue(t)');

% plot estimated heat production
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(N), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:N),hest1,'k-','LineWidth',2);
xlabel('time, days');
ylabel('hest(t)');

