...
 
Commits (5)
......@@ -45,7 +45,7 @@ namespace ITAPropagationPathSim
VistaVector3D position, wavefrontNormal;
propagationPath.push_back( std::make_shared<ITAGeo::CEmitterInhomogeneous>(ray.SourcePoint(), ray.InitialDirection()) );
for (int idx = 1; idx < ray.size()-1; idx++)
for (int idx = 1; idx < ray.size(); idx++)
{
position = ray[idx].position;
wavefrontNormal = ray[idx].wavefrontNormal;
......@@ -55,7 +55,13 @@ namespace ITAPropagationPathSim
wavefrontNormal[Vista::Z] = -wavefrontNormal[Vista::Z];
}
if (ray.IsReflectionIdx(idx))
if (idx == ray.size() - 1)
{
auto pReceiver = std::make_shared<ITAGeo::CSensorInhomogeneous>(position, wavefrontNormal, ray[idx].timeStamp);
pReceiver->dSpreadingLoss = ray.SpreadingLoss();
propagationPath.push_back(pReceiver);
}
else if (ray.IsReflectionIdx(idx))
{
auto anchor = std::make_shared<ITAGeo::CSpecularReflectionInhomogeneous>(position, wavefrontNormal, ray[idx].timeStamp);
anchor->v3FaceNormal = VistaVector3D(0, 0, 1);
......@@ -64,9 +70,6 @@ namespace ITAPropagationPathSim
else
propagationPath.push_back(std::make_shared<ITAGeo::CInhomogeneity>(position, wavefrontNormal, ray[idx].timeStamp));
}
auto pReceiver = std::make_shared<ITAGeo::CSensorInhomogeneous>(position, wavefrontNormal, ray.LastTimeStamp());
pReceiver->dSpreadingLoss = ray.SpreadingLoss();
propagationPath.push_back(pReceiver);
return propagationPath;
}
......@@ -74,6 +77,7 @@ namespace ITAPropagationPathSim
ITAGeo::CPropagationPathList ToPropagationPath(const std::vector<std::shared_ptr<CRay>>& vpRays)
{
ITAGeo::CPropagationPathList propagationPathList;
propagationPathList.reserve(vpRays.size());
for each (std::shared_ptr<CRay> pRay in vpRays)
{
if (!pRay)
......
......@@ -7,23 +7,30 @@ using namespace ITAPropagationPathSim::AtmosphericRayTracing;
#pragma region EigenraySearchWatcher
EigenraySearch::CEigenraySearchWatcher::CEigenraySearchWatcher(const V3DConstPtrVector& vv3ReceiverPositions, const double& tMax, const int& maxReflOrder)
: m_vpReceiverPositions(vv3ReceiverPositions)
EigenraySearch::CEigenraySearchWatcher::CEigenraySearchWatcher(const ReceiverDataVector& vReceiverData, const double& tMax)
: m_vReceiverData(vReceiverData)
, m_dTMax(tMax)
, m_iMaxReflectionOrder(maxReflOrder)
{
int maxReflOrder = 0;
for each (const CReceiverData& receiverData in m_vReceiverData)
std::max(maxReflOrder, receiverData.iReflectionOrder);
m_iMaxReflectionOrderRTAbort = maxReflOrder + 1; //Abort Ray Tracing only if the reflection order is 2 above the maximum eigenray reflection order
}
bool EigenraySearch::CEigenraySearchWatcher::AbortRequested(const std::shared_ptr<CRay> pRay) const
{
//TODO: For CAdaptiveWorker: Track whether distance increased and use this for abortion
return pRay->LastTimeStamp() >= m_dTMax || pRay->ReflectionOrder() > m_iMaxReflectionOrder;
return pRay->LastTimeStamp() >= m_dTMax || pRay->ReflectionOrder() > m_iMaxReflectionOrderRTAbort;
}
void EigenraySearch::CEigenraySearchWatcher::ProcessRay(std::shared_ptr<CRay> pRay) const
{
for each (const VistaVector3D* pReceiverPos in m_vpReceiverPositions)
pRay->UpdateMinimumReceiverDistance(*pReceiverPos);
for each (const CReceiverData& receiverData in m_vReceiverData)
{
if (pRay->ReflectionOrder() == receiverData.iReflectionOrder || pRay->ReflectionOrder() == receiverData.iReflectionOrder-1)
pRay->UpdateMinimumReceiverDistance(receiverData.v3Position);
}
}
void EigenraySearch::CEigenraySearchWatcher::FinalizeRay(std::shared_ptr<CRay> pRay) const
......@@ -74,8 +81,8 @@ EigenraySearch::RayPtr EigenraySearch::CWorkerBase::FindMinimumDistanceRay(const
EigenraySearch::CInitialWorker::CInitialWorker(const VistaVector3D& sourcePosition, const VistaVector3D& receiverPosition, const Simulation::Settings& simSettings, const RayTracingAbortSettings& abortSettings)
: CWorkerBase(sourcePosition, receiverPosition, simSettings, abortSettings)
{
V3DConstPtrVector receiverPositions = { &ReceiverPosition(), &MirroredReceiverPosition() };
m_pSimulationWatcher = std::make_shared<CEigenraySearchWatcher>(receiverPositions, rayTracingAbortSettings.maxTime, rayTracingAbortSettings.maxReflectionOrder + 1);
ReceiverDataVector receiverData = { CReceiverData(ReceiverPosition(), 0) , CReceiverData(MirroredReceiverPosition(), 1) };
m_pSimulationWatcher = std::make_shared<CEigenraySearchWatcher>(receiverData, rayTracingAbortSettings.maxTime);
simulationEngine.pExternalWatcher = m_pSimulationWatcher;
}
......@@ -94,20 +101,21 @@ CRayGrid EigenraySearch::CInitialWorker::InitRayGrid(const ITAGeo::CStratifiedAt
}
void EigenraySearch::CInitialWorker::FindMinimumDistanceRays(const std::set<RayPtr>& rays)
{
int maxReflOrder = 0;
for each (const RayPtr & ray in rays)
{
const CRayReceiverData* rayReceiverDistanceData = ray->ReceiverDistanceData( ReceiverPosition() );
if (rayReceiverDistanceData)
maxReflOrder = std::max( rayReceiverDistanceData->reflectionOrder, maxReflOrder );
//NOTE: For now this Worker can only find eigenrays up to reflection order 1
//int maxReflOrder = 0;
//for each (const RayPtr & ray in rays)
//{
// const CRayReceiverData* rayReceiverDistanceData = ray->ReceiverDistanceData( ReceiverPosition() );
// if (rayReceiverDistanceData)
// maxReflOrder = std::max( rayReceiverDistanceData->reflectionOrder, maxReflOrder );
rayReceiverDistanceData = ray->ReceiverDistanceData( MirroredReceiverPosition() );
if (rayReceiverDistanceData)
maxReflOrder = std::max( rayReceiverDistanceData->reflectionOrder, maxReflOrder );
}
maxReflOrder = std::min(maxReflOrder, rayTracingAbortSettings.maxReflectionOrder);
// rayReceiverDistanceData = ray->ReceiverDistanceData( MirroredReceiverPosition() );
// if (rayReceiverDistanceData)
// maxReflOrder = std::max( rayReceiverDistanceData->reflectionOrder, maxReflOrder );
//}
//maxReflOrder = std::min(maxReflOrder, rayTracingAbortSettings.maxReflectionOrder);
//maxReflOrder = std::min(1U, rayTracingAbortSettings.maxReflectionOrder);
int maxReflOrder = std::min(1, rayTracingAbortSettings.maxReflectionOrder);
vpMinDistanceRays.resize(maxReflOrder + 1);
for (int reflOrder = 0; reflOrder <= maxReflOrder; reflOrder++)
......@@ -135,8 +143,8 @@ EigenraySearch::CAdaptiveWorker::CAdaptiveWorker(const CRayGrid& rayGrid, const
fReceiverRadius = fmin(rayAdaptationSettings.accuracy.maxReceiverRadius, fReceiverRadius);
//TODO: Track whether distance increased and use this for abortion
V3DConstPtrVector receiverPositions = { &VirtualReceiverPosition() };
m_pSimulationWatcher = std::make_shared<CEigenraySearchWatcher>(receiverPositions, rayTracingAbortSettings.maxTime, iActiveReflexionOrder + 1);
ReceiverDataVector receiverData = { CReceiverData(VirtualReceiverPosition(), iActiveReflexionOrder) };
m_pSimulationWatcher = std::make_shared<CEigenraySearchWatcher>(receiverData, rayTracingAbortSettings.maxTime);
simulationEngine.pExternalWatcher = m_pSimulationWatcher;
}
......
......@@ -38,7 +38,6 @@
// STD
#include <vector>
#include <set>
#include <map>
#include <memory>
......@@ -50,16 +49,28 @@ namespace ITAPropagationPathSim
typedef std::shared_ptr<CRay> RayPtr;
typedef std::vector< RayPtr > RayVector;
typedef std::vector<const VistaVector3D* > V3DConstPtrVector;
class CReceiverData {
public:
const VistaVector3D& v3Position;
const int iReflectionOrder;
inline CReceiverData(const VistaVector3D& pos, const int reflOrder): v3Position(pos), iReflectionOrder(reflOrder) {};
inline CReceiverData& operator=(const CReceiverData& other)
{
if(this != &other)
*this = CReceiverData(other.v3Position, other.iReflectionOrder);
return *this;
};
};
typedef std::vector<CReceiverData> ReceiverDataVector;
class CEigenraySearchWatcher : public Simulation::IExternalWatcher
{
private:
V3DConstPtrVector m_vpReceiverPositions;
const double m_dTMax;
const int m_iMaxReflectionOrder;
ReceiverDataVector m_vReceiverData;
double m_dTMax;
int m_iMaxReflectionOrderRTAbort;
public:
CEigenraySearchWatcher(const V3DConstPtrVector& vv3ReceiverPositions, const double& tMax, const int& maxReflOrder);
CEigenraySearchWatcher(const ReceiverDataVector& vReceiverData, const double& tMax);
public:
virtual bool AbortRequested(const std::shared_ptr<CRay> pRay) const;
//! Interface function to Simulation::Engine: Updates the minimum receiver distance of the ray after each processing step
......
......@@ -59,9 +59,9 @@ void runTest(const CStratifiedAtmosphere& atmosphere, const VistaVector3D& sourc
engine.eigenraySettings.rayAdaptation.accuracy.maxSourceReceiverAngle = 1;
engine.eigenraySettings.rayAdaptation.accuracy.maxAngleForGeomSpreading = 0.01;
engine.eigenraySettings.rayAdaptation.advancedRayZooming.threshold = 0.5;
engine.eigenraySettings.rayAdaptation.advancedRayZooming.threshold = 1.0;
engine.eigenraySettings.rayTracing.maxReflectionOrder = 3;
engine.eigenraySettings.rayTracing.maxReflectionOrder = 1;
engine.eigenraySettings.rayTracing.maxTime = 15;
......@@ -121,6 +121,6 @@ int main(int iNumInArgs, char* pcInArgs[])
CStratifiedAtmosphere inhomAtmosphere = GetInhomogeneousAtmosphere();
TestReceiverNearGround(homAtmosphere, "Homogeneous");
//TestReceiverNearGround(inhomAtmosphere, "Inhomogeneous");
TestReceiverNearGround(inhomAtmosphere, "Inhomogeneous");
}
......@@ -30,7 +30,7 @@ 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.dt = 0.1; %Integration variable [s]
art.tMax = 15; %Maximum time of tracing [s]
art.bAdaptiveDt = true;
......@@ -51,18 +51,21 @@ art.maxSourceReceiverAngle = 1; %Maximum allowed angle between direc
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.rayZoomThresh = 1.0; %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;
art.maxReflOrder = 1;
%% Run Tests
art.atmosphere = homAtmos;
runTestReceiverNearGround(art, 'Homogeneous');
%% Test Matlab vs Cpp
% art.atmosphere = homAtmos;
% runTestReceiverNearGround(art, 'Homogeneous');
art.atmosphere = inhomAtmos;
%runTestReceiverNearGround(art, 'Inhomogeneous');
% art.atmosphere = inhomAtmos;
% runTestReceiverNearGround(art, 'Inhomogeneous');
%% Check C++ data only
checkResultsReceiverNearGround('Homogeneous');
checkResultsReceiverNearGround('Inhomogeneous');
function runTestReceiverNearGround(art, atmosphereSuffix)
sourceReceiverDistance = 2000;
......@@ -71,5 +74,11 @@ for sourceRecAngle = -90:45:90
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));
runMatlabVsCppEigenrayTest(art, cppFilename, fileparts(cppFilename), true);
end
function checkResultsReceiverNearGround(atmosphereSuffix)
for sourceRecAngle = -90:45:90
cppFilename = ['ReceiverNearGround_' atmosphereSuffix '_SourceAngle' num2str(sourceRecAngle) '.json'];
checkCppEigenrays(cppFilename, fileparts(cppFilename));
end
function checkCppEigenrays(cppFile, plotTitle)
%% Load C++ Results
cppRays = itaAtmosRay.load(cppFile);
%% Plot rays
receiverPos = [0 0 1.8];
nRays = numel(cppRays);
f = figure;
ax = axes(f);
hold(ax, 'on');
legendEntries = cell(size(cppRays));
for idx=1:nRays
receiverDistance = norm(cppRays(idx).r.cart(end,:) - receiverPos);
reflOrder = cppRays(idx).numReflections;
legendEntries{idx} = ['order: ' num2str(reflOrder) ', dist: ' num2str(receiverDistance) 'm'];
cppRays(idx).pr(f, {'linestyle', 'linewidth'}, {'--', 2.25});
end
plot3(cppRays(1).r0.x, cppRays(1).r0.y, cppRays(1).r0.z, ' ok', 'markerfacecolor', 'b')
plot3(receiverPos(1), receiverPos(2), receiverPos(3), ' ok', 'markerfacecolor', [0 0.6 0])
legend(legendEntries, 'location', 'northwest');
grid on
view(0,0)
title(plotTitle)
\ No newline at end of file
......@@ -35,7 +35,7 @@ hold(axRelativeDist, 'on')
for idx=1:nRays
%% Error Calculation
cppPoints = cppRays(idx).interpolateToRay(matlabRays(idx), 'nearest');
cppPoints = cppRays(idx).interpolateToRay(matlabRays(idx), 'linear');
matlabPoints = matlabRays(idx).applyGroundReflection();
deltaR = cppPoints - matlabPoints;
......@@ -71,15 +71,15 @@ 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});
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')
plot3(art.source.x, art.source.y, art.source.z, ' ok', 'markerfacecolor', 'b')
plot3(art.receiver.x, art.receiver.y, art.receiver.z, ' ok', 'markerfacecolor', [0 0.6 0])
legend({'Direct path - Matlab', 'Direct path - Cpp', '1st-order - Matlab', '1st-order - Cpp'}, 'location', 'northwest');
grid on
......