home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / PLOTXY.DI$ / FPLOT.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.8 KB  |  92 lines
  1. function [x0,y0] = fplot(fname,lims,npts,angl,subdiv)
  2. %FPLOT    Plot a function.
  3. %    FPLOT(FNAME,LIMS) plots the function specified by the text variable
  4. %    FNAME between the x-axis limits specified by LIMS = [XMIN XMAX].
  5. %    For example, FPLOT('sin',[0 10]) is a plot of the sine function.
  6. %    FPLOT(FNAME,LIMS,N) plots the function with a minimum of N samples.
  7. %    The default value for N is 25.
  8. %    FPLOT(FNAME,LIMS,N,ANGLE) plots the function with a minimum of N
  9. %    samples and samples more if contiguous segments are at an angle
  10. %    larger than ANGLE (in degrees).  The default for ANGLE is 10 degrees.
  11. %    FPLOT(FNAME,LIMS,N,ANGLE,SUBDIV) plots the function with a minimum
  12. %    of N samples, no more than ANGLE degrees bends in the curve, but no
  13. %    more than SUBDIV iterations to fill in the rapidly changing portions.
  14. %    The default for SUBDIV is 20.
  15. %    [X,Y] = fplot(FNAME,LIMS,...) returns the abscissae and ordinates for
  16. %    FNAME in the column vectors X and Y.  Plot X versus Y to see the picture.
  17.  
  18. %    L. Shure, 12-22-88
  19. %    Modified 1-2-92
  20. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  21.  
  22. if nargin < 2
  23.     error('X-axis limits must be specified.');
  24. elseif nargin < 5
  25.     subdiv = 20;
  26.     if nargin < 4
  27.         angl = 10;
  28.     end
  29.     if nargin < 3
  30.         npts = 25;
  31.     end
  32. end
  33. if isstr(lims) | isstr(npts) | isstr(angl) | isstr(subdiv)
  34.     error('Optional arguments must be numeric.');
  35. end
  36. iter = 0;
  37. angl = angl*pi/180;
  38. angl = abs(angl);
  39. ang = angl + 1;
  40. xmin = min(lims);
  41. xmax = max(lims);
  42. x = xmin + (0:npts-1)'*(xmax-xmin)/(npts-1);
  43. y = feval(fname,x);
  44. ynan = y(~isnan(y));
  45. yscal = max(max(abs(ynan)));
  46. y = y/yscal;
  47. close_enough = 0;
  48.  
  49. [m,n] = size(y);
  50. isvector = (n~=1);
  51.  
  52. % see if need to fill in values for 'smoothness'
  53. while any(ang > angl) & (iter <= subdiv) & (~close_enough)
  54.     [m,n] = size(y);
  55.     yi = diff(y);
  56.     xi = diff(x)*ones(1,n);
  57.     nx = length(xi) - 1;
  58.     radi = sqrt(xi.^2+yi.^2);
  59.     num = xi(1:nx,:).*xi(2:nx+1,:) + yi(1:nx,:).*yi(2:nx+1,:);
  60.     den = radi(1:nx,:).*radi(2:nx+1,:);
  61.     ang = abs(acos(num./den));
  62.     if isvector, ang = max(ang.').'; end
  63.     iter = iter + 1;
  64.     ii = find(ang > angl);    % need to interpolate more
  65.     if isempty(ii)
  66.         break;
  67.     else
  68.     % Choose new points 1/3 to each side of points with large angle.
  69.         xt = [(x(ii) + 2*x(ii+1))/3; (2*x(ii+1) + x(ii+2))/3]; 
  70.         yt = feval(fname,xt)/yscal;
  71.     if all(all((abs([y(ii,:);y(ii+1,:)]-yt)<.01) & (abs([y(ii+1,:);y(ii+2,:)]-yt)<.01)))
  72.         close_enough = 1;
  73.     end
  74.         [x, ind] = sort([x;xt]);
  75.         y = [y;yt];
  76.         y(:) = y(ind,:);
  77.     end
  78. end
  79. y = y*yscal;
  80. if iter >= subdiv
  81.     disp('WARNING: Maximum number of iterations has been reached.');
  82.     if nargout == 0
  83.         disp('         Plot may be inaccurate.');
  84.     end
  85. end
  86. if nargout == 0
  87.     plot(x,y)
  88. else
  89.     x0 = x;
  90.     y0 = y;
  91. end
  92.