Atmospheric Ray Tracing

- simulation engine now hold rays as shared pointers

include/ITAPropagationPathSim/AtmosphericRayTracing/Simulation/Engine.h

@@ -29,6 +29,10 @@
 #include <ITAPropagationPathSim/AtmosphericRayTracing/Rays.h>
 #include <ITAGeo/Atmosphere/StratifiedAtmosphere.h>
 
+// STD
+#include <vector>
+#include <memory>
+
 namespace ITAPropagationPathSim
 {
 	namespace AtmosphericRayTracing
@@ -44,7 +48,7 @@ namespace ITAPropagationPathSim
 				VistaVector3D v3SourcePosition;
 				std::vector<VistaVector3D> v3RayDirections;
 			private:
-				std::vector<CRay> rays;
+				std::vector<std::shared_ptr<CRay>> rays;
 				//const IAbortCriterion& iAbortCriterion; //!< Reference to externally defined abort criterion.
 
 			public:
@@ -52,7 +56,7 @@ namespace ITAPropagationPathSim
 				virtual ~CEngine() {}
 
 			public:
-				std::vector<CRay> Run(const ITAGeo::CStratifiedAtmosphere& atmosphere);
+				std::vector<std::shared_ptr<CRay>> Run(const ITAGeo::CStratifiedAtmosphere& atmosphere);
 			private:
 				void InitRays();
 			};

src/ITAPropagationPathSim/AtmosphericRayTracing/Simulation/Engine.cpp

@@ -18,12 +18,12 @@ namespace ITAPropagationPathSim { namespace AtmosphericRayTracing { namespace Si
 class CWorker
  {
  private:
-	 CRay& rRay;
+	 std::shared_ptr<CRay> pRay;
 	 const IAbortCriterion& rAbortCriterion; //!< Reference to externally defined abort criterion.
 	 CAdaptiveSolver mSolver;
  public:
-	 CWorker(CRay& ray, const SimulationSettings& simSettings, const IAbortCriterion& abortCriterion) :
-		 rRay(ray), rAbortCriterion(abortCriterion), mSolver(CAdaptiveSolver(simSettings)) {}
+	 CWorker(std::shared_ptr<CRay> ray, const SimulationSettings& simSettings, const IAbortCriterion& abortCriterion) :
+		 pRay(ray), rAbortCriterion(abortCriterion), mSolver(CAdaptiveSolver(simSettings)) {}
 
  private:
 	 bool GroundReflectionOccured(const double& rz1, const double& rz2) const
@@ -41,50 +41,50 @@ class CWorker
 	 }
 	 void ExtendRayByOnePeriod()
 	 {
-		 std::vector<unsigned int> iReflectionIndices = rRay.RelectionIndices();
-		 int reflectionOrder = rRay.ReflectionOrder();
+		 const std::vector<unsigned int> iReflectionIndices = pRay->RelectionIndices();
+		 const int reflectionOrder = pRay->ReflectionOrder();
 
 		 if (reflectionOrder < 2)
 			 ITA_EXCEPT_INVALID_PARAMETER("Cannot extend a ray with a reflection order below 2.");
 
-		 int idxStartReflection = iReflectionIndices[reflectionOrder-2];
-		 int idxEndReflection = iReflectionIndices[reflectionOrder - 1];
+		 const int idxStartReflection = iReflectionIndices[reflectionOrder-2];
+		 const int idxEndReflection = iReflectionIndices[reflectionOrder - 1];
 
-		 double tOffset = rRay.TimeStamps()[idxEndReflection] - rRay.TimeStamps()[idxStartReflection];
+		 double tOffset = pRay->TimeStamps()[idxEndReflection] - pRay->TimeStamps()[idxStartReflection];
 
 		 for (int idx = idxStartReflection+1; idx < idxEndReflection; idx++)
 		 {
-			 const double t = rRay.TimeStamps()[idx] + tOffset;
-			 const VistaVector3D r = rRay.SamplingPoints()[idx];
-			 const VistaVector3D n = rRay.WavefrontNormals()[idx];
-			 rRay.Append(r, n, t);
+			 const double t = pRay->TimeStamps()[idx] + tOffset;
+			 const VistaVector3D r = pRay->SamplingPoints()[idx];
+			 const VistaVector3D n = pRay->WavefrontNormals()[idx];
+			 pRay->Append(r, n, t);
 			 
-			 if (rAbortCriterion.AbortRequested(rRay))
+			 if (rAbortCriterion.AbortRequested(*pRay))
 				 return;
 		 }
 
-		 const double t = rRay.TimeStamps()[idxEndReflection] + tOffset;
-		 const VistaVector3D r = rRay.SamplingPoints()[idxEndReflection];
-		 const VistaVector3D n = rRay.WavefrontNormals()[idxEndReflection];
-		 rRay.AppendReflection(r, n, t);
+		 const double t = pRay->TimeStamps()[idxEndReflection] + tOffset;
+		 const VistaVector3D r = pRay->SamplingPoints()[idxEndReflection];
+		 const VistaVector3D n = pRay->WavefrontNormals()[idxEndReflection];
+		 pRay->AppendReflection(r, n, t);
 	 }
 	 void ExtendRayPeriodically()
 	 {
-		 while (!rAbortCriterion.AbortRequested(rRay))
+		 while (!rAbortCriterion.AbortRequested(*pRay))
 			 ExtendRayByOnePeriod();
 	 }
  public:
 	 void TraceRay(const ITAGeo::CStratifiedAtmosphere& atmosphere)
 	 {
-		 VistaVector3D r = rRay.LastPoint();
+		 VistaVector3D r = pRay->LastPoint();
 		 double rz = r[Vista::Z];
-		 VistaVector3D n = rRay.LastWavefrontNormal();
-		 double time = rRay.LastTimeStamp();
+		 VistaVector3D n = pRay->LastWavefrontNormal();
+		 double time = pRay->LastTimeStamp();
 
 		 VistaVector3D s = ODESolver::NormalToSlowness(n, rz, atmosphere);
 
 		 VistaVector3D rNew, sNew;
-		 while ( !rAbortCriterion.AbortRequested(rRay) )
+		 while ( !rAbortCriterion.AbortRequested(*pRay) )
 		 {
 			 mSolver.Process(r, s, atmosphere, rNew, sNew);
 			 n = ODESolver::SlownessToNormal(s);
@@ -95,12 +95,12 @@ class CWorker
 				 VistaVector3D rGround;
 				 double dtGround;
 				 InterpolateToReflectionPoint(r, rNew, dt, rGround, dtGround);
-				 rRay.AppendReflection(rGround, n, time + dtGround);
+				 pRay->AppendReflection(rGround, n, time + dtGround);
 
-				 if ( rAbortCriterion.AbortRequested(rRay) )
+				 if ( rAbortCriterion.AbortRequested(*pRay) )
 					 return;
 
-				 if (rRay.ReflectionOrder() >= 2) //Ray is periodic
+				 if (pRay->ReflectionOrder() >= 2) //Ray is periodic
 				 {
 					 ExtendRayPeriodically();
 					 return;
@@ -110,20 +110,20 @@ class CWorker
 			 r = rNew;
 			 s = sNew;
 			 time += dt;			 
-			 rRay.Append(r, n, time);
+			 pRay->Append(r, n, time);
 		 }
 	 }
  };
 }}}
 
 
-std::vector<CRay> CEngine::Run(const ITAGeo::CStratifiedAtmosphere& atmosphere)
+std::vector<std::shared_ptr<CRay>> CEngine::Run(const ITAGeo::CStratifiedAtmosphere& atmosphere)
 {
 	InitRays();
 	for (int idx = 0; idx < rays.size(); idx++)
 	{
-		CRay& ray = rays[idx];
-		auto worker = CWorker(ray, settings, iAbortCriterion);
+		//CRay& ray = rays[idx];
+		auto worker = CWorker(rays[idx], settings, iAbortCriterion);
 		worker.TraceRay(atmosphere);
 	}
 	return rays;
@@ -132,5 +132,5 @@ void CEngine::InitRays()
 {
 	rays.clear();
 	for each (const VistaVector3D & v3Direction in v3RayDirections)
-		rays.push_back(CRay(v3SourcePosition, v3Direction));
+		rays.push_back( std::make_shared<CRay>(v3SourcePosition, v3Direction) );
 }
\ No newline at end of file