Code covered by the BSD License

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

### Dominik Hofer (view profile)

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

gauss2Plane(lightsource, surface, options)
```function surface = gauss2Plane(lightsource, surface, options)
% gauss2plane - calculates the complex E field vector for every point on a
% grid if it is illuminated by a gaussian beam.
%
% http://dx.doi.org/10.1016/j.measurement.2013.05.003
% Section 3.1, Fig. 2

% map to more convenient variables
A		= surface.grid;					% surface grid, [xi; yi; zi]
As		= surface.resolution;			% surface grid spacing, [dx dy]
An		= surface.normal;				% main surface normal
Ad		= surface.normalDistance;		% minimal distance from the surface to the origin
ir		= lightsource.translations;		% displacement of the laser aperture
in		= lightsource.normal;			% laser beam direction
iLambda = lightsource.wavelength;				% laser wavelength
lx		= lightsource.coordinateSystem(:,1);	% light source x-axis
ly		= lightsource.coordinateSystem(:,2);	% light source y-axis

% find intersection laser beam/surface
% solve Hesse normal form: Anx*x + Any*y + Anz*z - Ad == 0
% where ir + lambda * in = [x y z]
% intersection point = aperture + lambda * propagation vector
lambda = (Ad - An' * ir) / (An' * in);
p = ir + lambda .* in;

% distance grid points - intersection point
d = A - p(:,ones(1,size(A,2)));

% project on aperture plane, get distances, too much data for pure matrix
% component of 'd' in 'in'-direction
dn = sum(d.*in(:,ones(1,size(d,2)))) / (in' * in);

% component of 'd' parallel to aperture, shortest distance to laser beam
da = d - dn(ones(3,1),:).*in(:,ones(1,size(d,2)));

% projection of 'da' on the light source' x- and y-axis
% 'lx' and 'ly' are assumed to be unit vectors
da_x = (da.' * lx).';	% magnitude of (da' * lx) / (lx' * lx) * lx
da_y = (da.' * ly).';

% its norm, length of the vectors
da_abs = sqrt(da(1,:).^2+da(2,:).^2+da(3,:).^2);

% illumination, simple gauss
g = lightsource.function;

% irradiance * pixelarea * gauss * exp(1i * dphi)
surface.E = iMag * (As(1)*As(2)) * g(da_x, da_y, da_abs, ...
lightsource.parameters) .* exp(1i * 2 * pi .* dn ./ iLambda);

% compute vectors of ideal reflection for every point, mirror light
% propagation vector on surface normals
if isempty(surface.mainReflectionVectors)
displayMessage(toc(options.tStart), 2, [' ' surface.name ': main reflection vectors']);
surface.mainReflectionVectors = getReflection(in, surface.gridNormals);
end```