% [N,M,avgMpt,maxMpt,x,y,maxf,maxD,minD] = ...
%    ZeroCurve(func,x1,y1,stepD,stepcontrol,range,tol,maxN,maxM) 
%
% Uses Newton's method to plot the solution of a 2-d implicit equation
% f(x,y) = 0.
%
% Inputs: func: a string with the name of the function to use,
%                of the form [f,dfdx,dfdy] = func(x,y)
%         x1,y1: initial coordinates, with f(x1,y1)==0
%         stepD: step size - in general we look for points separated
%                by a distance of about stepD
%         stepcontrol : 1 if we want to use step-control, else 0
%         range : [xmin,xmax,ymin,ymax] plot stays in this range
%         tol : the tolerance for stopping Newton's algorithm
%         maxN : the maximum number of points to plot on the curve
%         maxM : the maximum number of Newton iterations to allow
%
% Outputs: N : actual number of points found
%          M : actual number of Newton steps used
%          avgMpt : average number of Newton steps per point
%          maxMpt : maximum number of Newton steps used for a point
%          x,y : 1-by-n arrays : coordinates of points on the curve
%          maxf : maximum value of abs(f) at points found
%          maxD : maximum distance between consecutive points
%          minD : minimum distance between consecutive points
%

function [N,M,avgMpt,maxMpt,x,y,maxf,maxD,minD] = ...
    ZeroCurve(func,x1,y1,stepD,stepcontrol,range,tol,maxN,maxM);

% note on variable names: anything with "N" is a count of points
% anything with "M" is a count of Newton steps, anything with
% "D" is a distance between points.

% set things up
thisx = x1; thisy = y1; x = x1; y = y1;
M = 0; % count of total number of Newton steps
Mpt = 0; % list of num Newton steps per each point
stepDpt = stepD; % list of actual stepsizes used for each point
signedstepD = abs(stepD);
[f1,dfdx,dfdy] = eval([func,'(x1,y1)']);
f = f1;
stop = 0; % stop flag
Nsincex1y1 = 0; % number of points since last at (x1,y1)

% initial error check:
if (abs(f1)>tol)
  disp('ERROR: f(x1,y1) not equal to zero');
  return;
end

% main loop
while (stop==0)
   if (stepcontrol==1)
       [nextf,nextx,nexty,Mthispt,Mtoobig,sing,newsignedstepD] = ...
           NextZeroC(func,thisx,thisy,tol,maxM-M,signedstepD,2*stepD);
       stepDpt = [stepDpt,abs(newsignedstepD)];
       signedstepD = newsignedstepD*2;
       if abs(signedstepD) > abs(stepD)
           signedstepD = sign(signedstepD)*abs(stepD);
       end
   else
       [nextf,nextx,nexty,Mthispt,Mtoobig,sing] = ...
            NextZero(func,thisx,thisy,tol,maxM-M,signedstepD);
   end
   M = M + Mthispt;
   Nsincex1y1 = Nsincex1y1 + 1;
   if (sing==1)
       stop = 1; disp('Singularity encountered, stopping')
       Mpt = [Mpt, Mthispt]; 
   elseif (Mtoobig==1)
       stop = 1; disp('Too many Newton steps, stopping')
       Mpt = [Mpt, Mthispt]; 
   elseif (nextx<range(1)) | (nextx>range(2)) | ...
	    (nexty<range(3)) | (nexty>range(4))
       if (signedstepD>0)
           signedstepD = -abs(stepD);
           thisx = x1; thisy = y1;
           x = fliplr(x); y = fliplr(y); 
           f = fliplr(f); Mpt = fliplr(Mpt);
           Nsincex1y1 = 0; disp('Out of range, back to (x1,y1)')
       else 
           stop = 1; disp('Out of range, stopping')
       end
   elseif ((Nsincex1y1>50) & (Dist(nextx,nexty,x1,y1)<2*abs(signedstepD)))
       x = [x,nextx,x1]; y = [y,nexty,y1]; f = [f,nextf,f1];
       Mpt = [Mpt,Mthispt,0];
       stop = 1; disp('Closed curve, stopping')
   elseif ((size(x,2)==maxN-1) | (M==maxM))
       x = [x,nextx]; y = [y,nexty]; f = [f,nextf]; 
       Mpt = [Mpt, Mthispt];
       stop = 1; disp('Just enough points or Newton steps, stopping')
   else
       x = [x,nextx]; y = [y,nexty]; f = [f,nextf]; 
       Mpt = [Mpt, Mthispt]; 
       thisx = nextx; thisy = nexty;
   end
end

% setup return values
N = size(x,2); 
avgMpt = M/N;  
maxMpt = max(Mpt);
maxf = max(abs(f));
maxD = 0; minD = inf;
for cnt=1:N-1
  D = Dist(x(cnt),y(cnt),x(cnt+1),y(cnt+1));
  maxD = max(maxD,D); minD = min(minD,D);
end

% generate plots
hold off, clf
if (stepcontrol==1 & N>1)
  figure(3)
  semilogy(stepDpt,'k-');
  title('Actual stepsizes used for each point');
  rng = [0,size(stepDpt,2),minD/10,maxD*10];
  axis(rng);
  xlabel([' N = ',int2str(N),',  minD = ',num2str(minD),',  maxD = ',num2str(maxD)])
end
[X,Y] = meshgrid(range(1):(range(2)-range(1))/50:range(2),...
                 range(4):(range(3)-range(4))/50:range(3));
[f,dfdx,dfdy] = eval([func,'(X,Y)']);
figure(1)
hold off, clf, spy(f<0,'g'); hold on, spy(f>=0,'r');
title('Crude expensive plot of f(x,y)=0');

figure(2)
clf,plot(x,y,'kx-',x1,y1,'ko');
title(['f(x,y)=0, drawn with Newtons method and stepcontrol=',int2str(stepcontrol)]);
xlabel(['N = ',int2str(N),',  M = ',int2str(M), ...
        ',  avgMpt = ',num2str(avgMpt),'  maxMpt = ',int2str(maxMpt), ...
        ',  maxf = ',num2str(maxf)]),
axis(range);axis('square');