8

MATLAB を使用して、コンピューターで生成されたホログラムを作成しようとしています。等間隔のメッシュ グリッドを使用して空間グリッドを初期化すると、次の画像が得られました。

ここに画像の説明を入力

このパターンは、中央の領域を除いて、私が必要としているものです。フリンジはシャープだがぼやけている必要があります。メッシュグリッドの問題かもしれません。pol2cart極座標でグリッドを生成し、MATLAB の関数を使用してデカルト座標にマップしようとしました。残念ながら、それもうまくいきません。細かいグリッドを使用することをお勧めします。それもうまくいきません。スパイラル メッシュ グリッドを生成できれば、おそらく問題は解決できると思います。さらに、一般に、スパイラルアームの数は任意である可能性があります。誰かがこれについてヒントをくれますか?

コードを添付しました (私の最終的なプロジェクトはまったく同じではありませんが、同様の問題があります)。

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)
4

2 に答える 2

1

申し訳ありませんが、数値を投稿できません。しかし、これは役立つかもしれません。振幅空間変調器の実験用に書きました...

R=70;           % radius of curvature of fresnel lens (in pixel units)
A=0;             % oblique incidence by linear grating (1=oblique 0=collinear)
B=1;             % expanding by fresnel lens (1=yes 0=no)
L=7;            % topological charge
Lambda=30;       % linear grating fringe spacing (in pixels)
aspect=1/2;      % fraction of fringe period that is white/clear
xsize=1024;       % resolution (xres x yres number data pts calculated)
ysize=768;        % 

% define the X and Y ranges (defined to skip zero)
xvec = linspace(-xsize/2, xsize/2, xsize);     % list of x values
yvec = linspace(-ysize/2, ysize/2, ysize);     % list of y values

% define the meshes - matrices linear in one dimension
[xmesh, ymesh] = meshgrid(xvec, yvec);

% calculate the individual phase components
vortexPh = atan2(ymesh,xmesh);       % the vortex phase
linPh = -2*pi*ymesh;           % a phase of linear grating
radialPh = (xmesh.^2+ymesh.^2);     % a phase of defocus

% combine the phases with appropriate scales (phases are additive)
% the 'pi' at the end causes inversion of the pattern
Ph = L*vortexPh + A*linPh/Lambda + B*radialPh/R^2;

% transmittance function (the real part of exp(I*Ph))
T = cos(Ph);
% the binary version
binT = T > cos(pi*aspect);

% plot the pattern
% imagesc(binT)
imagesc(T)
colormap(gray)
于 2013-07-10T13:57:40.633 に答える