Atmospheric Ray Tracing - Simulation Engine

- Renaming of abort function
- Filled class for integration using time step adaptation with live

include/ITAPropagationPathSim/AtmosphericRayTracing/SimulationEngine.h

@@ -38,7 +38,7 @@ namespace ITAPropagationPathSim
 			{
 			public:
 				//! Returns true if the abort criterion for tracing the given ray is reached
-				virtual bool Check(const CRay& ray) const = 0; //TODO: Rename to Reached() ?
+				virtual bool AbortRequested(const CRay& ray) const = 0;
 			};
 			class ITA_PROPAGATION_PATH_SIM_API AbortAtMaxTime : public IAbortCriterion
 			{
@@ -46,13 +46,14 @@ namespace ITAPropagationPathSim
 				double dTMax;
 			public:
 				AbortAtMaxTime(double tMax = 30) : dTMax(tMax) {}
-				bool Check(const CRay& ray) const { return ray.GetMaxTime() >= dTMax; }
+				bool AbortRequested(const CRay& ray) const { return ray.GetMaxTime() >= dTMax; }
 			};
 
 			struct AdaptiveIntegrationSettings {
 				bool bActive = true;	//!< If this is set to false, the adaptation is bypassed and the integration step size is therefore constant
 				double dMaxError = 0.015;
 				double dUncriticalError = 0.005;
+				unsigned int iMaxAdaptationLevel = 32; //! Maximum times, the time step is halfed
 			};
 
 			struct SimulationSettings {
@@ -67,7 +68,7 @@ namespace ITAPropagationPathSim
 			{
 			private:
 				class CWorker;
-				class CIntegrationStepAdaptor;
+				class CAdaptiveIntegrator;
 
 			public:
 				SimulationSettings settings;

src/ITAPropagationPathSim/AtmosphericRayTracing/SimulationEngine.cpp

@@ -9,14 +9,112 @@
 
 using namespace ITAPropagationPathSim::AtmosphericRayTracing::Simulation;
 
-class CEngine::CIntegrationStepAdaptor
+class CEngine::CAdaptiveIntegrator
 {
 private:
-	const AdaptiveIntegrationSettings settings;
+	const SimulationSettings rSettings;
+	double dCurrentDt = rSettings.dIntegrationTimeStep;
+	unsigned int iDtAdaptationLevel = 0;
+
+	const unsigned int maxDtPortion = 1 << rSettings.AdaptiveIntegration.iMaxAdaptationLevel; //! = 2 ^ iMaxAdaptionLevel
+	unsigned int iTimeFramePortion = 0; //! Portion of processed time of an unadapted time step. Integer between 0 and maxDtPortion.
+	bool bIncreaseDtInNextStep = false; //! Indicates whether the error is small enough that the timestep can be increased during next step.
+	
 public:
-	CIntegrationStepAdaptor(const AdaptiveIntegrationSettings& settings) : settings(settings) {}
+	CAdaptiveIntegrator(const SimulationSettings& settings) : rSettings(settings) {}
+private:
+
+	void IntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere) const
+	{
+		switch (rSettings.solverMethod)
+		{
+		case ODESolver::EULER:
+			ODESolver::Euler(r, s, dCurrentDt, atmosphere);
+			break;
+		case ODESolver::RUNGE_KUTTA:
+			ODESolver::RungeKutta(r, s, dCurrentDt, atmosphere);
+			break;
+		default:
+			ITA_EXCEPT_INVALID_PARAMETER("Unknown solver method.");
+		}
+	}
+	void UpdateTimeFramePortion()
+	{
+		iTimeFramePortion += 1 << (rSettings.AdaptiveIntegration.iMaxAdaptationLevel - iDtAdaptationLevel);
+		if (iTimeFramePortion == maxDtPortion)
+			iTimeFramePortion = 0;
+	}
+	bool IncreasedDtFitsInRestOfTimeFrame() const
+	{
+		const unsigned int newDtPortion = 1 << (rSettings.AdaptiveIntegration.iMaxAdaptationLevel - (iDtAdaptationLevel - 1));
+		return ( (maxDtPortion-iTimeFramePortion) % newDtPortion) == 0;
+	}
+	bool TryIncreaseTimeStep()
+	{
+		if (iDtAdaptationLevel <= 0 || !IncreasedDtFitsInRestOfTimeFrame())
+			return false;
+
+		dCurrentDt *= 2;
+		iDtAdaptationLevel--;
+		return true;
+	}
+	bool TryDecreaseTimeStep()
+	{
+		if (iDtAdaptationLevel >= rSettings.AdaptiveIntegration.iMaxAdaptationLevel)
+			return false;
+
+		dCurrentDt /= 2;
+		iDtAdaptationLevel++;
+		return true;
+	}
+	double EstimateIntegrationError(const double& rz1, const double& rz2, const VistaVector3D& n2, const ITAGeo::CStratifiedAtmosphere& atmosphere) const
+	{
+		const VistaVector3D windVectorError = atmosphere.WindVector(rz2) - atmosphere.WindVector(rz1) - atmosphere.WindVectorGradient(rz1) * (rz2 - rz1);
+		const double error = windVectorError[Vista::X] * n2[Vista::X] + windVectorError[Vista::Y] * n2[Vista::Y]; //Same as .Dot since z-component is zero
+		return std::abs(error);
+	}
+	void AdaptiveIntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere)
+	{
+		if (bIncreaseDtInNextStep)
+			TryIncreaseTimeStep();
+
+		VistaVector3D rNew;
+		VistaVector3D sNew;
+		while (1)
+		{
+			rNew = VistaVector3D(r);
+			sNew = VistaVector3D(s);
+			IntegrationStep(rNew, sNew, atmosphere);
+
+			const double rz = abs(r[Vista::Z]);
+			const double rzNew = abs(rNew[Vista::Z]);
+			const VistaVector3D nNew = ODESolver::SlownessToNormal(sNew);
+
+			const double error = EstimateIntegrationError(rz, rzNew, nNew, atmosphere);
+			if (error <= rSettings.AdaptiveIntegration.dMaxError)
+			{
+				bIncreaseDtInNextStep = error < rSettings.AdaptiveIntegration.dUncriticalError;
+				break;
+			}
+
+			if (!TryDecreaseTimeStep())
+				break;
+		}
+		r = rNew;
+		s = sNew;
+
+		UpdateTimeFramePortion();
+	}
+
 public:
-	//TODO: Functionality of adapting dt according to settings, current ray values (and atmosphere?)
+	double GetCurrentTimeStep() const { return dCurrentDt; }
+	void Process(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere)
+	{
+		if (rSettings.AdaptiveIntegration.bActive)
+			AdaptiveIntegrationStep(r, s, atmosphere);
+		else
+			IntegrationStep(r, s, atmosphere);
+	}
 };
 
 
@@ -24,13 +122,26 @@ class CEngine::CWorker
  {
  private:
 	 CRay& rRay;
-	 const SimulationSettings& rSettings;
-	 const IAbortCriterion& iAbortCriterion; //!< Reference to externally defined abort criterion.
-	 CIntegrationStepAdaptor integrationStepAdaptor;
+	 const IAbortCriterion& rAbortCriterion; //!< Reference to externally defined abort criterion.
+	 CAdaptiveIntegrator mIntegrator;
  public:
 	 CWorker(CRay& ray, const SimulationSettings& simSettings, const IAbortCriterion& abortCriterion) :
-		 rRay(ray), rSettings(simSettings), iAbortCriterion(abortCriterion), integrationStepAdaptor(CIntegrationStepAdaptor(simSettings.AdaptiveIntegration)) {}
+		 rRay(ray), rAbortCriterion(abortCriterion), mIntegrator(CAdaptiveIntegrator(simSettings)) {}
+
  private:
+	 bool ReflectionOccured(const double& rz1, const double& rz2) const
+	 {
+		 return signbit(rz1) != signbit(rz2);
+	 }
+	 void GetReflectionPoint(const VistaVector3D& r1, const VistaVector3D& r2, const double& dt, VistaVector3D& rGround, double& dtGround)
+	 {
+		 const VistaVector3D dr = r2 - r1;
+		 const double timePortion = -r1[Vista::Z] / dr[Vista::Z];
+
+		 rGround = r1 + dr * dt;
+		 rGround[Vista::Z] = 0;
+		 dtGround = dt * timePortion;
+	 }
  public:
 	 void TraceRay(const ITAGeo::CStratifiedAtmosphere& atmosphere)
 	 {
@@ -41,25 +152,22 @@ class CEngine::CWorker
 
 		 VistaVector3D s = ODESolver::NormalToSlowness(n, rz, atmosphere);
 
-		 double dt = rSettings.dIntegrationTimeStep;
-
-		 while ( !iAbortCriterion.Check(rRay) )
+		 while ( !rAbortCriterion.AbortRequested(rRay) )
 		 {
-			 switch (rSettings.solverMethod)
-			 {
-			 case ODESolver::EULER:
-				 ODESolver::Euler(r, s, dt, atmosphere);
-				 break;
-			 case ODESolver::RUNGE_KUTTA:
-				 ODESolver::RungeKutta(r, s, dt, atmosphere);
-				 break;
-			 default:
-				 ITA_EXCEPT_INVALID_PARAMETER("Unknown solver method.");
-			 }
-			 time += dt;
-
-			 rz = abs(r[Vista::Z]);
-			 n = ODESolver::SlownessToNormal(s, rz, atmosphere);
+			 mIntegrator.Process(r, s, atmosphere);
+			 n = ODESolver::SlownessToNormal(s);
+			 const double dt = mIntegrator.GetCurrentTimeStep();
+
+			 //if (ReflectionOccured(1, 1))
+			 //{
+				// VistaVector3D rGround;
+				// double dtGround;
+				// GetReflectionPoint(r, r, dt, rGround, dtGround);
+				// rRay.AddReflection(rGround, n, time + dtGround);
+			 //}
+
+
+			 time += dt;			 
 			 rRay.Append(r, n, time);
 		 }
 	 }