function [npoint,nval,ngrad]=linesearch(fun,cpoint,dir,maxa,cval,cgrad,c1,c2);
%[npoint,nval,ngrad]=linesearch(@fun,cpoint,dir,maxa,cval,cgrad,c1,c2)


%compute directional derivative at npoint
dirderiv=dir'*cgrad;        



if dirderiv>=0 
  disp('Error: no descent direction');
end

% set lower bound on search interval
mina =0.;

% first steplength (heuristic)
a=min([1,maxa]);

% compute first point
npoint=cpoint+a*dir;
[nval,ngrad]=feval(fun,npoint);
ndirderiv = dir'*ngrad;


if (nval<cval+c1*a*dirderiv)  %Is Armijo satisfied? 
  if (ndirderiv>c2*dirderiv)  %Is curvature satified?
    disp(' first point is already ok ');
    return;                   % done
  else
    % "a" seems to be too short, find a valid right end point
    % we do this by interpolation on the last two values
    % but step length should double at least in each iteration
    left = 0;
    fl=cval;
    dfl=dirderiv;
    right = a;
    fr=nval;
    dfr=ndirderiv;
    while (maxa-right>eps(maxa))
      disp(['right of ' num2str(right)]);
      % compute interpolation  coeff0*dx^3+coeff1*dx^2+dfl*dx+fl
      dx=right-left;
      coeff0=((dfl-dfr)*dx^2+2*(fr-dx*dfl-fl)*dx)/(-dx^4);
      coeff1=(dfr-dfl-3*coeff0*dx^2)/2/dx;
      % find the relevant stationary points
      if (abs(coeff0)<eps)  % interpolation is quadratic 
        if (coeff1<eps)     
          disp('linear interpolation');
          a=min([right+2*dx, maxa]);
        else
          disp('quadratic interpolation');
          a=max([right+2*dx, left-dfl/(2*coeff1)]);
          a=min([a, maxa]);
        end
      else
        discr=coeff1^2-3*coeff0*dfl;
        if (discr<0)
          disp('no stationary points');
          a=min([right+2*dx, maxa]);
        else
          a=max([right+2*dx, left+(-coeff1+sqrt(discr))/(3*coeff0)]);
          a=min([a, maxa]);
        end
      end
      % store the gradient of the left end of the new search inteval
      lgrad=ngrad;
      % compute next point
      npoint=cpoint+a*dir;
      [nval,ngrad]=feval(fun,npoint);
      ndirderiv = dir'*ngrad;          
      % before checking the point set the new search interval
      left=right;
      fl=fr;
      dfl=dfr;
      right=a;
      fr=nval;
      dfr=ndirderiv;
      % check the new point for Wolfe
      if (nval<cval+c1*a*dirderiv)  %Is Armijo satisfied? 
        if (ndirderiv>c2*dirderiv)  %Is curvature satified?
          disp(' found point in right end point search ');
          return;
        end
        % still not the right endpoint, go on
      else
        % Armijo is not satified, the right end point is found
        break;
      end
    end  %while (maxa-right>eps(maxa))
    %if the code gets here and Armijo but not curvature is satisfied
    %then maximum stepsize is reached and we cannot go further
    if (nval<cval+c1*a*dirderiv)  %Is Armijo satisfied? 
       disp(' reached maximum step size within Armijo but without curvature');
       return;           
    end
  end 
else
  % Armijo is not satifisfied, the right end point is found; 
  left = 0;
  fl=cval;
  dfl=dirderiv;
  lgrad=cgrad;
  right =a;
  fr=nval;
  dfr=ndirderiv;
end

% search the interval between left and right with the invariant:
% the left point satisfies Armijo but not curvature
% the right point does not satisfy Armijo
while (right-left>eps)
  disp(['interval [' num2str(left) ',' num2str(right) ']']);
  % compute interpolation  coeff0*dx^3+coeff1*dx^2+dfl*dx+fl
  dx=right-left;
  coeff0=((dfl-dfr)*dx^2+2*(fr-dx*dfl-fl)*dx)/(-dx^4);
  coeff1=(dfr-dfl-3*coeff0*dx^2)/2/dx;
  % find the relevant stationary points
  if (abs(coeff0)<eps)  % interpolation is quadratic
    if (coeff1<eps)     
      disp('linear interpolation');
      a=(left+right)/2;
    else
      disp('quadratic interpolation');
      a=left-dfl/(2*coeff1);
    end
  else
    discr=coeff1^2-3*coeff0*dfl;
    if (discr<0)
      disp('no stationary points');
      a=(left+right)/2;
    else
      a=left+(-coeff1+sqrt(discr))/(3*coeff0);
    end
  end
  % ensure minimal progress within the interval
  a=min([a right-dx/10]);
  a=max([a left+dx/10]);
  % compute next point
  npoint=cpoint+a*dir;
  [nval,ngrad]=feval(fun,npoint);
  ndirderiv = dir'*ngrad;          
  % check the new point for Wolfe
  if (nval<cval+c1*a*dirderiv)  %Is Armijo satisfied? 
    if (ndirderiv>c2*dirderiv)  %Is curvature satified?
      disp(' point in interval is ok ');
      return;
    end
    % the left end point moves to a
    left=a;
    fl=nval;
    dfl=ndirderiv;
    lgrad=ngrad;
  else
    % Armijo is not satified, the right end point moves to a
    right=a;
    fr=nval;
    dfr=ndirderiv;
  end
end

%it is unlikely that the code gets here, because the interval
%is smaller than machine precision and no Wolfe point was found
%the left end point should at least satify Armijo if it is not
%the starting point itself. In both cases the left endpoint seems
%to be the better choice
    

disp(' interval too small, return left interval end point ');
npoint=cpoint+left*dir;
nval=fl;
ngrad=lgrad;

