% Semi-infinite body
% Inverse input problem
% Integral Equation model
% Conjugate gradient algorithm
clear
dt=0.1;nf=50;
% the discretized heat flux
for n=1:nf
   axet(n)=dt*n;
  if n<11
      q(n)=n;
  elseif  n<31
    q(n)=20-n;
 elseif n<41
     q(n)=-40+n;
   else
      q(n)=0;
   end
end
qexact=q;

% the output measurement
% with additive noise
for n=1:nf
   s=0;
   for i=1:n-1
     s=s+impulse2(n-i)*q(i);
   end
   Y(n)=s;
end
A=1.;
rng(129);  % always use the same noise
noise=A*randn(size(Y));Y=Y+noise;
eps=norm(noise)^2

% conjugate gradient algorithm

q(1:nf)=0;ls=1e9;iter=1;

while ls>eps
   % compute the LS-criterion
   for n=1:nf
      s=0;
      for i=1:n-1
         s=s+impulse2(n-i)*q(i);
      end
      T(n)=s;psi(n)=Y(n)-T(n);
   end
   ls=0.5*norm(psi)^2;
   lsw(iter)=ls;
   fprintf('iter %2d S=%g\n',iter, ls);

   
   % compute the gradient
   for n=1:nf
      s=0;
      for i=n+1:nf
         s=s+impulse2(i-n)*psi(i);
      end
      g(n)=s;
   end
   ng=norm(g)^2;
   if iter==1
      w=-g;
   else
      w=-g+ng/ng1*w1;
   end
   %compute the descent length
   % quadratic case
   for n=1:nf
      s=0;
      for i=1:n-1
         s=s+impulse2(n-i)*w(i);
      end
      dT(n)=s;
   end
   
   ro=(g*w')/(norm(dT)^2);
   q=q+ro*w;
   g1=g;w1=w;ng1=ng;
   iter=iter+1;
   
end

figure(1); clf;
semilogy([0:iter-2],lsw, 'linewidth', 1.0)
xlabel('iteration'); ylabel('S=0.5||Y-T||^2')
xticks([0:4])
grid on
figure(2); clf;
plot(axet,qexact,'k-',axet,q,'ro', 'linewidth', 1.0)
legend('exact','estimate')
xlabel('time'); ylabel('Heat Flux')
figure(3); clf;
plot(axet,Y,'k-',axet,T,'ro', 'linewidth', 1.0)
legend('exact','estimate')
xlabel('time'); ylabel('Temperature')

% fig=figure(1);
% print(fig,fullfile(pwd, 'Figure6_26.tif'),'-dtiff','-r600')
% fig=figure(2);
% print(fig,fullfile(pwd, 'Figure6_27.tif'),'-dtiff','-r600')
% fig=figure(3);
% print(fig,fullfile(pwd, 'Figure6_28.tif'),'-dtiff','-r600')

function y=impulse2(n)
% temperature response 
% to the heat flux impulse 
% the semi-infinite heat conduction
dt=0.1;
if n<1
      y=0;
else
     tt=1/(n*dt); y=sqrt(tt)*exp(-tt);
end
end