Plot of wave interference from interferometry.

Hi All,

I am looking for a way to simulate an interference image from a Twyman-Green interferometer. In the interferometer there is a flat mirror and a slightly curved mirror. Therefore the interference fringes are tighter together the farther from the center you get. I have yet not been able to replicate that decrease in interference period. I think I would need a for loop that iterates through the phase differences caused by the path difference. Yet I am unsure where to start with that. Or possibly a way to mathematically express a spherical wavefront and have that interfere with the plane wavefront? Any suggestions? Here is what I have so far.

[X,Y] = meshgrid(x,y); % create rectangullar mesh
x=-30:0.5:100; %axis
y=-30:0.5:30; %axis
phi = 4; %phase
k = 1; % wave vector
R = sqrt(X.^2+Y.^2); %rasius
f1 = sin(k*R - phi)
f2 = sin(k + phi);
Z = f1 + f2;
surf(X,Y,Z);
axis equal;

Thanks for anyone's help.

Answer 1

[shukry]

clc;

clear;

close;

%% Input parameters

% screen properties

D1 = 2e-1; % diameter of the source aperture [m]

D2 = 2e-1; % diameter of the observation aperture [m]

N = 512; % sampling numbers

deltan = D2/N; % observation-plane spacing

delta1 = D1/N; % source-plane spacing

% coordinates

[x1, y1] = meshgrid((-N/2 : N/2-1) * delta1); % coordinates

[theta1, r1] = cart2pol(x1, y1);

wvl = 0.532e-6; % optical wavelength [m] % beam parameters

k = 2*pi / wvl; % optical wavenumber [rad/m]

% Beam parameter

w0 = 0.001; % beam width (m)

l = 3; % OAM modes

p = 2;

nscr = 100; % number of phase screens

alpha = 11/3; % power specturm index

L0 = 5; % outer scale [m]

l0 = 0.005; % inner scale [m]

% Spherical wave parameters

R = .0001; % radius of curvature [m] (positive for diverging, negative for converging)

% input fields (x/y polarizations)

normFactor = sqrt(2*factorial(p) / (pi * factorial(abs(l) + p)))/w0;

u0 = normFactor .* exp(-1.*r1.^2./w0.^2) .* (sqrt(2).*r1./w0).^(l) .* Lagureer(p,l, 2.*r1.^2./w0.^2) .* exp(1i*l*theta1);

u2=exp(-1i * k*r1^2/R);

% Add spherical wave component

%spherical_phase = exp(-1i * k * r1.^2 / (2 * R));

u0 = u0+u2;

Cn2 = 1e-15;

dz = 10; % distance interval it should be greater than L0

Tz = dz * nscr; % propagation distance [m]

r0 = (0.423 .* k.^2 .* Cn2 * Tz)^(-3/5); % fried parameter

n = nscr; % planes for masks

% weighting for r0 across screens (same as your script)

sig = 0.01;

r0t1 = r0 * nscr;

r0t2 = r0 * nscr * nscr;

sigt = sig * nscr;

an = zeros(1, nscr);

parfor ii = 1:nscr

an(ii) = 1/(exp(-((ii-1).^2)/(sigt^2))*r0t1 + 1) * r0t1;

end

r0scrn = (r0t2/sum(an)) * an;

d1 = 10e-3;

%% simulate vacuum propagation free space without turbulence propagation Vector

z = (1 : n-1) * Tz / (n-1);

sg = exp(-(x1/(0.47*N*d1)).^16) .* exp(-(y1/(0.47*N*d1)).^16);

t = repmat(sg, [1 1 n]);

[x1, y1, Uvac] = ang_spec_multi_prop(u0, wvl, delta1, deltan, z, t);

% Intensity plot

subplot(2,2,1);

surf(x1, y1, abs(Uvac).^2);

view (2); shading interp; colormap("hot"); grid off; axis off;

title('Vacuum - Intensity');

set(gca,'LooseInset',get(gca,'TightInset'));

% Phase plot

%subplot(2,2,2);

%surf(x1, y1, angle(Uvac));

%view (2); shading interp; colormap("hot"); grid off; axis off;

%title('Vacuum - Phase');

%set(gca,'LooseInset',get(gca,'TightInset'));

%% -----------------------Step three: turbulence propagation-------------------------------------

zt = [0 z]; % propagation plane locations

ad = zt / zt(n);

delta = (1-ad) * delta1 + ad * deltan;

% loading random phase screen

phz = zeros(N, N, n); % initialize array for phase screens

sg_repeated = repmat(sg, [1 1 n]); % Repeat sg to match the dimensions

% loop over screens

parfor idxscr = 1:n

[phz_lo, phz_hi] = ft_sh_phase_screen(r0scrn(idxscr), N, delta(idxscr), alpha, L0, l0);

phz(:,:,idxscr) = phz_lo + phz_hi;

end

% ------------------------- propagate ---------------------------

[x1, y1, Uout] = ang_spec_multi_prop(u0, wvl, delta1, deltan, z, sg_repeated .* exp(1i * phz));

% Intensity plot

subplot(2,2,3);

surf(x1, y1, abs(Uout).^2);

view (2); shading interp; colormap("hot"); grid off; axis off;

title('Turbulence - Intensity');

set(gca,'LooseInset',get(gca,'TightInset'));

% Phase plot

subplot(2,2,4);

surf(x1, y1, angle(Uout));

view (2); shading interp; colormap("hot"); grid off; axis off;

title('Turbulence - Phase');

set(gca,'LooseInset',get(gca,'TightInset'));

%% -----------------------------------Basic Functions-------------------------------------

function [xn, yn, Uout] = ang_spec_multi_prop(Uin, wvl, delta1, deltan, z, t)

N = size(Uin, 1); % number of grid points

[nx, ny] = meshgrid((-N/2 : 1 : N/2 - 1));

k = 2*pi/wvl; % optical wavevector

% super-Gaussian absorbing boundary

nsq = nx.^2 + ny.^2;

w = 0.47*N;

sg = exp(-nsq.^8/w^16); clear('nsq', 'w');

z = [0, z]; % propagation plane locations

n = length(z);

Delta_z = z(2:n) - z(1:n-1); % propagation distances

ad = z / z(n); % grid spacings

delta = (1-ad) * delta1 + ad * deltan;

m = delta(2:n) ./ delta(1:n-1);

x1 = nx * delta(1);

y1 = ny * delta(1);

r1sq = x1.^2 + y1.^2;

Q1 = exp(1i*k/2*(1-m(1))/Delta_z(1)*r1sq);

Uin = Uin .* Q1 .* t(:,:,1);

for idx = 1 : n-1

deltaf = 1 / (N*delta(idx)); % spatial frequencies (of i^th plane)

fX = nx * deltaf;

fY = ny * deltaf;

fsq = fX.^2 + fY.^2;

Z = Delta_z(idx); % propagation distance

Q2 = exp(-1i*pi^2*2*Z/m(idx)/k*fsq); % quadratic phase factor

Uin = sg .* t(:,:,idx+1) .* ift2(Q2 .* ft2(Uin / m(idx), delta(idx)), deltaf); % compute the propagated field

end

xn = nx * delta(n); yn = ny * delta(n); % observation-plane coordinates

rnsq = xn.^2 + yn.^2;

Q3 = exp(1i*k/2*(m(n-1)-1)/(m(n-1)*Z)*rnsq);

Uout = Q3 .* Uin;

end

function [phz_lo, phz_hi] = ft_sh_phase_screen(r0, N, delta, alpha, L0, l0)

persistent last_data;

if isempty(last_data)

last_data = struct();

end

D = N*delta;

phz_hi = ft_phase_screen(r0, N, delta, alpha, L0, l0);

[x,y] = meshgrid((-N/2 : N/2-1) * delta);

phz_lo = zeros(size(phz_hi));

for p = 1:3

del_f = 1 / (3^p*D);

fx = (-1 : 1) * del_f;

[fx, fy] = meshgrid(fx);

fm = 2.68934/l0/(2*pi);

f0 = 4*pi/L0/(2*pi);

PSD_phi = 0.023.*r0.^(2-alpha).*(fx.^2 +fy.^2+f0^2).^(-alpha/2).*exp(-((fx.^2 +fy.^2)./fm.^2)).*(1-0.061.*((fx.^2 +fy.^2).^(1/2)./fm)+2.836.*((fx.^2 +fy.^2).^(1/2)./fm).^(3-alpha./2));

PSD_phi(2,2) = 0;

cn = (randn(3) + 1i*randn(3)).*sqrt(PSD_phi)*del_f;

SH = zeros(N);

for ii = 1:9

SH = SH + cn(ii)*exp(1i*2*pi*(fx(ii)*x+fy(ii)*y));

end

phz_lo = phz_lo + SH;

end

phz_lo = real(phz_lo) - mean(real(phz_lo(:)));

end

function phz = ft_phase_screen(r0, N, delta, alpha, L0, l0)

persistent last_data;

if isempty(last_data)

last_data = struct();

end

del_f = 1/(N*delta);

fx = (-N/2 : N/2-1) * del_f;

[fx, fy] = meshgrid(fx);

fm = 2.68934/l0/(2*pi);

f0 = 4*pi/L0/(2*pi);

PSD_phi = 0.023.*r0.^(2-alpha).*(fx.^2 +fy.^2+f0^2).^(-alpha/2).*exp(-((fx.^2 +fy.^2)./fm.^2)).*(1-0.061.*((fx.^2 +fy.^2).^(1/2)./fm)+2.836.*((fx.^2 +fy.^2).^(1/2)./fm).^(3-alpha./2));

PSD_phi(N/2+1,N/2+1) = 0;

cn = (randn(N) + 1i*randn(N)) .* sqrt(PSD_phi)*del_f;

phz = real(ift2(cn, 1));

end

function G = ft2(g, delta)

G = fftshift(fft2(fftshift(g))) * delta^2;

end

function g = ift2(G, delta_f)

N = size(G, 1);

g = ifftshift(ifft2(ifftshift(G))) * (N * delta_f)^2;

end

function y = Lagureer(p, l, x)

fact = @(x)arrayfun(@factorial, x);

n = p;

k = l;

m = 0:n;

a = factorial(n+k)*ones(1,length(m));

b = fact(n-m);

c = fact(k+m);

d = fact(m);

e = (-1).^m;

y = zeros(size(x));

for s = 1:n+1

y = y + a(s) ./ b(s) ./ c(s) ./ d(s) .* e(s) .* x.^m(s);

end

end

[shukry]

i need a code to study interference between vortex beam and sphrical wave