function ezplot(f,varargin)
%EZPLOT Easy to use function plotter.
%   EZPLOT(f) plots the expression f = f(x) over the default
%   domain -2*pi < x < 2*pi.
%
%   EZPLOT(f, [a,b]) plots f = f(x) over a < x < b
%
%   For implicitly defined functions, f = f(x,y)
%   EZPLOT(f) plots f(x,y) = 0 over the default domain
%   -2*pi < x < 2*pi and -2*pi < x < 2*pi
%   EZPLOT(f, [xmin,xmax,ymin,ymax]) plots f(x,y) = 0 over
%   xmin < x < xmax and  ymin < y < ymax.
%   EZPLOT(f, [a,b]) plots f(x,y) = 0 over a < x < b and a < y < b.
%   If f is a function of the variables u and v (rather than x and y), then
%   the domain endpoints a, b, c, and d are sorted alphabetically. Thus,
%   EZPLOT('u^2 - v^2 - 1',[-3,2,-2,3]) plots u^2 - v^2 - 1 = 0 over
%   -3 < u < 2, -2 < v < 3.
%
%   EZPLOT(x,y) plots the parametrically defined planar curve x = x(t)
%   and y = y(t) over the default domain 0 < t < 2*pi.
%   EZPLOT(x,y, [tmin,tmax]) plots x = x(t) and y = y(t) over
%   tmin < t < tmax.
%
%   EZPLOT(f, [a,b], FIG), EZPLOT(f, [xmin,xmax,ymin,ymax], FIG), or
%   EZPLOT(x,y, [tmin,tmax], FIG) plots the given function over the
%   specified domain in the figure window FIG.
%
%   Examples:
%     ezplot('cos(x)')
%     ezplot('cos(x)', [0, pi])
%     ezplot('1/y-log(y)+log(-1+y)+x - 1')
%     ezplot('x^2 - y^2 - 1')
%     ezplot('x^2 + y^2 - 1',[-1.25,1.25]); axis equal
%     ezplot('x^3 + y^3 - 5*x*y + 1/5',[-3,3])
%     ezplot('x^3 + 2*x^2 - 3*x + 5 - y^2')
%     ezplot('sin(t)','cos(t)')
%     ezplot('sin(3*t)*cos(t)','sin(3*t)*sin(t)',[0,pi])
%     ezplot('t*cos(t)','t*sin(t)',[0,4*pi])

%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.29 $  $Date: 1998/09/09 20:26:01 $

if ~ischar(f) & ~isa(f,'inline')
  error('Input must be a string expression, function name, or INLINE object.');
end
% Place f into "function" form (inline).
if nargin > 1 & (ischar(varargin{1}) | isa(varargin{1},'inline'))
   % Determine whether the two input functions have the same
   % independent variable.  That is, in the case of ezplot(x,y),
   % check that x = x(t) and y = y(t).  If not (x = x(p) and
   % y = y(q)), reject the plot.
   varx = symvar(f); vary = symvar(varargin{1});
   % This is the case of ezplot('2','f(q)') or ezplot('f(p)','3').
   if isempty(varx) | isempty(vary)
      T = union(varx,vary);
      if isempty(T), T = {'x'}; end
      [f,vars] = getfunc(f,T{1});
   elseif isequal(varx,vary)
      [f,vars] = getfunc(f,varx{1});
   else
      [f,vars] = getfunc(f,'t');
   end
else
   [f,vars] = getfunc(f,'x');
end
nvars = length(vars);
ninputs = length(varargin);

switch ninputs
case 0
   % ezplot(f) or ezplot(f(x,y) = 0) -> use default setting.
case 1
   % ezplot(x,y)
   if isa(varargin{1},'char') | length(varargin{1}) == 1
      f = {f, getfunc(varargin{1},vars{:})};
   % ezplot(f(x,y) = 0,[xmin,xmax,ymin,ymax])
   elseif length(varargin{1}) == 4 & nvars == 2
      V = varargin;
      varargin{1} = [V{1}(1),V{1}(2)];
      varargin{2} = [V{1}(3),V{1}(4)];
   end
   % ezplot(f,[xmin,ymin]) covered in the default setting.
case 2
   % ezplot(x,y,domain)
   if (isa(varargin{1},'char') & length(varargin{2}) > 1)
      f = {f,getfunc(varargin{1},vars{:})};
   end
   % ezplot(f,xmin,ymin), ezplot(f,[xmin,xmax],fig), and 
   % ezplot(f(x,y) = 0,[xmin,xmax],[ymin,ymax]) are covered
   % in the default setting.
otherwise
   f = {f,getfunc(varargin{1},'x')};
end

if length(f) == 1
   switch nvars
      case 1
         % Account for variables of [char] length > 1
         if ~isa(f,'inline')
            F = vectorize(inline(char(f),vars{1}));
         else
            F = f;
         end
         ezplot1(F,varargin{:});
         flabel = strrep(char(f),'.','');
         title(texlabel(flabel),'interpreter','tex');
         if ninputs>0 & isa(varargin{1},'double') & length(varargin{1}) == 4
            axis(varargin{1});
         elseif ninputs > 1 & isa(varargin{2},'double') & ...
                length(varargin{2}) == 4
            axis(varargin{2});
         end
      case 2
         ezimplicit(f,varargin{:});
      otherwise
         error([char(f),' cannot be plotted in the xy-plane.']);
      end
elseif length(f) == 2
   ezparam(f{1},f{2},varargin{2:end});
else
   error([f ' cannot be plotted in the xy-plane.']);
end

%------------------------

function ezimplicit(f,varargin)
% EZIMPLICIT Plot of an implicit function in 2-D.
%    EZIMPLICIT(f) plots the string expression f that defines
%    an implicit function f(x,y) = 0 for x0 < x < x1
%    and y0 < y < y1, whose default values are x0 = -2*pi = y0
%    and x1 = 2*pi = y1.
%
%   EZIMPLICIT(f,[x0,x1]) plots the implicit function f(x,y) = 0
%   for x0 < x < x1, x0 < y < x1.
%
%   EZIMPLICIT(f,[x0,x1],[y0,y1]) plots the implicit function
%   f(x,y) = 0 for x0 < x < x1, y0 < y < y1.
%   In the case that f is not a function of x and y 
%   (rather, say u and v), then the domain endpoints [u0,u1] 
%   [v0,v1] are given alphabetically.
%
%   Examples:
%      ezimplicit('x^2 - y^2 - 1',[-4,4])
%      ezimplicit('x^3 + y^3 - 5*x*y + 1/5',[-3,3])
%      ezimplicit('1/y-log(y)+log(-1+y)+x - 1')
%      ezimplicit('x^2 + y^2 = 1',[-1,1])
%      ezimplicit('x^3 - x*y + y^3 - 1',[-2,2])
%      ezimplicit('x^2 + y^2 - (x*y)^2')
%      ezimplicit('x^2 + y^2 = 1',[-1,1])
%      ezimplicit('x^4 + y^4 - 1',[-1,1])
%      ezimplicit('x + tan(x*y)',[-2,2],[-6,6])

% Find the 2 variables that define the implicit function
if isa(f,'cell')
   vars = symvar(char(f));
else
   vars = symvar(f);
end
% Turn f back into a string of chars.
f = char(f);
% If f is a string equation f(x,y) = g(x,y), change
% the equal sign '=' to a minus sign '-'
f = strrep(f,'=','-');
% Choose the number of points in the plot
npts = 250;

fig = [];
switch length(vars)
case 0
   x = 'x'; y = 'y';
case 1
   y = 'y';
case 2
   x = vars{1}; y = vars{2};
otherwise
   % If there are more than 2 variables, send an error message
   W = {vars{1},vars{2}};
   error(['ezplot requires numeric values for ' setdiff(vars,W)]);
end
% Define the computational space
switch nargin
case 1
   X = linspace(-2*pi,2*pi,npts);
   Y = X;
case 2
   if length(varargin{1}) == 1
      fig = varargin{1};
      X = linspace(-2*pi,2*pi,npts); Y = X;
   else 
      X = linspace(varargin{1}(1),varargin{1}(2),npts);
      Y = X;
   end
case 3
   if length(varargin{1}) == 1
      fig = varargin{1};
      X = linspace(varargin{1}(1),varargin{1}(2),npts);
      Y = X;
   elseif length(varargin{2}) == 1 & length(varargin{1}) == 2
      fig = varargin{2};
      X = linspace(varargin{1}(1),varargin{1}(2),npts);
      Y = X;
   elseif length(varargin{2}) == 1 & length(varargin{1}) == 4
      fig = varargin{2};
      X = linspace(varargin{1}(1),varargin{1}(2),npts);
      Y = linspace(varargin{1}(3),varargin{1}(4),npts);
   else
      X = linspace(varargin{1}(1),varargin{1}(2),npts);
      Y = linspace(varargin{2}(1),varargin{2}(2),npts);
   end
end

[X,Y] = meshgrid(X,Y);
u =  eval((idrep(idrep(f,x,'X'),y,'Y')));

% Determine u scale so that "most" of the u values
% are in range, but singularities are off scale.

u = real(u);
uu = sort(u(isfinite(u)));
N = length(uu);
if N > 16
   del = uu(fix(15*N/16)) - uu(fix(N/16));
   umin = max(uu(1)-del/16,uu(fix(N/16))-del);
   umax = min(uu(N)+del/16,uu(fix(15*N/16))+del);
elseif N > 0
   umin = uu(1);
   umax = uu(N);
else
   umin = 0;
   umax = 0;
end
if umin == umax, umin = umin-1; umax = umax+1; end

% Eliminate vertical lines at discontinuities.

ud = (0.5)*(umax - umin); umean = (umax + umin)/2;
[nr,nc] = size(u);
% First, search along the rows . . .
for j = 1:nr
   k = 2:nc;
   kc = find( abs(u(j,k) - u(j,k-1)) > ud );
   ki = find( max( abs(u(j,k(kc)) - umean), abs(u(j,k(kc)-1) - umean) ) );
   if any(ki), u(j,k(kc(ki))) = NaN; end
end
% . . . then search along the columns.
for j = 1:nc
   k = 2:nr;
   kr = find( abs(u(k,j) - u(k-1,j)) > ud );
   kj = find( max( abs(u(k(kr),j) - umean), abs(u(k(kr)-1,j) - umean) ) );
   if any(kj), u(k(kr(kj)),j) = NaN; end
end

if isempty(fig), fig = gcf; end
figure(fig);
contour(X(1,:),Y(:,1),u,[0,0],'-')
xlabel(texlabel(x)); ylabel(texlabel(y));
flabel = strrep(char(f),'.','');
title([texlabel(flabel),' = 0']);

%----------------------------------

function ezparam(varargin)
% EZPARAM Easy to use 2-d parametric curve plotter.
%   EZPARAM(x,y) plots the planar curves r(t) = (x(t),y(t)).
%   The default domain in t [0,2*pi].
%
%   EZPARAM(x,y,[tmin,tmax]) plots r(t) = (x(t),y(t)) for
%   tmin < t < tmax.
%
%   Examples
%       ezparam('sin(t)','cos(t)')
%       ezparam('sin(3*t)*cos(t)','sin(3*t)*sin(t)',[0,pi])
%       ezparam('2*cos(t)','3*sin(t)')
%       ezparam('csc(t)','cot(t)',[0,pi])
%       ezparam('t*cos(t)','t*sin(t)',[0,4*pi])
%       ezparam('t^2,sqrt(1 + t^4)')
%       ezparam('t^2','sqrt(1 + t^4)',[0,1])
%       ezparam('cosh(t)','sinh(t)',[-4,4])

fig = [];
% If x or y are M-files, build them as inline functions of t.
if exist(char(varargin{1})) == 2
   x = vectorize(inline([char(varargin{1}) '(t)'],'t'));
elseif isa(varargin{1},'inline')
   x = varargin{1};
else
   x = char(varargin{1});
end
if exist(char(varargin{2})) == 2
   y = vectorize(inline([char(varargin{2}) '(t)'],'t'));
elseif isa(varargin{2},'inline')
   y = varargin{2};
else
   y = char(varargin{2});
end

% Determine the variables in that define the parametrization
varsx = symvar(x);
varsy = symvar(y);
% Place commmon variables in a set.
vars = union(varsx,varsy);
% Determine the number of variables in r = (x,y)
N = length(vars);

Npts = 300;

% Determine the domains in t:
switch nargin
case 2
   T = linspace(0,2*pi,Npts);
case 3
   if length(varargin{3}) == 1
      fig = varargin{3};
      T = linspace(0,2*pi,Npts);
   else
      T = linspace(varargin{3}(1),varargin{3}(2),Npts);
   end
case 4
   if length(varargin{4}) == 1
      fig = varargin{4};
      T = linspace(varargin{3}(1),varargin{3}(2),Npts);
   elseif length(varargin{3}) == 1
      fig = varargin{3};
      T = linspace(varargin{4}(1),varargin{4}(2),Npts);
   else
      fig = gcf;
      T = linspace(varargin{3},varargin{4},Npts);
   end
end

switch N
   case 1 % planar curve
      t = vars{1};
      if isempty(t), t = 't'; end
   %if length(varsx) == 1 & length(varsy) == 1
      X =  eval((idrep(char(x),t,'T')));
      Y =  eval((idrep(char(y),t,'T')));
      if isempty(fig), fig = gcf; end
      figure(fig);
      plot(X,Y);
      xlabel('x'); ylabel('y');
      axis equal;
      xlabel = strrep(char(x),'.',''); ylabel = strrep(char(y),'.','');
      title(['x = ' texlabel(xlabel), ', y = ' texlabel(ylabel)]);
   otherwise
      error('Cannot plot parametrized surfaces.  Try ezsurf.')
end

%-------------------------

function ezplot1(f,xrange,fig)
%EZPLOT1 Easy to use function plotter.
%   EZPLOT1(f) plots a graph of f(x) where f is a string or a
%   symbolic expression representing a mathematical expression
%   involving a single symbolic variable, say 'x'.  
%   The range of the x-axis is approximately  [-2*pi, 2*pi]
%
%   EZPLOT1(f,xmin,xmax) or EZPLOT(f,[xmin,xmax]) uses the specified
%   x-range instead of the default [-2*pi, 2*pi].
%
%   EZPLOT1(f,[xmin xmax],fig) uses the specified figure number 
%   instead of the current figure.
%
%   Examples:
%       ezplot1('erf(x)')
%       ezplot1 erf(x)
%       ezplot1('tan(sin(x))-sin(tan(x))')
%       ezplot1 tan(sin(x))-sin(tan(x))

% Set defaults

if nargin==0, error('Not enough input arguments.'); end

[F,msg] = fcnchk(f,'vectorized'); % Make a vectorized feval-able function.
if ~isempty(msg), error(msg); end

if nargin < 2, xrange = [-2*pi 2*pi]; end
if isstr(xrange), xrange = eval(xrange); end
if nargin < 3, fig = gcf; end
if nargin == 3
   if length(xrange) == 1
      xrange = [xrange fig]; fig = gcf;
   elseif isstr(fig)
      xrange = [xrange eval(fig)]; fig = gcf;
   end
end
warns = warning;
warning('off');

% Sample on initial interval.

units = get(gca,'units');
set(gca,'units','pixels');
% npts = # of pixels in the axis width.
npts = get(gca,'position')*[0;0;1;0];
set(gca,'units',units);
t = (0:npts-1)/(npts-1);
xmin = min(xrange);
xmax = max(xrange);
x = xmin + t*(xmax-xmin);
y = feval(F,x);
if length(y) == 1, % Check for a scalar function
  for i=1:length(x),
    y(i) = feval(F,x(i));
  end
  y = reshape(y,size(x)); 
end
k = find(abs(imag(y)) > 1.e-6*abs(real(y)));
if any(k), x(k) = []; y(k) = []; end
npts = length(y);
if isempty(y) & npts == 0
   warning on;
   warning(['Cannot plot ' char(f) ':  This function has no real values.']);
   return
end
% Reduce to an "interesting" x interval.

if (npts > 1) & (nargin < 2)
   dx = x(2)-x(1);
   dy = diff(y)/dx;
   dy(npts) = dy(npts-1);
   k = find(abs(dy) > .01);
   if isempty(k), k = 1:npts; end
   xmin = x(min(k));
   xmax = x(max(k));
   if xmin < floor(4*xmin)/4 + dx, xmin = floor(4*xmin)/4; end
   if xmax > ceil(4*xmax)/4 - dx, xmax = ceil(4*xmax)/4; end
   x = xmin + t*(xmax-xmin);
   y = feval(F,x);
   if length(y) == 1, y = y(ones(size(x))); end
   k = find(abs(imag(y)) > 1.e-6*abs(real(y)));
   if any(k), y(k) = NaN; end
end

% Determine y scale so that "most" of the y values
% are in range, but singularities are off scale.

y = real(y);
u = sort(y(isfinite(y)));
npts = length(u);
if isempty(u)
   u = repmat(NaN,size(x));
   npts = prod(size(x));
end
ymin = u(1);
ymax = u(npts);
if npts > 4
   del = u(fix(7*npts/8)) - u(fix(npts/8));
   ymin = max(u(1)-del/8,u(fix(npts/8))-del);
   ymax = min(u(npts)+del/8,u(fix(7*npts/8))+del);
end
 
% Eliminate vertical lines at discontinuities.
 
k = 2:length(y);
k = find( ((y(k) > ymax/2) & (y(k-1) < ymin/2)) | ...
          ((y(k) < ymin/2) & (y(k-1) > ymax/2)) );
if any(k), y(k) = NaN; end

% Plot the function

figure(fig);
plot(x,y)
if ymax > ymin
   axis([xmin xmax ymin ymax])
else
   axis([xmin xmax get(gca,'ylim')])
end

if isa(F,'inline'), v = argnames(F); else v = 'x'; end
xlabel(texlabel(v));

if isstr(F),
  title([f '(x)'],'Interpreter','none'),
elseif isa(f,'inline')
  title(formula(f),'Interpreter','none')
else
  title(f,'Interpreter','none')
end
warning(warns)

%-------------------------

function [F,V] = getfunc(f,defvar)
%GETFUNC Extracts a vectorized inline-function from the
%   input f.
%   [F,V] = GETFUNC(F,DEFVAR) uses the name DEFVAR for 1 input
%   functions with no explicit variables.

% Case 1:  f is a builtin function or existing M-file.
if (isequal(symvar(f),{f}) & exist(f))
   if exist(f)==5
     N = 1; % Assume all built-in functions take one input like sin(x).
   else
     N = nargin(f);
   end
   switch N
   case 1
      V = {defvar};
      F = vectorize(inline([f '(' defvar ')'],defvar));
   case 2
      V = {'x','y'};
      F = vectorize(inline([f '(x,y)'],'x','y'));
   end
% Case 2:  f is a function string (i.e., f = 'x*cos(x)').
elseif ischar(f)
   V = symvar(f);
   if isempty(V), V = {defvar}; end
   F = vectorize(inline(f,V{:}));
% Case 2:  f is an inline funciton.
elseif isa(f,'inline')
   % If f is already vectorized, do nothing . . .
   if findstr(char(f),'.')
      F = f;
   % . . . otherwise, vectorize f.
   else
      F = vectorize(f);
   end
   V = symvar(f);
else
   error([char(f) ' is an invalid function.']);
   V = {};
end

%-------------------------

function s = idrep(s1,s2,s3)
%IDREP Replaces all occurences of s2 in the string s1 by s3.
%   For example, idrep(x*exp(x*y),x,W) returns W*exp(W*y).

% Take care of trivial cases:
%   s2 > s1
%   s1 is empty
%   s2 is empty
if (length(s2) > length(s1)) | isempty(s2) | isempty(s1),
  s=s1;
  return;
end;

% Place non-identifiers at the beginning and end of the string.
s=[]; s1 = [ '#' s1 '#' ];
% Find all occurrences of s2 in s1
s2len=length(s2);
s2pos=findstr(s1,s2);

% Build resulting string
s1pos=1;
for i=1:length(s2pos),
  if ~isalphanum(s1(s2pos(i)-1)) & ~isalphanum(s1(s2pos(i)+s2len))
     s=[s,s1(s1pos:s2pos(i)-1),s3];
     s1pos=s2pos(i)+s2len;
  end
end;
s=[s,s1(s1pos:length(s1))];
% Remove the #'s at the beginning and end of s.
s(1) = []; s(end) = [];

%-------------------------

function B = isalphanum(S)
%ISALPHANUM True for alpha-numeric characters.
%   For a string S, ISALPHANUM returns 1 for alpha-numeric
%   characters (or underscores) and 0 otherwise.
%   Example:
%      S = 'x*exp(x - y) + cosh(x*z^2)';
%      isalphanum(S)  returns
%      (1,0,1,1,1,0,1,0,0,0,1,0,0,0,0,1,1,1,1,0,1,0,1,0,1,0)

B = isletter(S) | (S >= '0' & S <= '9') | (S == '_');
