%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  function [ H ] = hessian( x, t, y )
%
%       Calculates the Hessian corresponding to Example 2.1
%       in J. Nocedal and S. J. Wright, "Numerical Optimization",
%       Springer-Verlag, 1999
%
%  Parameters
%
%  On entry:
%
%     x        vector nx1 with approximated parameters.
%
%     t        vector mx1 with times of the observations.
%
%     y        vector mx1 with measured observations. 
%
%  On return:
%
%     H        matrix nxn with Hessian of the objective function.
%
%
%  Written by L. N. Vicente.
%  February 21, 2000.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  function [ H ] = hessian( x, t, y )
%
% Obtain dimension m.
  m = size(t,1);
%
% Calculate the residual vector at current x.
  r = y - x(1) - x(2)*exp(-(x(3)-t).^2/x(4)) - x(5)*cos(x(6).*t);
%
% Calculate the Jacobian of the residual vector.
  J = [ -ones(m,1) ...
        -exp(-(x(3)-t).^2/x(4)) ...
        2*x(2)/x(4)*exp(-(x(3)-t).^2/x(4)).*(x(3)-t) ...
        x(2)/(x(4)^2)*exp(-(x(3)-t).^2/x(4)).*((x(3)-t).^2) ...
        -cos(x(6).*t) ...
        x(5)*sin(x(6).*t).*t ];
%
% Calculate the Hessian of the residual vector.
%
  H2(1,:) = 0;
  H2(:,1) = 0;
%
  H2(2,2) = 0;
  H2(2,3) = sum( 2/x(4)*exp(-(x(3)-t).^2/x(4)).*(x(3)-t).*r );
  H2(3,2) = H2(2,3);
  H2(2,4) = sum( -(x(4)^2)*exp(-(x(3)-t).^2/x(4)).*(x(3)-t).^2.*r );
  H2(4,2) = H2(2,4);
  H2(2,5) = 0;
  H2(5,2) = H2(2,5);
  H2(2,6) = 0;
  H2(6,2) = H2(2,6);
%
  H2(3,3) = sum( -4*x(2)/(x(4)^2)*exp(-(x(3)-t).^2/x(4)).*((x(3)-t).^2).*r + ...
                  2*x(2)/x(4)*exp(-(x(3)-t).^2/x(4)).*r );
  H2(3,4) = sum( -2*x(2)/(x(4)^2)*exp(-(x(3)-t).^2/x(4)).*(x(3)-t).*r + ...
                  2*x(2)/(x(4)^3)*exp(-(x(3)-t).^2/x(4)).*((x(3)-t).^3).*r );
  H2(4,3) = H2(3,4);
  H2(3,5) = 0;
  H2(5,3) = H2(3,5);
  H2(3,6) = 0;
  H2(6,3) = H2(3,6);
%
  H2(4,4) = sum( -2*x(2)/(x(4)^3)*exp(-(x(3)-t).^2/x(4)).*(x(3)-t).^2.*r + ...
                 x(2)/(x(4)^3)*exp(-(x(3)-t).^2/x(4)).*(x(3)-t).^4.*r );
  H2(4,5) = 0;
  H2(5,4) = H2(4,5);
  H2(4,6) = 0;
  H2(6,4) = H2(4,6);
%
  H2(5,5) = 0;  
  H2(5,6) = sum( sin(x(6).*t).*t.*r );  
  H2(6,5) = H2(5,6);  
  H2(6,6) = sum( x(5)*cos(x(6).*t).*(t.^2).*r );
%
  H1 = J'*J; 
  H  = H1 + H2;
%
% End of function hessian.
