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.

huygensFresnel(objA, objB, options, index)
function B = huygensFresnel(objA, objB, options, index)
% huygensFresnel - calculate sum of phasors
%
% objA ... source
% objB ... destination
% options
% index ... first element to process, next element = index+objA.stride


%% define some stats
saves    = 0;                    % number of saves accomplished
err      = 0;                    % number of array indexing errors
dispCnt_ = options.display;
saveCnt_ = options.save;


%% check presence of old 'data' files, resume their computation
fileList = dir([options.dataFolder filesep '*data*.mat']);
Btemp = [];
if length(fileList) >= 1
    displayMessage(2, [num2str(length(fileList)) ' save files found, begin rebuilding old solution']); 
    fileNames = {fileList.name}';
	
	for m = 1:length(fileNames)
		load([options.dataFolder filesep fileNames{m}]);
	end
	
    % rebuild result, clear partial data
    for m = 1:length(fileNames)
		saves = saves + 1;
        eval(['Btemp = [Btemp B' num2str(m) '];']);
        clear(['B' num2str(m)]);
    end
    
    % set first index to process, replace E2
    options.startIndex = length(Btemp) + 1;
    displayMessage(2, 'start to resume old job'); 
end


%% define GPU data
%  considering aperture indices and stride

% source field, only use values within the aperture
if ~isempty(objA.aperture)
	AX = objA.grid(:,objA.aperture);
	A_ = objA.E(:,objA.aperture);	
else
	AX = objA.grid;
	A_ = objA.E;	
end
AXx	= toGpu(AX(1,index:objA.stride:end), options);
AXy	= toGpu(AX(2,index:objA.stride:end), options);
AXz	= toGpu(AX(3,index:objA.stride:end), options);
if ~isempty(objA.E)
	A	= toGpu(A_(index:objA.stride:end), options);
else
	displayMessage(-1, [objA.name ': ERROR missing source E-field']);
end

% use indicatrix main lobe direction or surface normals
if ~isempty(objA.mainReflectionVectors)
	if ~isempty(objA.aperture)
		mainReflectionVectors = objA.mainReflectionVectors(:,objA.aperture);
	else
		mainReflectionVectors = objA.mainReflectionVectors;
	end
	ANx = toGpu(mainReflectionVectors(1,index:objA.stride:end), options);
	ANy = toGpu(mainReflectionVectors(2,index:objA.stride:end), options);
	ANz = toGpu(mainReflectionVectors(3,index:objA.stride:end), options);
else
	if ~isempty(objA.aperture)
		gridNormals = objA.gridNormals(:,objA.aperture);
	else
		gridNormals = objA.gridNormals;
	end
	ANx = toGpu(gridNormals(1,index:objA.stride:end), options);
	ANy = toGpu(gridNormals(2,index:objA.stride:end), options);
	ANz = toGpu(gridNormals(3,index:objA.stride:end), options);
end

	% save memory
	clear 	AX A_ mainReflectionVectors gridNormals;

% destination field
if ~isempty(objB.aperture)
	aperture = objB.aperture;
else
	aperture = ':';
end
BXx	= objB.grid(1,aperture);
BXy	= objB.grid(2,aperture);
BXz	= objB.grid(3,aperture);
B = zeros(1,size(BXx,2));
if ~isempty(Btemp)
	B(1:options.startIndex-1) = Btemp;			% copy over old data
end
B = toGpu(B, options);

if isempty(objA.wavelength)
	displayMessage(-1, [objA.name ': ERROR missing parameter ''wavelength''']);
end

%% compute sum, start looping
displayMessage(3, ['Loop starts: ' num2str(size(B,2)) ' points to go']); 

N		= size(B,2);						% number of points to compute
ANabs	= sqrt(ANx.*ANx+ANy.*ANy+ANz.*ANz);	% length of all normal vectors; .* numerically better than .^2
k		= 2 * pi / objA.wavelength * objA.refractiveIndex; % wave number & optical path

A = A ./ objA.refractiveIndex;				% optical path

% dummy variables
dabs = toGpu(zeros(1,length(ANabs)), options);
cosA = toGpu(zeros(1,length(ANabs)), options);

if options.dryRun == false
	if options.useGpu						% GPU specific code
		% more dummies
		d1 = toGpu(zeros(1,length(ANabs)), options);
		d2 = toGpu(zeros(1,length(ANabs)), options);
		d3 = toGpu(zeros(1,length(ANabs)), options);
		setComplex(d3);
		
		for n = options.startIndex : N

			% path lengths
			dx = BXx(n) - AXx;
			dy = BXy(n) - AXy;
			dz = BXz(n) - AXz;

			% dabs = sqrt(dx.*dx + dy.*dy + dz.*dz);
			GPUtimes(dx,dx,d1);
			GPUtimes(dy,dy,d2);
			GPUplus(d1,d2,d1);
			GPUtimes(dz,dz,d2);
			GPUplus(d1,d2,d1);
			GPUsqrt(d1,dabs);

			% compute cos(angle(r01,n)) for the aperture: cos() = a.b /(|a|*|b|)
			% cosA = (dx.*ANx + dy.*ANy + dz.*ANz) ./ (dabs .* ANabs);
			GPUtimes(dx,ANx,d1);
			GPUtimes(dy,ANy,d2);
			GPUplus(d1,d2,d1);
			GPUtimes(dz,ANz,d2);
			GPUplus(d1,d2,d1);
			GPUtimes(dabs,ANabs,d2);
			GPUrdivide(d1,d2,cosA);

			% and the obliquity factor 1/2*(1+cos(theta)) for Kirchhoff theory
			GPUplus(1,cosA,d1);
			GPUtimes(0.5,d1,cosA);

			% now the integral (but no dx*dy)
			% B(n) = sum(A .* exp(1i .* k .* dabs) .* cosA ./ dabs);
			d3 = 1i .* k .* dabs;
			GPUexp(d3,d3);
			GPUtimes(A,d3,d3);
			GPUtimes(d3,cosA,d3);
			GPUrdivide(d3,dabs,d3);
			B(n) = sum(d3);

			% display status, yep I'm still alive
			if dispCnt_ == 1
				displayMessage(3, [num2str(n) ' pixels of ' num2str(size(B,2))]);
				dispCnt_ = options.display;
			else
				dispCnt_ = dispCnt_ - 1;
			end

			% save some data, in case I'll die
			if saveCnt_ == 1
				saveCnt_ = options.save;
				saves = saves + 1;
				displayMessage(3, ['Saving data from ' num2str(n-options.save+1) ' to ' num2str(n)]);
				fname = [options.dataFolder filesep 'data' num2str(saves) '.mat'];
				varname = ['B' num2str(saves)];
				eval([varname '=double(B(n-saveCnt_+1:n));']);
				save(fname,varname,'n');
				eval(['clear ' varname]);
			else
				saveCnt_ = saveCnt_ - 1;
			end
		end % end for
	else
		for n = options.startIndex : N				% more general code for both CPU and GPU

			% path lengths
			dx = BXx(n) - AXx;
			dy = BXy(n) - AXy;
			dz = BXz(n) - AXz;
			dabs = sqrt(dx.*dx + dy.*dy + dz.*dz);

			% compute cos(angle(r01,n)) for the aperture: cos() = a.b /(|a|*|b|)
			% and the obliquity factor 1/2*(1+cos(theta)) for Kirchhoff theory
			cosA = (dx.*ANx + dy.*ANy + dz.*ANz) ./ (dabs .* ANabs);
			cosA = 1 ./ 2 .* (1 + cosA);

			% now the integral (but no dx*dy)
			B(n) = sum(A .* exp(1i .* k .* dabs) .* cosA ./ dabs);

			% display status, yep I'm still alive
			if dispCnt_ == 1
				displayMessage(3, [num2str(n) ' pixels of ' num2str(size(B,2))]);
				dispCnt_ = options.display;
			else
				dispCnt_ = dispCnt_ - 1;
			end

			% save some data, in case I'll die
			if saveCnt_ == 1
				saveCnt_ = options.save;
				saves = saves + 1;
				displayMessage(3, ['Saving data from ' num2str(n-options.save+1) ' to ' num2str(n)]);
				fname = [options.dataFolder filesep 'data' num2str(saves) '.mat'];
				varname = ['B' num2str(saves)];
				eval([varname '=double(B(n-saveCnt_+1:n));']);
				save(fname,varname,'n');
				eval(['clear ' varname]);
			else
				saveCnt_ = saveCnt_ - 1;
			end
		end % end for		
	end % end if
end % end if

% done
displayMessage(3, 'Loop ended'); 

% normalize B; B * (areaA) / (areaB)
Bnorm = double(B) ./ (1i * objA.wavelength) .* ...
    (objA.resolution(1) * objA.resolution(2)) ./ ...
	(objB.resolution(1) * objB.resolution(2));

% think of the aperture
B = zeros(1,numel(objB.grid(1,:)));
B(aperture) = Bnorm;

%% save results & delete old data
delete([options.dataFolder filesep '*data*.mat']);     % delete old data

Contact us