title: spiral meshgrid in matlab
body: <p>I'm trying to produce some computer generated holograms by using MATLAB. I used equally spaced mesh grid to initialize the spatial grid, and I got the following image</p>
  
  <p><img src="https://i.stack.imgur.com/XP9aV.jpg" alt="enter image description here"></p>
  
  <p>This pattern is sort of what I need except the center region. The fringe should be sharp but blurred. I think it might be the problem of the mesh grid. I tried generate a grid in polar coordinates and the map it into Cartesian coordinates by using MATLAB's <code>pol2cart</code> function. Unfortunately, it doesn't work as well. One may suggest that using fine grids. It doesn't work too. I think if I can generate a spiral mesh grid, perhaps the problem is solvable. In addition, the number of the spiral arms could, in general, be arbitrary, could anyone give me a hint on this?</p>
  
  <p>I've attached the code (My final projects are not exactly the same, but it has a similar problem).</p>
  
  <pre><code>clc; clear all; close all;
  %% initialization
  tic
  lambda = 1.55e-6;
  k0 = 2*pi/lambda;
  c0 = 3e8;
  eta0 = 377;
  scale = 0.25e-6;
  NELEMENTS = 1600;
  GoldenRatio = (1+sqrt(5))/2;
  g = 2*pi*(1-1/GoldenRatio);
  
  pntsrc = zeros(NELEMENTS, 3);
  phisrc = zeros(NELEMENTS, 1);
  for idxe = 1:NELEMENTS
    pntsrc(idxe, :) = scale*sqrt(idxe)*[cos(idxe*g), sin(idxe*g), 0];
    phisrc(idxe) = angle(-sin(idxe*g)+1i*cos(idxe*g));
  end
  phisrc = 3*phisrc/2; % 3 arms (topological charge ell=3)
  
  %% post processing
  sigma = 1;
  polfilter = [0, 0, 1i*sigma; 0, 0, -1; -1i*sigma, 1, 0]; % cp filter
  
  xboundl = -100e-6; xboundu = 100e-6;
  yboundl = -100e-6; yboundu = 100e-6;
  xf = linspace(xboundl, xboundu, 100);
  yf = linspace(yboundl, yboundu, 100);
  zf = -400e-6;
  [pntobsx, pntobsy] = meshgrid(xf, yf);
  % how to generate a right mesh grid such that we can generate a decent result?
  pntobs = [pntobsx(:), pntobsy(:), zf*ones(size(pntobsx(:)))];
  % arbitrary mesh may result in "wrong" results
  
  NPNTOBS = size(pntobs, 1);
  nxp = length(xf);
  nyp = length(yf);
  
  %% observation
  Eobs = zeros(NPNTOBS, 3);
  
  matlabpool open local 12
  parfor nobs = 1:NPNTOBS
    rp = pntobs(nobs, :);
    Erad = [0; 0; 0];
    for idx = 1:NELEMENTS
      rs = pntsrc(idx, :);
      p = exp(sigma*1i*2*phisrc(idx))*[1 -sigma*1i 0]/2; % simplified here
      u = rp - rs;
      r = sqrt(u(1)^2+u(2)^2+u(3)^2); %norm(u);
      u = u/r; % unit vector
      ut = [u(2)*p(3)-u(3)*p(2),...
        u(3)*p(1)-u(1)*p(3), ...
        u(1)*p(2)-u(2)*p(1)]; % cross product: u cross p
      Erad = Erad + ... % u cross p cross u, do not use the built-in func
        c0*k0^2/4/pi*exp(1i*k0*r)/r*eta0*...
        [ut(2)*u(3)-ut(3)*u(2);...
        ut(3)*u(1)-ut(1)*u(3); ...
        ut(1)*u(2)-ut(2)*u(1)]; 
    end
    Eobs(nobs, :) = Erad; % filter neglected here
  end
  matlabpool close
  Eobs = Eobs/max(max(sum(abs(Eobs), 2))); % normailized
  
  %% source, gaussian beam
  E0 = 1;
  w0 = 80e-6;
  theta = 0; % may be titled
  RotateX = [1, 0, 0; ...
    0, cosd(theta), -sind(theta); ...
    0, sind(theta), cosd(theta)];
  
  Esrc = zeros(NPNTOBS, 3);
  for nobs = 1:NPNTOBS
    rp = RotateX*[pntobs(nobs, 1:2).'; 0];
    z = rp(3);
    r = sqrt(sum(abs(rp(1:2)).^2));
    zR = pi*w0^2/lambda;
    wz = w0*sqrt(1+z^2/zR^2);
    Rz = z^2+zR^2;
    zetaz = atan(z/zR);
    gaussian = E0*w0/wz*exp(-r^2/wz^2-1i*k0*z-1i*k0*0*r^2/Rz/2+1i*zetaz);% ...
    Esrc(nobs, :) = (polfilter*gaussian*[1; -1i; 0]).'/sqrt(2)/2;
  end
  Esrc = [Esrc(:, 2), Esrc(:, 3), Esrc(:, 1)];
  Esrc = Esrc/max(max(sum(abs(Esrc), 2)));  % normailized
  toc
  
  %% visualization
  fringe = Eobs + Esrc; % I'll have a different formula in my code
  normEsrc = reshape(sum(abs(Esrc).^2, 2), [nyp nxp]);
  normEobs = reshape(sum(abs(Eobs).^2, 2), [nyp nxp]);
  normFringe = reshape(sum(abs(fringe).^2, 2), [nyp nxp]);
  
  close all;
  xf0 = linspace(xboundl, xboundu, 500);
  yf0 = linspace(yboundl, yboundu, 500);
  [xfi, yfi] = meshgrid(xf0, yf0);
  data = interp2(xf, yf, normFringe, xfi, yfi);
  figure; surf(xfi, yfi, data,'edgecolor','none');
  % tri = delaunay(xfi, yfi); trisurf(tri, xfi, yfi, data, 'edgecolor','none');
  xlim([xboundl, xboundu])
  ylim([yboundl, yboundu])
  % colorbar
  view(0,90)
  colormap(hot)
  axis equal
  axis off
  title('fringe thereo. ', ...
    'fontsize', 18)
  </code></pre>