1 function [xr,yr,zr,xi,yi,zi] = CmplxFacets(R,e,TypeDomain,TypeCut,n,StrFunc)
   2 
   3         //  A function to compute the facets for drawing a complex function
   4         //  on a square or a disk with branch(es) cut(s) on Ox or Oy
   5         //
   6         //  TypeDomain : "Square" or "Disk"
   7         //     TypeCut : "Ox" or "Oy"
   8         //           R : length of half a side of the square or radius of the disk
   9         //           e : thin layer to avoid the branch(es) cut(s)
  10         //           n : a scalar (for Square) or a 2-vector = [ntheta, nr]
  11         //               (for Disk) for discretization
  12         //     StrFunc : the string which names the complex function (this is
  13         //               because primitive don't pass as function argument)
  14         
  15         if TypeDomain == "Square" then
  16                 if TypeCut == "Ox" then
  17                         x1 = linspace(-R, R, n);
  18                         y1 = linspace( e, R, n/2);
  19                 else  // for TypeCut = "Oy" ...
  20                         x1 = linspace( e, R, n/2);
  21                         y1 = linspace(-R, R, n);
  22                 end
  23                 X1 = ones(y1')*x1 ; Y1 = y1'*ones(x1);
  24         else // for TypeDomain = "Disk"
  25                 r = linspace(0,R, n(2));
  26                 if TypeCut == "Ox" then
  27                         theta = linspace(0,%pi,n(1))';
  28                         X1 = cos(theta)*r;
  29                         Y1 = e + sin(theta)*r;
  30                 else // for TypeCut = "Oy"
  31                         theta = linspace(-%pi/2,%pi/2,n(1))';
  32                         X1 = e + cos(theta)*r;
  33                         Y1 = sin(theta)*r;
  34                 end
  35         end
  36         X2 = -X1 ; Y2 = -Y1;
  37         Z1 = evstr(StrFunc+"(X1 + %i*Y1)");
  38         Z2 = evstr(StrFunc+"(X2 + %i*Y2)");
  39         [xr1,yr1,zr1] = nf3d(X1,Y1,real(Z1));
  40         [xr2,yr2,zr2] = nf3d(X2,Y2,real(Z2));
  41         xr = [xr1 xr2]; yr = [yr1 yr2]; zr = [zr1 zr2];
  42         [xi1,yi1,zi1] = nf3d(X1,Y1,imag(Z1));
  43         [xi2,yi2,zi2] = nf3d(X2,Y2,imag(Z2));
  44         xi = [xi1 xi2]; yi = [yi1 yi2]; zi = [zi1 zi2];
  45 endfunction