function [points,t] = fnplt(f,varargin)
%FNPLT Plot a function.
%
%   FNPLT(F)  plots the function in F on its basic interval.
%
%   FNPLT(F,SYMBOL,INTERV,LINEWIDTH,JUMPS) plots the function F on the specified
%   INTERV = [a,b] (default is the basic interval), using the specified
%   plotting SYMBOL (default is '-'), and the specified 
%   LINEWIDTH (default is 1), and using NaNs in order to show any jumps as
%   actual jumps in case JUMPS is a string beginning with 'j'.
%
%   The four optional arguments may appear in any order, with INTERV
%   the one of size [1 2], SYMBOL and JUMPS strings, and LINEWIDTH the
%   scalar.
%
%   If the function in F is 2-vector-valued, the planar curve is
%   plotted.  If the function in F is d-vector-valued with d>2, the
%   space curve given by the first three components of F is plotted.
%
%   If the function is multivariate, it is plotted as a bivariate function
%   at the midpoint of its basic intervals in additional variables, if any.
%
%   POINTS = FNPLT(F,...)   does not plot, but returns instead the sequence 
%   of 2D-points or 3D-points it would have plotted.
%
%   [POINTS,T] = FNPLT(F,...)  also returns, for a vector-valued F, the 
%   corresponding vector T of parameter values.
%
%   For example, fnplt(f,'r',3,[1 3]) plots the graph of the function in F,
%   restricted to the interval [1 .. 3], in red with linewidth 3 (assuming
%   that F describes a univariate function).

%   Carl de Boor 25 feb 89
%   cb : December 1, 1990 (add interval specification)
%   cb : 27 aug 94 (replace ppcut; change so it plots
%                   only when there's no output argument.)
%   cb : 03 nov 95 (allow linewidth, unordered optional inputs, and curve plots)
%   cb : 24 jan 96 (v5 compatibility, include breaks)
%   cb : 21 feb 96 (B-form is plotted as zero outside its basic interval)
%   cb : 17 mar 96 (avoid round-off induced jump in B-form plotting)
%   cb : 29aug96 (add BBform, and optional output of parameter values used)
%   cb : 25oct97 (also plot multivariate functions)
%   cb : 05may98 (standardize the help; accept older version of the forms)
%   cb : 27jun98 (enable display of jump discontinuities across breaks)
%   cb : 09jan00 (enable plotting of rBform = nurb and rpform)
%   cb : 13feb00 (use NaNs for proper jumps only when explicitly required)
%   cb : 07may00 (correct error in handling rational splines)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.13 $

symbol=''; interv=[]; linewidth=[]; jumps=0;
for j=2:nargin
   arg = varargin{j-1};
   if isstr(arg)
      if arg(1)=='j', jumps = 1;
      else, symbol = arg;
      end
   else
      [ignore,d]=size(arg);
      if ignore~=1, error(['arg',num2str(j),' is incorrect.']), end
      if d==1
         linewidth=arg;
      else
         interv=arg;
      end
   end
end

if ~isstruct(f), f = fn2fm(f); end

if ~strcmp(f.form([1 2]),'pp')
   givenform = f.form; f = fn2fm(f,'pp'); basicint = ppbrk(f,'interval');
end

if ~isempty(interv), f = ppbrk(f,interv); end

[breaks,coefs,l,k,d] = ppbrk(f);
if iscell(breaks)
   m = length(breaks);
   for i=m:-1:3
      x{i} = (breaks{i}(1)+breaks{i}(end))/2;
   end
   npoints = 51;
   ii = [1]; if m>1, ii = [2 1]; end
   for i=ii
      x{i}= linspace(breaks{i}(1),breaks{i}(end),npoints);
   end
   v = ppual(f,x);
   if exist('basicint') % we converted from B-form to ppform, hence must now
                        % enforce the basic interval for the underlying spline.
      for i=ii
         temp = find(x{i}<basicint{i}(1)|x{i}>basicint{i}(2));
         if d==1
            if ~isempty(temp), v(:,temp,:) = 0; end
            v = permute(v,[2,1]);
         else
            if ~isempty(temp), v(:,:,temp,:) = 0; end
            v = permute(v,[1,3,2]);
         end
      end
   end
else     % we are dealing with a univariate spline
   npoints = 101;
   x = [breaks(2:l) linspace(breaks(1),breaks(l+1),npoints)];
   v = ppual(f,x); 
   if l>1 % make sure of proper treatment at jumps if so required
      if jumps
         tx = breaks(2:l); temp = zeros(d,l-1); temp(:) = NaN; 
      else
         tx = []; temp = zeros(d,0);
      end
      x = [breaks(2:l) tx x]; 
      v = [ppual(f,breaks(2:l),'left') temp v];
   end
   [x,inx] = sort(x); v = v(:,inx);

   if exist('basicint') % we converted from B-form to ppform, hence must now
                        % enforce the basic interval for the underlying spline.
      temp = find(x<basicint(1)|x>basicint(2));
      if ~isempty(temp), v(temp) = zeros(d,length(temp)); end
   end
end

if exist('givenform')&givenform(1)=='r' % we are dealing with a rational spline
                                        % so, need to divide by last component
   d = d-1;
   sizev = size(v); sizev(1) = d;
   v = reshape(v(1:d,:)./repmat(v(d+1,:),d,1),sizev);
end 
if nargout==0
   if isempty(symbol), symbol = '-'; end
   if isempty(linewidth), linewidth=1; end
   if iscell(x)
      switch d
         case 1
            [yy,xx] = meshgrid(x{2},x{1});
            surf(xx,yy,reshape(v,length(x{1}),length(x{2})))
         case 2
            v = squeeze(v); roughp = 1+(npoints-1)/5;
            vv = reshape(cat(1,...
                 permute(v(:,1:5:npoints,:),[3,2,1]),...
                 repmat(NaN,[1,roughp,2]),...
                 permute(v(:,:,1:5:npoints),[2,3,1]),...
                 repmat(NaN,[1,roughp,2])), ...
                       [2*roughp*(npoints+1),2]);
            plot(vv(:,1),vv(:,2))
         case 3
            v = permute(reshape(v,[3,length(x{1}),length(x{2})]),[2 3 1]);
            surf(v(:,:,1),v(:,:,2),v(:,:,3))
      otherwise
      end
   else
      switch d
         case 1, plot(x,v,symbol,'linew',linewidth)
         case 2, plot(v(1,:),v(2,:),symbol,'linew',linewidth)
         otherwise, plot3(v(1,:),v(2,:),v(3,:),symbol,'linew',linewidth)
      end
   end
else
   if iscell(x)
      switch d
      case 1
         [yy,xx] = meshgrid(x{2},x{1});
         points = {xx,yy,reshape(v,length(x{1}),length(x{2}))};
      case 2
         [yy,xx] = meshgrid(x{2},x{1});
         points = {xx,yy,reshape(v,[2,length(x{1}),length(x{2})])};
      case 3
         points = {squeeze(v(1,:)),squeeze(v(2,:)),squeeze(v(3,:))};
         t = {x{1:2}}
      otherwise
      end
   else
      if d==1, points = [x;v];
      else, t = x; points = v([1:min([d,3])],:); end
   end
end
