home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / NUMGRID.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.7 KB  |  52 lines
  1. function G = numgrid(R,n)
  2. %NUMGRID Number the grid points in a two dimensional region.
  3. %    G = NUMGRID('R',n) numbers the points on an n-by-n grid in
  4. %    the subregion of -1<=x<=1 and -1<=y<=1 determined by 'R'.
  5. %    SPY(NUMGRID('R',n)) plots the points.
  6. %    DELSQ(NUMGRID('R',n)) generates the 5-point discrete Laplacian.
  7. %    The regions currently available are:
  8. %       'S' - the entire square.
  9. %       'L' - the L-shaped domain made from 3/4 of the entire square.
  10. %       'C' - like the 'L', but with a quarter circle in the 4-th square.
  11. %       'D' - the unit disc.
  12. %       'A' - an annulus.
  13. %       'H' - a heart-shaped cardioid.
  14. %       'B' - the exterior of a "Butterfly".
  15. %       'N' - a nested dissection ordering of the square.
  16. %    To add other regions, edit toolbox/matlab/demos/numgrid.m.
  17. %
  18. %    See also: DELSQ, DELSQSHOW, DELSQDEMO.
  19.  
  20. %    C. Moler, 7-16-91.
  21. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  22.  
  23. if R == 'N'
  24.    G = nested(n);
  25. else
  26.    x = ones(n,1)*(2*(1:n)-n-1)/(n-1);
  27.    y = flipud(x');
  28.    if R == 'S'
  29.       G = (x > -1) & (x < 1) & (y > -1) & (y < 1);
  30.    elseif R == 'L'
  31.       G = (x > -1) & (x < 1) & (y > -1) & (y < 1) & ( (x > 0) | (y > 0));
  32.    elseif R == 'C'
  33.       G = (x > -1) & (x < 1) & (y > -1) & (y < 1) & ((x+1).^2+(y+1).^2 > 1);
  34.    elseif R == 'D'
  35.       G = x.^2 + y.^2 < 1;
  36.    elseif R == 'A'
  37.       G = ( x.^2 + y.^2 < 1) & ( x.^2 + y.^2 > 1/3);
  38.    elseif R == 'H'
  39.       RHO = .75; SIGMA = .75;
  40.       G = (x.^2+y.^2).*(x.^2+y.^2-SIGMA*y) < RHO*x.^2;
  41.    elseif R == 'B'
  42.       t = atan2(y,x);
  43.       r = sqrt(x.^2 + y.^2);
  44.       G = (r >= sin(2*t) + .2*sin(8*t)) & ...
  45.           (x > -1) & (x < 1) & (y > -1) & (y < 1);
  46.    else
  47.       error('Invalid region type.');
  48.    end
  49.    k = find(G);
  50.    G(k) = (1:length(k))';
  51. end
  52.