% eda09_06
% example of a filter problem
% using impulse repsonse of a heat-producing layer
% with solution via generalized least squares
% using bicongugate gradient method
% with the prior information of smoothness
% Note: symbol q for theta, temperture

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

clear a H;
global a H;

% 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(N,1);
g(1)=0.0; % manually set first value
g(2:N) = (1/(rho*cp)) * (Dtseconds/sqrt(2*pi)) .* ...
    (1./sqrt(2*kappa*tseconds(2:N))) .* ...
    exp( -0.5*(z^2)./ (2*kappa*tseconds(2:N)));

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

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

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

% solve for h given qobs and g using generalized least squares
% with a prior smoothness constraint using biconjugate gradient

% autocorrelation, a, of g. Just keep positive values
al = xcorr(g);
Na = length(al);
a = al((Na+1)/2: Na);

% cross-correlation, c, of qobs2 and g. Just keep positive values
cl = xcorr(qobs2, g);
Nc = length(cl);
c = cl((Nc+1)/2: Nc);

% set up H and h and f
e=1*max(abs(g));
K=M;
L=N+K;
H=spalloc(K,M,3*K);
H(1,1)=e; H(1,2)=-e;
for j = [2:K]
    H(j,j-1)=e;
end
for j = [1:K]
    H(j,j)=-2*e;
end
for j = [1:K-1]
    H(j,j+1)=e;
end
H(K,K-1)=e; H(K,K)=-e;
h=zeros(K,1);

% set up F'f = GT qobs + HT h
% GT qobs is c=qobs2*g
HTh = H'* h;
FTf = c + HTh;

% solve
hest3 = bicg( @autofun, FTf, 1e-10, 3*L);

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

% plot true heat production
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [tdays(1), tdays(M), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:M),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(M), -1.1*max(htrue), 1.1*max(htrue)] );
plot(tdays(1:M),hest3,'k-','LineWidth',2);
xlabel('time, days');
ylabel('hest(t)');
