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.

getGridNormals(X, Y, Z)
function gridNormals = getGridNormals(X, Y, Z)
% getGridNormals - compute the surface normals in every grid point
%
% X ... [M N] matrix of x-coordinates
% Y ... [M N] matrix of y-coordinates
% Z ... [M N] matrix of z-coordinates

% recover meshgrids
X = myPadArray(X);
Y = myPadArray(Y);
Z = myPadArray(Z);

% differences of neighbouring points
dxx = diff(X);          % differences in x-coordinates in x-direction
dxy = diff(X.').';		% differences in x-coordinates in y-direction
dyx = diff(Y);
dyy = diff(Y.').';
dzx = diff(Z);
dzy = diff(Z.').';

	% clear and save memory
	clear X Y Z;

% differences of next two neighbours,
% cut 2 elements along the longer matrix-dimension, i.e. boundary nodes do
% not have surface normals, hence the padding
dxx2 = dxx(1:end-1,2:end-1) + dxx(2:end,2:end-1);
dxy2 = dxy(2:end-1,1:end-1) + dxy(2:end-1,2:end);
clear dxx dxy;

dyx2 = dyx(1:end-1,2:end-1) + dyx(2:end,2:end-1);
dyy2 = dyy(2:end-1,1:end-1) + dyy(2:end-1,2:end);
clear dyx dyy;

dzx2 = dzx(1:end-1,2:end-1) + dzx(2:end,2:end-1);
dzy2 = dzy(2:end-1,1:end-1) + dzy(2:end-1,2:end);
clear dzx dzy;

% build vectors of differences between x=left/right y=top/down neighbours
vx = zeros(3, numel(dxx2));
vx(1,:) = reshape(dxx2, 1, []);
vx(2,:) = reshape(dyx2, 1, []);
vx(3,:) = reshape(dzx2, 1, []);

vy = zeros(3, numel(dxx2));
vy(1,:) = reshape(dxy2, 1, []);
vy(2,:) = reshape(dyy2, 1, []);
vy(3,:) = reshape(dzy2, 1, []);

	% clear and save memory
	clear dxx2 dxy2 dyz2 dyy2 dzx2 dzy2;

% cross product yields surface normals on the center point
N = cross(vx, vy, 1);

	% clear and save memory
	clear vx vy;

% finally normalize them
absN = sqrt(sum(N.^2));
gridNormals(1,:) = N(1,:) ./ absN;
gridNormals(2,:) = N(2,:) ./ absN;
gridNormals(3,:) = N(3,:) ./ absN;

Contact us