home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / DATAFUN.DI$ / GRADIENT.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.6 KB  |  59 lines
  1. function [xx,yy] = grad(a,xax,yax)
  2. %GRADIENT Approximate gradient.
  3. %    [PX,PY] = GRADIENT(Z,DX,DY) returns the numerical partial derivatives
  4. %    of matrix Z in matrices PX = dZ/dx and PY = dZ/dy.   DX and DY
  5. %    may be scalars containing the sample spacing in the X and Y
  6. %    directions, or they may be vectors containing all the explicit
  7. %    locations.
  8. %
  9. %    [PX,PY] = GRADIENT(Z) assumes DX = DY = 1.
  10. %
  11. %    If Y is a vector, GRADIENT(Y) and GRADIENT(Y,DX) return the one
  12. %    dimensional numerical derivative dY/dX.
  13. %
  14. %    For example, try
  15. %       [x,y] = meshgrid(-2:.2:2, -2:.2:2);
  16. %       z = x .* exp(-x.^2 - y.^2);
  17. %       [px,py] = gradient(z,.2,.2);
  18. %       contour(z),hold on, quiver(px,py), hold off
  19. %
  20. %    See also DIFF, DEL2, QUIVER, CONTOUR.
  21.  
  22. %    Charles R. Denham, MathWorks 3-20-89
  23. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  24.  
  25. [m,n] = size(a);
  26. if nargin == 1, xax = 1; yax = 1; end
  27. if nargin == 2, yax = xax; end
  28. if length(xax) == 1, xax = xax .* (0:n-1); end
  29. if length(yax) == 1, yax = yax .* (0:m-1); end
  30.  
  31. y = [];
  32. ax = xax(:).';
  33. for i = 1:2
  34.    x = y;
  35.    [m,n] = size(a);
  36.    y = zeros(m, n);
  37.    j = 1:m;
  38.    if n > 1
  39.       d = ax(2) - ax(1);
  40.       y(j, 1) = (a(j, 2) - a(j, 1)) ./ d;     % Left edge.
  41.       d = ax(n) - ax(n-1);
  42.       y(j, n) = (a(j, n) - a(j, n-1)) ./ d;   % Right edge.
  43.    end
  44.    if n > 2
  45.       k = 1:n-2;
  46.       d = ones(m, 1) * (ax(k+2) - ax(k));
  47.       y(j, k+1) = (a(j, k+2) - a(j, k)) ./ d;   % Middle.
  48.    end
  49.    a = a.';
  50.    ax = yax(:).';
  51. end
  52. z = (x + sqrt(-1) .* y.');
  53. if nargout < 2
  54.    xx = z;
  55.  else
  56.    xx = real(z); yy = imag(z);
  57. end
  58.  
  59.