Commit 99fa56f8 authored by Philipp Schäfer's avatar Philipp Schäfer
Browse files

ART - Eigenray Engine Test

- now exports cpp rays via propagation path list
- added matlab function to compare results
parent 9ab6e525
......@@ -41,6 +41,8 @@
#include <stdio.h>
#include <math.h>
#include <omp.h>
using namespace std;
using namespace ITAGeo;
using namespace ITAGeo::Utils;
......@@ -68,8 +70,8 @@ void runTest(const CStratifiedAtmosphere& atmosphere, const VistaVector3D& sourc
std::vector<std::shared_ptr<CRay>> eigenrays = engine.Run(atmosphere, sourcePosition, receiverPosition);
cout << "Starting export..." << endl;
CPropagationPath ray = ToPropagationPath(*eigenrays[0]);
JSON::Export(ray, filename + ".json");
CPropagationPathList raysAsPropPathList = ToPropagationPath(eigenrays);
JSON::Export(raysAsPropPathList, filename + ".json");
cout << "Finished" << endl << endl;
}
......@@ -112,10 +114,13 @@ CStratifiedAtmosphere GetInhomogeneousAtmosphere()
int main(int iNumInArgs, char* pcInArgs[])
{
//Disable multi-threading for debugging purposes
omp_set_num_threads(1);
CStratifiedAtmosphere homAtmosphere = GetHomogeneousAtmosphere();
CStratifiedAtmosphere inhomAtmosphere = GetInhomogeneousAtmosphere();
TestReceiverNearGround(homAtmosphere, "Homogeneous");
TestReceiverNearGround(inhomAtmosphere, "Inhomogeneous");
//TestReceiverNearGround(inhomAtmosphere, "Inhomogeneous");
}
function EigenrayEngineTest()
%% Define Atmospheres
inhomAtmos = itaAtmosphere;
%Settings
%(These are the standard settings, itaAtmosphere is already initialized with these values)
inhomAtmos.windProfile = WindProfile.Log; %Enum: Zero, Const, Log
inhomAtmos.tempProfile = TempProfile.ISA; %Enum: Const, ISA
inhomAtmos.z0 = 0.1; %Surface Roughness for Log Wind Profile [m]
inhomAtmos.kappa = 0.4; %Karman constant for Log Wind Profile []
inhomAtmos.vStar = 0.6; %Friction velocity for Log Wind Profile [m/s]
inhomAtmos.vDirConst = [1 0 0]; %Normal in wind direction []
inhomAtmos.R = 1.4; %Air Gas Constant, used to calculate c
inhomAtmos.gamma = 402.8/1.4; %Ratio of Specific Heats, used to calculate c
homAtmos = inhomAtmos;
homAtmos.windProfile = WindProfile.Zero;
homAtmos.tempProfile = TempProfile.Const;
homAtmos.TConst = 293.15;
homAtmos.constStaticPressure = 101325;
%% Define Ray Tracer
art = itaART;
art.source = [0 0 50]; %Can also be a single itaCoordinates
art.receiver = [0 0 1.8]; %Can also be a single itaCoordinates
art.dt = 0.01; %Integration variable [s]
art.tMax = 15; %Maximum time of tracing [s]
art.bAdaptiveDt = true;
art.maxErrorV = 0.015;
art.uncriticalErrorV = 0.005;
%art: maximum level for dt adaptation = 32
art.bCopyPeriodicRefl = true;
art.bUseSlowness = true;
art.bUseWaitbar = false;
%% Eigenray search properties
art.bLimitEigenraySearchDirections = false; %If true, directions in which an eigenray is looked for will be limited
art.maxReceiverRadius = 1; %Maximum receiver radius used to calculate receiverRadius
art.maxSourceReceiverAngle = 1; %Maximum allowed angle between direct connection and eigenray used to calculate receiverRadius
art.maxEigenraySearchIterations = 30; %Abort criterion for eigenray search: maximum number of adaptations.
art.minEigenraySearchAngle = 0.001; %Abort criterion for eigenray search: minimum angle between rays.
art.rayZoomThresh = 0.5; %Threshold used for improved ray zooming. If this is undershot, the algorithm does not make a decision.
art.angleForGeomSpreading = 0.01; %Angle between rays used to calculate the geometrical spreading
art.maxReflOrder = 3;
%% Run Tests
art.atmosphere = homAtmos;
runTestReceiverNearGround(art, 'Homogeneous');
art.atmosphere = inhomAtmos;
%runTestReceiverNearGround(art, 'Inhomogeneous');
function runTestReceiverNearGround(art, atmosphereSuffix)
sourceReceiverDistance = 2000;
for sourceRecAngle = -90:45:90
art.source.x = sourceReceiverDistance*sind(sourceRecAngle) + art.receiver.x;
art.source.y = art.receiver.y;
art.source.z = sourceReceiverDistance*cosd(sourceRecAngle) + art.receiver.z;
cppFilename = ['ReceiverNearGround_' atmosphereSuffix '_SourceAngle' num2str(sourceRecAngle) '.json'];
runMatlabVsCppEigenrayTest(art, cppFilename, fileparts(cppFilename));
end
function runMatlabVsCppEigenrayTest(art, cppFile, plotTitleBase, doPlotRays)
if nargin == 3
doPlotRays = false;
end
%% Load C++ Results
cppRays = itaAtmosRay.load(cppFile);
%% Ray Tracing
matlabRays = art.findEigenrays;
%% Number of rays
if numel(matlabRays) ~= numel(cppRays)
warning('Number of eigenrays of Matlab and C++ simulation does not match')
end
%% Plot rays
if doPlotRays
plotRays(art, matlabRays, cppRays)
title(plotTitleBase)
end
%% Plot distance
plotAbsoluteAndRelativeDistance(matlabRays, cppRays, plotTitleBase)
function plotAbsoluteAndRelativeDistance(matlabRays, cppRays, titleStr)
%% Init
nRays = min(numel(matlabRays), numel(cppRays));
axAbsoluteDist = axes(figure);
axRelativeDist = axes(figure);
hold(axAbsoluteDist, 'on')
hold(axRelativeDist, 'on')
for idx=1:nRays
%% Error Calculation
cppPoints = cppRays(idx).interpolateToRay(matlabRays(idx), 'nearest');
matlabPoints = matlabRays(idx).applyGroundReflection();
deltaR = cppPoints - matlabPoints;
deltaRAbs = deltaR.r;
pathLength = [1 matlabRays(idx).pathLength(2:matlabRays(idx).numPoints)]';
relDeltaR = deltaRAbs ./ pathLength * 100;
plot(axAbsoluteDist, matlabRays(idx).t, deltaRAbs)
plot(axRelativeDist, matlabRays(idx).t, relDeltaR)
end
%%
title(axAbsoluteDist, [titleStr ' - C++ vs Matlab'])
xlabel(axAbsoluteDist, 't [s]')
ylabel(axAbsoluteDist, 'd [m]')
legend(axAbsoluteDist, {'Direct path', '1st-order reflection'}, 'location', 'northwest');
grid(axAbsoluteDist, 'on')
title(axRelativeDist, [titleStr ' - C++ vs Matlab'])
xlabel(axRelativeDist, 't [s]')
ylabel(axRelativeDist, '\eta d [%]')
legend(axRelativeDist, {'Direct path', '1st-order reflection'}, 'location', 'northwest');
grid(axRelativeDist, 'on')
function plotRays(art, matlabRays, cppRays)
nRays = min(numel(matlabRays), numel(cppRays));
f = figure;
baseLineWidth = 1.5;
for idx=1:nRays
matlabRays(idx).pr(f, {'linestyle', 'linewidth'}, {'-', baseLineWidth*1.25});
hold on
cppRays(idx).pr(f, {'linestyle', 'linewidth'}, {':', baseLineWidth*1.5});
end
% matlabRays.pr(f, {'linestyle', 'color', 'linewidth'}, {'-', [0 0.6 0], baseLineWidth*1.25});
% matlabRungeRay.pr(f, {'linestyle', 'color', 'linewidth'}, {'-.',[0 0 0.6], baseLineWidth});
% cppRays.pr(f, {'linestyle', 'color', 'linewidth'}, {':',[1.0 0.6 0], baseLineWidth*1.5});
% cppRungeRay.pr(f, {'linestyle', 'color', 'linewidth'}, {'--', [0.6 0 0], baseLineWidth*1.25});
art.source.plot('markerfacecolor', 'b')
art.receiver.plot('markerfacecolor', 'g')
legend({'Direct path - Matlab', 'Direct path - Cpp', '1st-order - Matlab', '1st-order - Cpp'}, 'location', 'northwest');
grid on
view(0,0)
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment