% 2-D steady-state / inverse conduction
% Adjoint Method
clear
global A B Y C bs bair w T
 Tair=0;w=[0; 0];
 bs=[-2; -1; -1; 0; -1; 0; -1; 0];
 bair=[0; 0; 0; 0; 0; 0; Tair; Tair];
 Y=[0.9258; 0.6847; 0.1953];
 C=[0 1 0 0 0 0 0 0;
    0 0 0 0 1 0 0 0;
    0 0 0 0 0 0 0 1];
 
 for i=1:41
    xx(i)=0.+(i-1)*0.05;
    for j=1:41
      yy(j)=1+(j-1)*0.1;
      B=[xx(i);yy(j)];
      Z(j,i)=phi(0);
end
end
figure(1); clf;
contour(xx,yy,Z,10);hold on

clear xx,clear yy;


B=[2; 1];
xx(1)=B(1);yy(1)=B(2);

for k=2:20
% compute the LS criterion

s(k-1)=phi(0);

% compute the gradient
 psi=inv(A')*C'*(Y-C*T);
 g=-[bs'*psi;2*(psi(7)*(T(7)-Tair)+psi(8)*(T(8)-Tair))];
 ng=g'*g;
 
% determine the descent direction
 if mod(k,2)==0
    w=-g;
 else  
    w=-g+ng/ng1*w1;
 end
 
% determine the optimal descent length (line search)
 ro=fminbnd(@phi,0,1e5);
 
% new iterate
 B=B+ro*w;
 g1=g;w1=w;ng1=ng;
 xx(k)=B(1);yy(k)=B(2);
 
 end

plot(xx,yy,'lineWidth',2)
xtickformat('%0.1f')
ytickformat('%0.1f')
xlabel('\beta_1=T_s')
ylabel('\beta_2 = Bi')

figure(2); clf;
semilogy(s)
xlabel('iteration')
ylabel('S')
grid on

% fig=figure(2);
% print(fig,fullfile(pwd, 'Figure6_24.tif'),'-dtiff','-r600')
% fig=figure(1);
% print(fig,fullfile(pwd, 'Figure6_25.tif'),'-dtiff','-r600')


function y=phi(r)
global A B Y C bs bair w  T
BW=B+r*w;
b=BW(1)*bs-2*BW(2)*bair;
A=  ...
[-4   1  1  0  0  0  0  0; ...
  2  -4  0  1  0  0  0  0; ...
  1   0 -4  1  1  0  0  0; ...
  0   1  2 -4  0  1  0  0; ...
  0   0  1  0 -4  1  1  0; ...
  0   0  0  1  2 -4  0  1; ...
  0   0  0  0  2  0 -4-2*BW(2) 1; ...
  0   0  0  0  0  2  2 -4-2*BW(2)];
 T=A\b; e=Y-C*T;
 y=0.5*(e'*e);

end