Code covered by the BSD License  

Highlights from
Huygens-Fresnel integral approximation, free-form apertures and rough surfaces

image thumbnail

Huygens-Fresnel integral approximation, free-form apertures and rough surfaces

by

 

Simulate wave propagation through free-form apertures, or off rough surfaces.

fig05_rectangular_aperture.m
% fig05_rectangular_aperture.m
%-----------------------------
% Compute the diffraction of a 1x1 mm quadratic aperture projected onto a
% screen at a distance of 1000 m.
%
% http://dx.doi.org/10.1016/j.measurement.2013.05.003
% Section 4.2, Fig. 5
%
% wx = wy = 0.5 mm;
% lambda  = 658e-9 nm;
% z       = 1000 m;
%
% NEEDS GPUmat (http://sourceforge.net/projects/gpumat/)

% prepare workspace
clear; clc;									% clear up

addpath('gpuSpeckle');						% add function folder to PATH variable
GPUstart();									% start GPUmat

%% define optical components in the path of light

% general options
options = createStruct('options', ...
	'display', 500, ...						% display a message every x frames
	'save', 100000, ...						% save a snapshot every y frames
	'title','fig05_rect_aperture_1', ...	% title, 1st part of folder name
	'description', 'Fig 6: rect aperture (wx = wy = 0.5 mm, z = 1000 m)', ...	% description text
	'singlePrecisionIsError', true, ...     % throw an error for older/cheaper graphics cards
    'overwrite',true ...					% overwrite existing data
	);

% light source
% 'function' ... left out to use standard 2d Gaussian beam
lightsource = createStruct('source', ...
	'subtype', 'planeWave', ...
	'name', 'light source', ...
	'wavelength', 658e-9, ...				% wavelength
	'irradiance', 1, ...                    % power in watts/m^2
	'normalDistance', 1, ...				% normal distance to base coordinate rigin
	'rotations', [0 180 0]*pi/180 ...		% rotations of aperture with respect to the origin [alpha beta gamma]
	);

% reflecting surface/aperture
surface = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'rough surface', ...
	'size', [1 1] * 1e-3, ...				% size of illuminated area [x y]
	'resolution', [.5 .5] * 1e-6, ...			% source grid spacing [dx dy]
	'rotations', [0 0 0]*pi/180, ...		% rotations of plane with respect to the origin [alpha beta gamma]
	'translations', [0 0 0]' ...			% translation of the plane with respect to the origin [x y z]
	);

% observation screen
camera = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'camera', ...
	'nGrid', [64 3], ...					% number imaging grid points [m n]
	'size', [4 4], ...						% size of active area [x y]
	'normalDistance', 1000, ...				% normal distance to base coordinate rigin
	'rotations', [0 180 0]*pi/180 ...		% rotations with respect to the origin [alpha beta gamma]
	);

%% calculate path of light
surface = light2Plane(lightsource, surface, options);
[camera, options] = plane2Plane(surface, camera, options);

%% show results
% prepare empty figure, aspect ratio 4:3
set(0,'defaulttextinterpreter','none');
h=figure('units','centimeters', 'Position',[10 10 10 10/4*3]); % position: [llx lly width height]

% map to more convenient variables
s = camera.resolution;		% spatial resolution [dx dy]
m = camera.nGrid;			% number of points [m n]

% coordinate vectors
x = -s(1) * (m(1)-1)/2 : s(1) : s(1) * (m(1)-1)/2;
y = -s(2) * (m(2)-1)/2 : s(2) : s(2) * (m(2)-1)/2;

% camera E field is in (W/m)/m -> scale by pixel area, compute intensity
C = reshape(camera.E, camera.nGrid(1),[]).*(s(1)*s(2));
C = C.*conj(C);

% Fraunhofer diffraction: (wx = wy = 0.5e-3, z = 1000)
% f = I(x,y) = ((2 wx 2 wy)/(lambda z))^2 sinc(2 wx x /(lambda z))^2
% f1 ... scaled by 1e6 -> unit: \mu W / m^2, nicer tick values
f = @(x) ((2*.5e-3)^2/(658e-9*1000))^2*sinc(2*.5e-3*x/(658e-9*1000)).^2;
f1 = @(x) 1e6.*f(x);

% plot simulated function values and errors with respect to f(x)
[AX,H1,H2] = plotyy(x,1e6.*C(:,floor(end/2+1)),x,1e9.*(C(:,floor(end/2+1))-f(x).'));
grid on; hold on;
% analytical solution f1(x)
H3 = ezplot(f1,[-2,2,0,f1(0)]);

% annotate figure
title('Rectangular Aperture ($2w_x$=$2w_y$=1e-3, z=1000)');
xlabel('Screen Coordinate in m');
ylabel(AX(1),'Irradiance $I$ in $\upmu$W/m$^2$');
H4 = ylabel(AX(2),' $\Delta I$ in nW/m$^2$');
set(H4, 'Position', get(H4,'Position')+[.1 0 0]);
set(AX(1),'YLim',[0, 2.5]);
set(AX(2),'YLim',[-5, 5]);
set(AX(1),'YTick',[0 .5 1 1.5 2 2.5]);
set(AX(2),'YTick',[-5 -3 -1 1 3 5]);
set(H1,'LineStyle','none','Marker','o','Color','red');
set(H1,'LineWidth',2); set(H2,'LineWidth',2); set(H3,'LineWidth',2);
set(H2,'LineStyle','--');
set(gca,'looseinset',get(gca,'tightinset'));

% mymlf2pdf(h,'05rectAperture', 'fontsize', 9, 'packages',{{'','upgreek'}})

Contact us