ART - Simulation Engine

- now if a point is directly on the ground (e.g. source point) no additional point is added for the reflection.

include/ITAPropagationPathSim/AtmosphericRayTracing/Rays.h

@@ -101,6 +101,8 @@ namespace ITAPropagationPathSim
 			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);
+			//! Adds the last point to the list of reflection points
+			void AddLastPointToReflectionList();
 			inline void SetSpreadingLoss(const double& spreadingLoss) { m_dSpreadingLoss = spreadingLoss; };
 
 			//! Returns true, if both rays have the same initial direction

src/ITAPropagationPathSim/AtmosphericRayTracing/Rays.cpp

@@ -54,6 +54,10 @@ void CRay::AppendReflection(const VistaVector3D& position, const VistaVector3D&
 	Append(position, wavefrontNormal, timeStamp);
 	m_viReflectionIndices.push_back(NumPoints() - 1);
 }
+void CRay::AddLastPointToReflectionList()
+{
+	m_viReflectionIndices.push_back(NumPoints() - 1);
+}
 
 
 bool CRay::SameDirection(const CRay& other) const

src/ITAPropagationPathSim/AtmosphericRayTracing/Simulation/Engine.cpp

@@ -24,26 +24,28 @@ 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)
+	 inline 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
+	 inline 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)
+	 inline void InterpolateToReflectionPoint(const VistaVector3D& r1, const VistaVector3D& r2, const double& dt, VistaVector3D& rReflection, double& dtReflection)
 	 {
 		 const VistaVector3D dr = r2 - r1;
+		 if (dr[Vista::Z] == 0)
+			 ITA_EXCEPT_INVALID_PARAMETER("Both points are already at same altitude.");
 		 const double portion = -r1[Vista::Z] / dr[Vista::Z];
 
 		 rReflection = r1 + dr * portion;
 		 dtReflection = dt * portion;
 	 };
-	 void ExtendRayByOnePeriod()
+	 inline void ExtendRayByOnePeriod()
 	 {
 		 const std::vector<int>& iReflectionIndices = pRay->ReflectionIndices();
 		 const int reflectionOrder = pRay->ReflectionOrder();
@@ -80,12 +82,12 @@ class CWorker
 		 pRay->AppendReflection(r, n, t);
 		 rExternalWatcher.ProcessRay(pRay);
 	 };
-	 void ExtendRayPeriodically()
+	 inline void ExtendRayPeriodically()
 	 {
 		 while ( !rExternalWatcher.AbortRequested(pRay) )
 			 ExtendRayByOnePeriod();
 	 };
-	 void CalculateRay(const ITAGeo::CStratifiedAtmosphere& atmosphere)
+	 inline void CalculateRay(const ITAGeo::CStratifiedAtmosphere& atmosphere)
 	 {
 		 VistaVector3D r = pRay->LastPoint();
 		 double rz = r[Vista::Z];
@@ -103,12 +105,17 @@ class CWorker
 
 			 if (GroundReflectionOccured(r[Vista::Z], rNew[Vista::Z]))
 			 {
-				 VistaVector3D rGround;
-				 double dtGround;
-				 InterpolateToReflectionPoint(r, rNew, dt, rGround, dtGround);
-				 rGround[Vista::Z] = 0; //Making sure z-component is truely zero
-				 pRay->AppendReflection(rGround, n, time + dtGround);
-				 rExternalWatcher.ProcessRay(pRay);
+				 if ( std::abs(r[Vista::Z]) < DBL_EPSILON ) //Last point was excactly at z = 0: => Interpolation would lead to a second point at ground
+					 pRay->AddLastPointToReflectionList();
+				 else
+				 {
+					 VistaVector3D rGround;
+					 double dtGround;
+					 InterpolateToReflectionPoint(r, rNew, dt, rGround, dtGround);
+					 rGround[Vista::Z] = 0; //Making sure z-component is truely zero
+					 pRay->AppendReflection(rGround, n, time + dtGround);
+					 rExternalWatcher.ProcessRay(pRay);
+				 }
 
 				 if (rExternalWatcher.AbortRequested(pRay))
 					 return;
@@ -128,7 +135,7 @@ class CWorker
 		 }
 	 };
  public:
-	 void TraceRay(const ITAGeo::CStratifiedAtmosphere& atmosphere)
+	 inline void TraceRay(const ITAGeo::CStratifiedAtmosphere& atmosphere)
 	 {
 		 CalculateRay(atmosphere);
 		 rExternalWatcher.FinalizeRay(pRay);