%
% function [f2,x2,y2,M,Mtoobig,sing] =
%           NextZero(func,x1,y1,tol,maxM,stepD)
%
% This function uses Newton's method to look for a zero of a function
% f(x,y): It assumes you start with a particular zero [x1,y1]
% and it looks for another zero about a distance stepD away.
%
% inputs: func:	a string with the name of a function of the
%	        form [f,dfdx,dfdy] = func(x,y)
%         x1: initial x-value
%         y1: initial y-value
%         tol: tolerance - stops when it finds a point [x,y] with
%              |f(x,y)| <= tol
%         maxM: maximum number of Newton steps to allow - stops
%              even if it hasn't found a root if it exceeds maxM
%         stepD: the approximate distance away from 
%              [x1,y1] to look for the next zero - this
%              is a signed distance - if > 0 we look in the
%              direction of grad(f) rotated counter-clockwise by
%              pi/2, and if < 0 we look in the direction of grad(f)
%              rotated clockwise by pi/2
%
% outputs: f2: value of f at the (approximate) zero it found
%          x2: x-coordinate of the zero it found
%          y2: y-coordinate of the zero it found
%          M: number of Newton steps used
%          Mtoobig: 1 if it stopped because M exceeded
%              maxM, 0 otherwise
%          sing: 1 if it stopped because it encountered
%              a singularity: dfdx = 0 = dfdy (0 otherwise)

function [f2,x2,y2,M,Mtoobig,sing] = ...
        NextZero(func,x1,y1,tol,maxM,stepD)

% store initial values
[f1,dfdx1,dfdy1] = eval([func,'(x1,y1)']);
normgrad = sqrt(dfdx1^2 + dfdy1^2); % norm of gradient vector

% set flags and exit if singularity at start
if (normgrad==0)
   Mtoobig = 0; sing = 1;
   f2 = 0; x2 = 0; y2 = 0; M = 0;
   return;
end

% scale gradient vector to have length 1
unitdfdx = dfdx1/normgrad;
unitdfdy = dfdy1/normgrad;

% find point distance stepD away from (x1,y1) in
% direction perpendicular to gradient vector
x = x1 + stepD*(-unitdfdy);
y = y1 + stepD*unitdfdx;

% main Newton loop
M = 0; s = 0;
[f,dfdx,dfdy] = eval([func,'(x,y)']);
x2 = x; y2 = y;
while (M<maxM) & (abs(f)>tol) & ((dfdx~=0) | (dfdy~=0))
    s = s - f/(unitdfdx*dfdx + unitdfdy*dfdy);
    x2 = x + s*unitdfdx;
    y2 = y + s*unitdfdy;
    [f,dfdx,dfdy] = eval([func,'(x2,y2)']);
    M = M + 1;
end

% prepare return values and flags
f2 = f;
if (M>=maxM) & (abs(f)>tol)
   Mtoobig = 1;
else
   Mtoobig = 0;
end
if (dfdx==0) & (dfdy==0) & (abs(f)>tol)
   sing = 1;
else
   sing = 0;
end