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.

fig07_rect_aperture_curved_surface_jitter.m
% fig07_rect_aperture_curved_surface_jitter.m
%--------------------------------------------
% Compute the diffraction of a 2.56 x 2.56 mm quadratic aperture projected
% onto the curved surface of a LINOS lens with diameter d = 20 mm, focal
% length f = 50, and distance L = 2*f = 100 mm.
%
% http://dx.doi.org/10.1016/j.measurement.2013.05.003
% Section 4.2, Fig. 7
%
% wx = wy = 1.28 mm;
% lambda  = 633e-9 nm;
% z       = 100 mm;
%
% 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', 5000, ...					% display a message every x frames
	'save', 50000, ...						% save a snapshot every y frames
    'title','fig07_diffraction', ...		% title of this simulation
    'folder','objects', ...					% folder for this simulation
	'description', 'Fig. 7: rect aperture, projected onto a lens surface, difference w/o jitter', ...	% description text
    'overwrite', false, ...					% overwrite existing data
    'singlePrecisionIsError', true, ...     % throw an error for older/cheaper graphics cards
    'dryRun', false ...
	);

% no light source, include starting conditions in surface generation,
% add surface parameter 'E'

% rectangular aperture, no jitter added
surface1 = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'surface 1 without jitter', ...
	'nGrid', [512 512], ...                 % number imaging grid points [m n]
	'resolution', [5 5] * 1e-6, ...         % imaging 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]
	'E', ones(1, 512 * 512), ...			% evenly illuminated    
	'wavelength', 633e-9 ...				%    with light of this wavelelength
	);

% rectangular aperture, jitter added
surface2 = createStruct('surface', ...
	'subtype', 'plane', ...
	'name', 'surface 2 with jitter', ...
	'nGrid', [512 512], ...                 % number imaging grid points [m n]
	'resolution', [5 5] * 1e-6, ...         % imaging 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]
	'E', ones(1, 512 * 512), ...			% evenly illuminated
	'wavelength', 633e-9, ...				%    with light of this wavelelength
	'xyJitter', [5 5] * 1e-6 ...
	);

% curved lens surfaces:
% LINOS/Qioptiq G063852000, biconvex, focal length 50 mm, NBK-7 material,
D  = 25.4e-3;								% optics diameter,
dm = 7e-3;									% center thickness,
dr = 3.8e-3	;								% edge thickness,
s  = 48.24e-3;								% image distance,
f  = 50e-3;									% focal length,
h  = (dm-dr)/2;								% sagital height,
R  = (4*h^2+s^2)/(8*h);						% lens radius, 
z0 = f + s + dm/2;							% distance to lens center plane at dm/2,

z = @(x,y,P1,P2) -P1(1) + P1(2).*(1-sqrt(1-(x.^2+y.^2)./P1(2).^2));
% z  = -dm/2 + R-sqrt(R^2+x^2+y^2) .......... curved surface height over (x,y)
% P1 ... 'functionParameters', 1xN list of parameters
% P2 ... 'scaledFunctionParameters', 2xN list of scaled parameters, line 1
%			is scaled by 'nGrid(1)', and line 2 by 'nGrid(2)', respectively

r = 10e-3;									% radius of a circular aperture

% lens1 ... curved screen for surface 1 (no jitter)
lens1 = createStruct('surface', ...
	'subtype', 'function', ...
	'name', 'lens front surface', ...
	'size', [20 20] * 1e-3, ...				% size of illuminated area [x y]
	'resolution', [10 10] * 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 z0]', ...			% translation of the plane with respect to the origin [x y z]				
	'function', z, ...						% curved surface height over (x,y)
	'functionParameters', [dm/2 R], ...		% P1(1) = max.thickness/2; P1(2) = R = (4h^2+s^2)/(8h) circular segment
	'apertureFunction',@(x,y) (x.^2+y.^2-r.^2 < 0) & (x < 0), ...	% circular aperture, upper half -x, logical array
	'refractiveIndex', 1.51508 ...			% http://refractiveindex.info/?group=SCHOTT&material=N-BK7; lambda=633nm
	);

% lens2 ... curved screen for surface 2 (with jitter)
lens2 = createStruct('surface', ...
	'subtype', 'function', ...
	'name', 'lens front surface', ...
	'size', [20 20] * 1e-3, ...				% size of illuminated area [x y]
	'resolution', [10 10] * 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 z0]', ...			% translation of the plane with respect to the origin [x y z]				
	'function', z, ...						% curved surface height over (x,y)
	'functionParameters', [dm/2 R], ...		% P1(1) = max.thickness/2; P1(2) = R = (4h^2+s^2)/(8h) circular segment
	'apertureFunction',@(x,y) (x.^2+y.^2-r.^2 < 0) & (x >= 0), ...	% circular aperture, lower half +x, logical array
	'refractiveIndex', 1.51508 ...			% http://refractiveindex.info/?group=SCHOTT&material=N-BK7; lambda=633nm
	);


%% calculate path of light
[lens1, options] = plane2Plane(surface1, lens1, options);
save('lens1.mat','lens1');
[lens2, options] = plane2Plane(surface2, lens2, options);
save('lens2.mat','lens2');


%% show results
% reshape to rectangular grids
E1=reshape(lens1.E,lens1.nGrid(1),[]);
E2=reshape(lens2.E,lens2.nGrid(1),[]);

% compute intensities
I1 = E1 .* conj(E1);							% left, no jitter
I2 = E2 .* conj(E2);							% right, with jitter

% combine half images, add a white separating line
I = I2;
I(1:end/2,:) = I1(1:end/2,:);
I(end/2-15:end/2+15,:) = zeros(31,size(I,1));	% image divider

% compute coordinate grids
X = reshape(lens1.grid(1,:), lens1.nGrid(1),[]);
Y = reshape(lens1.grid(2,:), lens1.nGrid(1),[]);
x = X(:,1);
y = Y(1,:);

% plot and annotate half-images side by side,
% x are vertical columns -> rot90 for horizontal display
h = figure('units','centimeters','Position',[10 10 8 8]);
A = imagesc(1e3*x,1e3*y,rot90(log(I)));
axis equal tight;
set(gca,'YDir','normal')
set(get(A,'Parent'),'XLim',[-10 10])
set(get(A,'Parent'),'YLim',[-10 10])
xlabel('$x$ in mm');
ylabel('$y$ in mm');
set(gca,'LooseInset',get(gca,'TightInset'))

% mymlf2pdf(h,'07diffraction')

Contact us