function [x0,y0] = fplot(fname,lims,npts,angl,subdiv)
%FPLOT	Plot a function.
%	FPLOT(FNAME,LIMS) plots the function specified by the text variable
%	FNAME between the x-axis limits specified by LIMS = [XMIN XMAX].
%	For example, FPLOT('sin',[0 10]) is a plot of the sine function.
%	FPLOT(FNAME,LIMS,N) plots the function with a minimum of N samples.
%	The default value for N is 25.
%	FPLOT(FNAME,LIMS,N,ANGLE) plots the function with a minimum of N
%	samples and samples more if contiguous segments are at an angle
%	larger than ANGLE (in degrees).  The default for ANGLE is 10 degrees.
%	FPLOT(FNAME,LIMS,N,ANGLE,SUBDIV) plots the function with a minimum
%	of N samples, no more than ANGLE degrees bends in the curve, but no
%	more than SUBDIV iterations to fill in the rapidly changing portions.
%	The default for SUBDIV is 20.
%	[X,Y] = fplot(FNAME,LIMS,...) returns the abscissae and ordinates for
%	FNAME in the column vectors X and Y.  Plot X versus Y to see the picture.

%	L. Shure, 12-22-88
%	Modified 1-2-92
%	Copyright (c) 1984-93 by The MathWorks, Inc.

if nargin < 2
    error('X-axis limits must be specified.');
elseif nargin < 5
    subdiv = 20;
    if nargin < 4
        angl = 10;
    end
    if nargin < 3
        npts = 25;
    end
end
if isstr(lims) | isstr(npts) | isstr(angl) | isstr(subdiv)
	error('Optional arguments must be numeric.');
end
iter = 0;
angl = angl*pi/180;
angl = abs(angl);
ang = angl + 1;
xmin = min(lims);
xmax = max(lims);
x = xmin + (0:npts-1)'*(xmax-xmin)/(npts-1);
y = feval(fname,x);
ynan = y(~isnan(y));
yscal = max(max(abs(ynan)));
y = y/yscal;
close_enough = 0;

[m,n] = size(y);
isvector = (n~=1);

% see if need to fill in values for 'smoothness'
while any(ang > angl) & (iter <= subdiv) & (~close_enough)
    [m,n] = size(y);
    yi = diff(y);
    xi = diff(x)*ones(1,n);
    nx = length(xi) - 1;
    radi = sqrt(xi.^2+yi.^2);
    num = xi(1:nx,:).*xi(2:nx+1,:) + yi(1:nx,:).*yi(2:nx+1,:);
    den = radi(1:nx,:).*radi(2:nx+1,:);
    ang = abs(acos(num./den));
    if isvector, ang = max(ang.').'; end
    iter = iter + 1;
    ii = find(ang > angl);    % need to interpolate more
    if isempty(ii)
        break;
    else
	% Choose new points 1/3 to each side of points with large angle.
        xt = [(x(ii) + 2*x(ii+1))/3; (2*x(ii+1) + x(ii+2))/3]; 
        yt = feval(fname,xt)/yscal;
	if all(all((abs([y(ii,:);y(ii+1,:)]-yt)<.01) & (abs([y(ii+1,:);y(ii+2,:)]-yt)<.01)))
		close_enough = 1;
	end
        [x, ind] = sort([x;xt]);
        y = [y;yt];
        y(:) = y(ind,:);
    end
end
y = y*yscal;
if iter >= subdiv
    disp('WARNING: Maximum number of iterations has been reached.');
    if nargout == 0
        disp('         Plot may be inaccurate.');
    end
end
if nargout == 0
    plot(x,y)
else
    x0 = x;
    y0 = y;
end
