function [xnew,gradnew,alpha,p] = newton(searchmtd,fnhess,fngrad,fn,gradold,xold)
% usage: [xnew,gradnew,alpha,p] = newton(searchmtd,fnhess,fngrad,fn,gradold,xold)
% description: does a single step of Newton's method. Returns
% new x, its gradient, step length and search direction.
% Input parameters:
% searchmtd: line search method
% fnhess: function hessian name
% fngrad: function grad name
% fn: function name
% gradold: gradient at starting point
% xold: starting point

% global variables:
global parms;

% local variables:
% xnew: new x value
% gradnew: gradf(xnew)
% alpha: step length used
% p: descent direction used

p = -feval(fnhess,xold)\gradold;
if (p'*gradold >= 0) % not a descent direction? so restart
  disp('Restart due to non-descent direction');
  p = -gradold;
end
alpha = feval(searchmtd,fngrad,fn,xold,p);
xnew = xold + alpha*p;
gradnew = feval(fngrad,xnew);