...
 
Commits (4)
......@@ -64,10 +64,26 @@ namespace ITAPropagationPathSim
else
propagationPath.push_back(std::make_shared<ITAGeo::CInhomogeneity>(position, wavefrontNormal, ray[idx].timeStamp));
}
propagationPath.push_back( std::make_shared<ITAGeo::CSensorInhomogeneous>(position, wavefrontNormal, ray.LastTimeStamp()) );
auto pReceiver = std::make_shared<ITAGeo::CSensorInhomogeneous>(position, wavefrontNormal, ray.LastTimeStamp());
pReceiver->dSpreadingLoss = ray.SpreadingLoss();
propagationPath.push_back(pReceiver);
return propagationPath;
}
ITAGeo::CPropagationPathList ToPropagationPath(const std::vector<std::shared_ptr<CRay>>& vpRays)
{
ITAGeo::CPropagationPathList propagationPathList;
for each (std::shared_ptr<CRay> pRay in vpRays)
{
if (!pRay)
continue;
propagationPathList.push_back( ToPropagationPath(*pRay) );
}
return propagationPathList;
}
}
}
}
......
......@@ -154,7 +154,6 @@ EigenraySearch::RayPtr EigenraySearch::CAdaptiveWorker::Run(const ITAGeo::CStrat
PostProcessEigenray(atmosphere);
return pMinDistanceRay;
//TODO: Postprocessing: Calculate spreading loss factor
}
......
......@@ -178,11 +178,11 @@ void CRay::InterpolateToRealMinimumPosition(const VistaVector3D& receiverPos)
VistaVector3D rBefore = rTmpMin;
if (iMinReceiverDistance > 0)
rBefore = ClosestPointOnLineSegmentToReceiver(rTmpMin, at(iMinReceiverDistance + 1).position, receiverPos);
rBefore = ClosestPointOnLineSegmentToReceiver(rTmpMin, at(iMinReceiverDistance - 1).position, receiverPos);
dMin = (rAfter - receiverPos).GetLength();
const float dMinBefore = (rAfter - receiverPos).GetLength();
const float dMinBefore = (rBefore - receiverPos).GetLength();
if (dMinBefore < dMin)
{
iMinReceiverDistance--;
......
......@@ -24,14 +24,17 @@ class CWorker
const IExternalWatcher& rExternalWatcher; //!< Reference to externally defined abort criterion.
CAdaptiveSolver solver;
public:
CWorker(std::shared_ptr<CRay> ray, const Settings& simSettings, const IExternalWatcher& externalWatcher) :
pRay(ray), rExternalWatcher(externalWatcher), solver(CAdaptiveSolver(simSettings)) {}
CWorker(std::shared_ptr<CRay> ray, const Settings& simSettings, const IExternalWatcher& externalWatcher)
: pRay(ray)
, rExternalWatcher(externalWatcher)
, solver( CAdaptiveSolver(simSettings) )
{};
private:
bool GroundReflectionOccured(const double& rz1, const double& rz2) const
{
return signbit(rz1) != signbit(rz2);
}
};
void InterpolateToReflectionPoint(const VistaVector3D& r1, const VistaVector3D& r2, const double& dt, VistaVector3D& rReflection, double& dtReflection)
{
const VistaVector3D dr = r2 - r1;
......@@ -39,7 +42,7 @@ class CWorker
rReflection = r1 + dr * portion;
dtReflection = dt * portion;
}
};
void ExtendRayByOnePeriod()
{
const std::vector<int>& iReflectionIndices = pRay->ReflectionIndices();
......@@ -48,7 +51,7 @@ class CWorker
if (reflectionOrder < 2)
ITA_EXCEPT_INVALID_PARAMETER("Cannot extend a ray with a reflection order below 2.");
const int idxStartReflection = iReflectionIndices[reflectionOrder-2];
const int idxStartReflection = iReflectionIndices[reflectionOrder - 2];
const int idxEndReflection = iReflectionIndices[reflectionOrder - 1];
double tOffset = pRay->at(idxEndReflection).timeStamp - pRay->at(idxStartReflection).timeStamp;
......@@ -76,12 +79,12 @@ class CWorker
n[Vista::Z] = -n[Vista::Z];
pRay->AppendReflection(r, n, t);
rExternalWatcher.ProcessRay(pRay);
}
};
void ExtendRayPeriodically()
{
while (!rExternalWatcher.AbortRequested(pRay))
while ( !rExternalWatcher.AbortRequested(pRay) )
ExtendRayByOnePeriod();
}
};
void CalculateRay(const ITAGeo::CStratifiedAtmosphere& atmosphere)
{
VistaVector3D r = pRay->LastPoint();
......@@ -107,7 +110,7 @@ class CWorker
pRay->AppendReflection(rGround, n, time + dtGround);
rExternalWatcher.ProcessRay(pRay);
if ( rExternalWatcher.AbortRequested(pRay) )
if (rExternalWatcher.AbortRequested(pRay))
return;
if (pRay->ReflectionOrder() >= 2) //Ray is periodic
......@@ -119,17 +122,17 @@ class CWorker
r = rNew;
s = sNew;
time += dt;
time += dt;
pRay->Append(r, n, time);
rExternalWatcher.ProcessRay(pRay);
}
}
};
public:
void TraceRay(const ITAGeo::CStratifiedAtmosphere& atmosphere)
{
CalculateRay(atmosphere);
rExternalWatcher.FinalizeRay(pRay);
}
};
};
}}}
......
......@@ -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