Commit 136f974c authored by Dipl.-Ing. Jonas Stienen's avatar Dipl.-Ing. Jonas Stienen
Browse files
parents c1238de1 e39cdc5b
......@@ -29,8 +29,8 @@ namespace ITAPropagationPathSim
namespace EigenraySearch {
struct ITA_PROPAGATION_PATH_SIM_API RayAdaptationSettings {
struct {
unsigned int maxNAdaptations = 15; //< Abort after X adaptations of the ray resolution
//double minAngleResolutionDeg = 0.001; //< []
int maxNAdaptations = 15; //< Abort after X adaptations of the ray resolution
double minAngleResolutionDeg = 0.001; //< []
} abort;
struct {
......@@ -44,7 +44,7 @@ namespace ITAPropagationPathSim
} advancedRayZooming;
};
struct ITA_PROPAGATION_PATH_SIM_API RayTracingAbortSettings {
unsigned int maxReflectionOrder = 3; //< Maximum considered order of reflections
int maxReflectionOrder = 3; //< Maximum considered order of reflections
double maxTime = 30; //< [s]
};
struct ITA_PROPAGATION_PATH_SIM_API Settings {
......
......@@ -53,20 +53,23 @@ namespace ITAPropagationPathSim
public:
float distance;
int idxMinDist;
int reflectionOrder;
VistaVector3D posMinDist;
inline CRayReceiverData() : idxMinDist(-1), distance(-1), posMinDist(VistaVector3D()) {}
inline CRayReceiverData(const int iMin, const float dMin, const VistaVector3D& rMin = VistaVector3D()) : idxMinDist(iMin), distance(dMin), posMinDist(rMin) {}
inline CRayReceiverData() : idxMinDist(-1), distance(-1), reflectionOrder(-1), posMinDist(VistaVector3D()) {}
inline CRayReceiverData(int iMin, float dMin, int reflOrder, const VistaVector3D& rMin = VistaVector3D()) : idxMinDist(iMin), distance(dMin), reflectionOrder(reflOrder), posMinDist(rMin) {}
};
class ITA_PROPAGATION_PATH_SIM_API CRay : public std::vector<CRayElement>
{
private:
std::vector<unsigned int> iReflectionIndices;
std::vector<int> iReflectionIndices;
double dSpreadingLoss = -1; //!< Spreadingloss at last point of ray (receiver), calculated at the end of Eigenray search.
typedef const VistaVector3D* ReceiverPositionPtr;
std::map<ReceiverPositionPtr, CRayReceiverData> mReceiverDistanceMap;
typedef std::map<ReceiverPositionPtr, CRayReceiverData> ReceiverDistanceMap;
ReceiverDistanceMap mReceiverDistanceMap;
public:
CRay(const VistaVector3D& v3SourcePos, const double& thetaDeg, const double& phiDeg);
......@@ -75,9 +78,9 @@ namespace ITAPropagationPathSim
public:
#pragma region Get Functions
const std::vector<unsigned int>& RelectionIndices() const { return iReflectionIndices; }
const std::vector<int>& ReflectionIndices() const { return iReflectionIndices; }
unsigned int NumPoints() const { return size(); }
int NumPoints() const { return size(); }
const VistaVector3D& SourcePoint() const { return front().position; }
const VistaVector3D& InitialDirection() const { return front().wavefrontNormal; }
......@@ -85,13 +88,19 @@ namespace ITAPropagationPathSim
const VistaVector3D& LastWavefrontNormal() const { return back().wavefrontNormal; }
const double& LastTimeStamp() const { return back().timeStamp; }
unsigned int ReflectionOrder() const { return iReflectionIndices.size(); }
int ReflectionOrder() const { return iReflectionIndices.size(); }
//! Returns the reflection order of the ray element with given index
int ReflectionOrder(const int idx) const;
//! Returns the spreading loss at end point of receiver. If this has not been calculated yet, this returns -1.
inline double SpreadingLoss() const { return dSpreadingLoss; }
#pragma endregion
//! Appends a new element with position, wavefront normal and timestamp to the ray
void Append(const VistaVector3D& position, const VistaVector3D& wavefrontNormal, const double& timeStamp);
//! Appends a new element to the ray and adds its index to iReflectionIndices
void AppendReflection(const VistaVector3D& position, const VistaVector3D& wavefrontNormal, const double& timeStamp);
inline void SetSpreadingLoss(const double& spreadingLoss) { dSpreadingLoss = spreadingLoss; }
//! Returns true, if both rays have the same initial direction
bool SameDirection(const CRay& other) const;
......@@ -110,7 +119,8 @@ namespace ITAPropagationPathSim
//MINIMUM RECEIVER DISTANCE
public:
//! Returns a pointer to ray receiver distance data for the given receiver position. Returns nullptr if data does not exist.
const CRayReceiverData* ReceiverDistanceData(const VistaVector3D& receiverPos);
const CRayReceiverData* ReceiverDistanceData(const VistaVector3D& receiverPos) const;
//! Updates the distance to given receiver using the last point of the ray. Returns true if distance became smaller.
/**
* To be called after each integration step.
......
......@@ -34,14 +34,6 @@ const VistaVector3D& EigenraySearch::CWorkerBase::VirtualReceiverPosition(const
return v3MirroredReceiverPosition;
return v3ReceiverPosition;
}
const VistaVector3D& EigenraySearch::CWorkerBase::VectorToVirtualReceiver(const VistaVector3D& point, const int reflectionOrder) const
{
return VirtualReceiverPosition(reflectionOrder) - point;
}
float EigenraySearch::CWorkerBase::DistanceToVirtualReceiver(const VistaVector3D& point, const int reflectionOrder) const
{
return VectorToVirtualReceiver(point, reflectionOrder).GetLength();
}
inline VistaVector3D ClosestPointOnLineSegmentToReceiver(const VistaVector3D& segmentP1, const VistaVector3D& segmentP2, const VistaVector3D& receiverPoint)
{
......@@ -79,7 +71,7 @@ void EigenraySearch::CWorkerBase::InterpolateToRealMinimumPosition(const RayPtr&
}
}
EigenraySearch::RayPtr EigenraySearch::CWorkerBase::FindMinimumDistanceRay(const RayVector& rays, const int reflectionOrder)
EigenraySearch::RayPtr EigenraySearch::CWorkerBase::FindMinimumDistanceRay(const std::set<RayPtr>& rays, const int reflectionOrder)
{
float dMin = _FMAX;
const VistaVector3D& receiverPos = VirtualReceiverPosition(reflectionOrder);
......@@ -102,7 +94,7 @@ EigenraySearch::RayPtr EigenraySearch::CWorkerBase::FindMinimumDistanceRay(const
bool EigenraySearch::CInitialWorker::AbortRequested(const RayPtr& pRay) const
{
return pRay->LastTimeStamp() >= rayTracingAbortSettings.maxTime || pRay->ReflectionOrder() > rayTracingAbortSettings.maxReflectionOrder;
return pRay->LastTimeStamp() >= rayTracingAbortSettings.maxTime || pRay->ReflectionOrder() > rayTracingAbortSettings.maxReflectionOrder+1;
}
void EigenraySearch::CInitialWorker::UpdateMinimumReceiverDistance(RayPtr& pRay) const
{
......@@ -119,6 +111,7 @@ std::vector<CRayGrid> EigenraySearch::CInitialWorker::Run(const ITAGeo::CStratif
{
CRayGrid rayGrid = InitRayGrid(atmosphere);
simulationEngine.Run(atmosphere, rayGrid.UniqueRays());
FindMinimumDistanceRays(rayGrid.UniqueRays());
return FinalizeResult(rayGrid);
}
CRayGrid EigenraySearch::CInitialWorker::InitRayGrid(const ITAGeo::CStratifiedAtmosphere& atmosphere)
......@@ -126,14 +119,22 @@ CRayGrid EigenraySearch::CInitialWorker::InitRayGrid(const ITAGeo::CStratifiedAt
return CEquiangularRayDistribution(v3SourcePosition, 7, 10);
//TODO: Set limits for additional directions according to atmosphere if possible
}
void EigenraySearch::CInitialWorker::FindMinimumDistanceRays(const RayVector& rays)
void EigenraySearch::CInitialWorker::FindMinimumDistanceRays(const std::set<RayPtr>& rays)
{
unsigned int maxReflOrder = 0;
int maxReflOrder = 0;
for each (const RayPtr & ray in rays)
maxReflOrder = std::max( ray->ReflectionOrder(), maxReflOrder );
{
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);
maxReflOrder = std::min(1U, rayTracingAbortSettings.maxReflectionOrder);
//maxReflOrder = std::min(1U, rayTracingAbortSettings.maxReflectionOrder);
vpMinDistanceRays.resize(maxReflOrder + 1);
for (int reflOrder = 0; reflOrder <= maxReflOrder; reflOrder++)
......@@ -154,14 +155,14 @@ std::vector<CRayGrid> EigenraySearch::CInitialWorker::FinalizeResult(const CRayG
bool EigenraySearch::CAdaptiveWorker::AbortRequested(const RayPtr& pRay) const
{
return pRay->LastTimeStamp() >= rayTracingAbortSettings.maxTime || pRay->ReflectionOrder() > rayTracingAbortSettings.maxReflectionOrder;
return pRay->LastTimeStamp() >= rayTracingAbortSettings.maxTime || pRay->ReflectionOrder() > iActiveReflexionOrder+1;
}
EigenraySearch::CAdaptiveWorker::CAdaptiveWorker(const CRayGrid& rayGrid, const VistaVector3D& receiverPosition, const Simulation::Settings& simSettings, const Settings& eigenraySettings, const int activeReflectionOrder)
: CWorkerBase(rayGrid.SourcePosition(), receiverPosition, simSettings, eigenraySettings.rayTracing),
adaptiveRayGrid(rayGrid), iActiveReflexionOrder(activeReflectionOrder), rayAdaptationSettings(eigenraySettings.rayAdaptation)
{
const float sourceReceiverDistance = DistanceToVirtualReceiver(rayGrid.SourcePosition());
const float sourceReceiverDistance = (VirtualReceiverPosition() - rayGrid.SourcePosition()).GetLength();
fReceiverRadius = tan(rayAdaptationSettings.accuracy.maxSourceReceiverAngle) * sourceReceiverDistance;
fReceiverRadius = fmin(rayAdaptationSettings.accuracy.maxReceiverRadius, fReceiverRadius);
}
......@@ -169,14 +170,16 @@ EigenraySearch::CAdaptiveWorker::CAdaptiveWorker(const CRayGrid& rayGrid, const
EigenraySearch::RayPtr EigenraySearch::CAdaptiveWorker::Run(const ITAGeo::CStratifiedAtmosphere& atmosphere)
{
simulationEngine.Run(atmosphere, adaptiveRayGrid.UniqueRays());
pMinDistanceRay = FindMinimumDistanceRay(adaptiveRayGrid.UniqueRays(), iActiveReflexionOrder);
while (!EigenrayAccuracyReached())
{
adaptiveRayGrid.ZoomIntoRay(pMinDistanceRay);
simulationEngine.Run(atmosphere, adaptiveRayGrid.NewRaysOfLastAdaptation());
pMinDistanceRay = FindMinimumDistanceRay(adaptiveRayGrid.UniqueRays(), iActiveReflexionOrder);
nAdaptations++;
}
PostProcessEigenray(atmosphere);
return pMinDistanceRay;
//TODO: Postprocessing: Calculate spreading loss factor
}
......@@ -185,24 +188,17 @@ void EigenraySearch::CAdaptiveWorker::UpdateMinimumReceiverDistance(RayPtr& pRay
{
pRay->UpdateMinimumReceiverDistance( VirtualReceiverPosition() );
//const VistaVector3D deltaVec = VectorToVirtualReceiver(pRay->LastPoint());
//const float receiverDistance = deltaVec.GetLength();
//if (fMinReceiverDistance <= receiverDistance)
// return false;
//fMinReceiverDistance = receiverDistance;
//v3MinReceiverVector = deltaVec;
//iMinReceiverDistanceIdx = pRay->NumPoints() - 1;
//pMinDistanceRay = pRay;
//return true;
//TODO: Track whether distance increased and use this for abortion
}
bool EigenraySearch::CAdaptiveWorker::EigenrayAccuracyReached()
{
if (fMinReceiverDistance <= fReceiverRadius)
const CRayReceiverData* receiverData = pMinDistanceRay->ReceiverDistanceData(VirtualReceiverPosition());
if (receiverData && receiverData->distance <= fReceiverRadius)
return true;
bool abort = nAdaptations > rayAdaptationSettings.abort.maxNAdaptations; //TODO: also check delta angle of RayResolutionAdapter?
bool abort = nAdaptations > rayAdaptationSettings.abort.maxNAdaptations
|| adaptiveRayGrid.MaxAngularResolution() < rayAdaptationSettings.abort.minAngleResolutionDeg;
if (abort)
{
//TODO: Warning: Abort criterion reached!
......@@ -211,14 +207,31 @@ bool EigenraySearch::CAdaptiveWorker::EigenrayAccuracyReached()
}
void EigenraySearch::CAdaptiveWorker::PostProcessEigenray(const ITAGeo::CStratifiedAtmosphere& atmosphere)
{
CalculateEigenraySpreadingLoss(atmosphere);
SetEigenrayEndPoint();
CalculateEigenraySpreadingLoss(atmosphere);
}
void EigenraySearch::CAdaptiveWorker::SetEigenrayEndPoint()
{
const CRayReceiverData* receiverData = pMinDistanceRay->ReceiverDistanceData(VirtualReceiverPosition());
const int idxBeforeMin = receiverData->idxMinDist;
const VistaVector3D& rMin = receiverData->posMinDist;
//Interpolate to new point of minimum
const VistaVector3D& r1 = pMinDistanceRay->at(idxBeforeMin).position;
const VistaVector3D& r2 = pMinDistanceRay->at(idxBeforeMin + 1).position;
const float pathFraction = (rMin - r1).GetLength() / (r2 - r1).GetLength();
CRayElement newRayElement = pMinDistanceRay->at(idxBeforeMin).Interpolate(pMinDistanceRay->at(idxBeforeMin + 1), pathFraction);
//Remove all points before and add new point
pMinDistanceRay->resize(idxBeforeMin + 2);
pMinDistanceRay->back() = newRayElement;
}
void EigenraySearch::CAdaptiveWorker::CalculateEigenraySpreadingLoss(const ITAGeo::CStratifiedAtmosphere& atmosphere)
{
const double time = pMinDistanceRay->at(iMinReceiverDistanceIdx).timeStamp;
const double zReceiver = std::abs(pMinDistanceRay->at(iMinReceiverDistanceIdx).position[Vista::Z]);
const double zRef = std::abs(pMinDistanceRay->SourcePoint()[Vista::Z]);
const double& time = pMinDistanceRay->back().timeStamp;
const double zReceiver = std::abs( pMinDistanceRay->back().position[Vista::Z] );
const double zRef = std::abs( pMinDistanceRay->SourcePoint()[Vista::Z] );
if (adaptiveRayGrid.MaxAngularResolution() > rayAdaptationSettings.accuracy.maxAngleForGeomSpreading || !adaptiveRayGrid.Is2D())
{
......@@ -234,43 +247,6 @@ void EigenraySearch::CAdaptiveWorker::CalculateEigenraySpreadingLoss(const ITAGe
const double cReceiver = atmosphere.SpeedOfSound(zReceiver);
const double cRef = atmosphere.SpeedOfSound(zRef);
const double spreadingLoss = surfaceRef/surfaceReceiver * cReceiver/cRef;
}
void EigenraySearch::CAdaptiveWorker::SetEigenrayEndPoint()
{
const VistaVector3D& receiverPos = VirtualReceiverPosition();
const VistaVector3D& rTmpMin = pMinDistanceRay->at(iMinReceiverDistanceIdx).position;
VistaVector3D rAfter = rTmpMin;
if (iMinReceiverDistanceIdx < pMinDistanceRay->size() - 1)
rAfter = ClosestPointOnLineSegmentToReceiver(rTmpMin, pMinDistanceRay->at(iMinReceiverDistanceIdx + 1).position, receiverPos);
VistaVector3D rBefore = rTmpMin;
if (iMinReceiverDistanceIdx > 0)
rBefore = ClosestPointOnLineSegmentToReceiver(rTmpMin, pMinDistanceRay->at(iMinReceiverDistanceIdx + 1).position, receiverPos);
int idxBeforeMin = iMinReceiverDistanceIdx;
VistaVector3D rMin = rAfter;
if ((rBefore - receiverPos).GetLength() < (rAfter - receiverPos).GetLength())
{
idxBeforeMin--;
rMin = rBefore;
}
//TODO:
//Everything up to here should be done during finalization of EACH ray.
//Afterwards the ray of minimum distance has to be calculated.
//The part below should be executed during post processing.
//Interpolate to new point of minimum
const VistaVector3D& r1 = pMinDistanceRay->at(idxBeforeMin).position;
const VistaVector3D& r2 = pMinDistanceRay->at(idxBeforeMin + 1).position;
const float pathFraction = (rMin - r1).GetLength() / (r2-r1).GetLength();
CRayElement newRayElement = pMinDistanceRay->at(idxBeforeMin).Interpolate(pMinDistanceRay->at(idxBeforeMin + 1), pathFraction);
//Remove all points before and add new point
pMinDistanceRay->resize(idxBeforeMin + 2);
pMinDistanceRay->at(pMinDistanceRay->size() - 1) = newRayElement;
pMinDistanceRay->SetSpreadingLoss( surfaceRef/surfaceReceiver * cReceiver/cRef );
}
#pragma endregion
......@@ -37,6 +37,7 @@
// STD
#include <vector>
#include <set>
#include <map>
#include <memory>
......@@ -74,10 +75,8 @@ namespace ITAPropagationPathSim
const VistaVector3D& ReceiverPosition() const { return v3ReceiverPosition; }
const VistaVector3D& MirroredReceiverPosition() const { return v3MirroredReceiverPosition; }
const VistaVector3D& VirtualReceiverPosition(const int reflectionOrder) const;
const VistaVector3D& VectorToVirtualReceiver(const VistaVector3D& point, const int reflectionOrder) const;
float DistanceToVirtualReceiver(const VistaVector3D& point, const int reflectionOrder) const;
RayPtr FindMinimumDistanceRay(const RayVector& rays, const int reflectionOrder);
RayPtr FindMinimumDistanceRay(const std::set<RayPtr>& rays, const int reflectionOrder);
void InterpolateToRealMinimumPosition(const RayPtr& pRay, const VistaVector3D& receiverPos, int& iMinReceiverDistanceIdx, float& dMin, VistaVector3D& rMin);
};
......@@ -99,7 +98,7 @@ namespace ITAPropagationPathSim
std::vector<CRayGrid> Run(const ITAGeo::CStratifiedAtmosphere& atmosphere);
private:
CRayGrid InitRayGrid(const ITAGeo::CStratifiedAtmosphere& atmosphere);
void FindMinimumDistanceRays(const RayVector& rays);
void FindMinimumDistanceRays(const std::set<RayPtr>& rays);
std::vector<CRayGrid> FinalizeResult(const CRayGrid& initialRayGrid);
};
......@@ -110,14 +109,10 @@ namespace ITAPropagationPathSim
CAdaptiveRayGrid adaptiveRayGrid;
const int iActiveReflexionOrder;
RayAdaptationSettings rayAdaptationSettings;
float fReceiverRadius;
RayPtr pMinDistanceRay;
float fMinReceiverDistance;
VistaVector3D v3MinReceiverVector;
int iMinReceiverDistanceIdx;
int nAdaptations = 0;
float fReceiverRadius;
public:
//! Interface function to Simulation::Engine: Checks whether ray does not have to be traced anymore.
......@@ -130,13 +125,11 @@ namespace ITAPropagationPathSim
virtual void UpdateMinimumReceiverDistance(RayPtr& pRay) const;
inline const VistaVector3D& VirtualReceiverPosition() const { return CWorkerBase::VirtualReceiverPosition(iActiveReflexionOrder); }
inline const VistaVector3D& VectorToVirtualReceiver(const VistaVector3D& point) const { return CWorkerBase::VectorToVirtualReceiver(point, iActiveReflexionOrder); }
inline float DistanceToVirtualReceiver(const VistaVector3D& point) const { return CWorkerBase::DistanceToVirtualReceiver(point, iActiveReflexionOrder); }
bool EigenrayAccuracyReached();
void PostProcessEigenray(const ITAGeo::CStratifiedAtmosphere& atmosphere);
void CalculateEigenraySpreadingLoss(const ITAGeo::CStratifiedAtmosphere& atmosphere);
void SetEigenrayEndPoint();
void CalculateEigenraySpreadingLoss(const ITAGeo::CStratifiedAtmosphere& atmosphere);
};
}
}
......
......@@ -32,6 +32,19 @@ CRay::CRay(const VistaVector3D& v3SourcePos, const VistaVector3D& v3Direction)
Append(v3SourcePos, v3Direction.GetNormalized(), 0.0);
}
//int CRay::ReflectionOrder(const int idx) const
//{
// if (idx < 0 || idx >= NumPoints())
// return -1;
//
// auto reflOrderIt = std::lower_bound(iReflectionIndices.cbegin(), iReflectionIndices.cend(), idx);
// if (reflOrderIt == iReflectionIndices.cend())
// return NumPoints() - 1;
// return *reflOrderIt;
//}
void CRay::Append(const VistaVector3D& position, const VistaVector3D& wavefrontNormal, const double& timeStamp)
{
push_back( CRayElement(position, wavefrontNormal, timeStamp) );
......@@ -105,11 +118,13 @@ CRay::const_iterator CRay::IteratorAfterTime(const double& time) const
// ----RECEIVER DISTANCE----
// -------------------------
#pragma region Receiver distance
const CRayReceiverData* CRay::ReceiverDistanceData(const VistaVector3D& receiverPos)
const CRayReceiverData* CRay::ReceiverDistanceData(const VistaVector3D& receiverPos) const
{
if (mReceiverDistanceMap.find(&receiverPos) == mReceiverDistanceMap.cend()) //Receiver not yet in map
ReceiverDistanceMap::const_iterator receiverDataIt = mReceiverDistanceMap.find(&receiverPos);
if (receiverDataIt == mReceiverDistanceMap.cend()) //Receiver not yet in map
return nullptr;
return &mReceiverDistanceMap[&receiverPos];
return &receiverDataIt->second;
}
......@@ -118,15 +133,17 @@ bool CRay::UpdateMinimumReceiverDistance(const VistaVector3D& receiverPos)
float distance = (LastPoint() - receiverPos).GetLength();
if (mReceiverDistanceMap.find(&receiverPos) == mReceiverDistanceMap.cend()) //Receiver not yet in map
{
mReceiverDistanceMap[&receiverPos] = CRayReceiverData(NumPoints() - 1, distance);
mReceiverDistanceMap[&receiverPos] = CRayReceiverData(NumPoints() - 1, distance, ReflectionOrder());
return true;
}
if (distance >= mReceiverDistanceMap[&receiverPos].distance)
CRayReceiverData& receiverDistanceData = mReceiverDistanceMap[&receiverPos];
if (distance >= receiverDistanceData.distance)
return false;
mReceiverDistanceMap[&receiverPos].distance = distance;
mReceiverDistanceMap[&receiverPos].idxMinDist = NumPoints() - 1;
receiverDistanceData.distance = distance;
receiverDistanceData.idxMinDist = NumPoints() - 1;
receiverDistanceData.reflectionOrder = ReflectionOrder();
return true;
}
......
......@@ -42,7 +42,7 @@ class CWorker
}
void ExtendRayByOnePeriod()
{
const std::vector<unsigned int> iReflectionIndices = pRay->RelectionIndices();
const std::vector<int>& iReflectionIndices = pRay->ReflectionIndices();
const int reflectionOrder = pRay->ReflectionOrder();
if (reflectionOrder < 2)
......
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