field_propagation_Gerchberg_Saxton_NIST_0503.rar_GS algorithm_Ge

clc;clear; T = '645px-NIST_logo.svg.png'; z=4e-2; %propagation distance, meters l=364e-9; %wavelength, meters k=2*pi/l; %wavenumber k, m^-1 n=2^12; %number of grid points x = 6e-3; %width of grid in real space, m source_w = 1500e-6; %width of gaussian source W = 5e-3;%width of logo in target plane, m yoffset = -1.5e-3; %vertical shift (up for positive value) of the logo w0 = 270e-6; %width of metasurface, m g = 0.9; %diffraction efficiency, 0 <= g <= 1 cutoff = 3.5e5; %spatial frequency cutoff, 1/m x0 = linspace(-x/2,x/2,n); [X,Y]=meshgrid(x0,x0); df=1/x; f=n*df; f0=linspace(-f/2,f/2,n); [FX,FY]=meshgrid(f0,f0); rp = @(x,y,x0,y0)rect(x/x0).*(rect(y/y0)); %change logo to the correct physical size in space target = importdata(T); target = double(target.alpha); s = size(target); x1 = linspace(-W/2,W/2,s(2)); y1 = linspace(-W/2*s(1)/s(2)-yoffset,W/2*s(1)/s(2)-yoffset,s(1)); [X1,Y1]=meshgrid(x1,y1); target = interp2(X1,Y1,target,X,Y,'nearest',0); %a gaussian beam source_B = sqrt(exp(-(X.^2+Y.^2)/2/(source_w)^2)); %a square beam with the same size as the metasurface square source_A = rp(X,Y,w0,w0); s_filter = rp(FX,FY,cutoff,cutoff); %fourier transform the logo in the target plane C0 = jfft2(target); %propagate the fields back to the mask plane D0 = C0.*exp(2i*pi*-z*sqrt(1/l^2-FX.^2-FY.^2)).*circ(sqrt((l*FX).^2+(l*FY).^2)); D0(isnan(D0))=0; %inverse transform to receive the required amplitude and phase at the mask %plane to get a perfect hologram A = jifft2(D0); for j = 1:5 j %replace the amplitude in the mask plane with a square beam B = source_A .* ((1-g)+g*A./abs(A)); %propagate to the target plane B2 = jfft2(B); B3 = B2.*exp(2i*pi*z*sqrt(1/l^2-FX.^2-FY.^2)).*circ(sqrt((l*FX).^2+(l*FY).^2)); C = jifft2(B3); %replace the amplitude in the target plane with the target amplitude D = target .* C ./ abs(C); %propagate back to the mask plane D2 = jfft2(D); D3 = D2.*exp(2i*pi*-z*sqrt(1/l^2-FX.^2-FY.^2)).*circ(sqrt((l*FX).^2+(l*FY).^2)); D3(isnan(D3))=0; D3 = D3.* s_filter; A = jifft2(D3); end Q = A ./ abs(A); mask = Q(round(n/2*(1-w0/x)):round(n/2*(1+w0/x)),round(n/2*(1-w0/x)):round(n/2*(1+w0/x))); %cut the source phase down to a square, apply gaussian beam input actual = source_B.*((1-rp(X,Y,w0,w0))+(rp(X,Y,w0,w0).*((1-g)+g*A./abs(A)))); B12 = jfft2(actual); B13 = B12.*exp(2i*pi*z*sqrt(1/l^2-FX.^2-FY.^2)); C1 = jifft2(B13); jshow(C.^2,1,1);title('with square illumination'); jshow(C1.^2,2,1);title('with gaussian illumination'); %hologram phase p=angle(mask)+pi; %adding