Commit baf383e8 authored by Philipp Schäfer's avatar Philipp Schäfer
Browse files

Atmospheric Ray Tracing

- Added ODESolver function that do not overwrite the current value
- Simulation/AdaptiveSolver now does not overwrite current values but returns new ones
- Simulation/Engine now knows current and new value for each integration step
parent 5e6e7402
......@@ -37,7 +37,7 @@ namespace ITAPropagationPathSim
//! Converts a wavefront normal to a slowness vector for given altitude in stratified atmosphere
VistaVector3D ITA_PROPAGATION_PATH_SIM_API NormalToSlowness(const VistaVector3D& n, const double& rz, const ITAGeo::CStratifiedAtmosphere& atmosphere);
//! Converts a slowness vector to a wavefront normal in an inhomogeneous atmosphere
//! Converts a slowness vector to a wavefront normal
/**
* The wavefront normal has the same direction as the slowness vector (n = s / norm(s))
*/
......@@ -50,6 +50,18 @@ namespace ITAPropagationPathSim
* @param[in,out] s Current slowness vector. Will be updated.
* @param[in] dt Time step for integration
* @param[in] atmosphere Stratified atmosphere
*/
void ITA_PROPAGATION_PATH_SIM_API Euler(VistaVector3D& r, VistaVector3D& s, const double& dt, const ITAGeo::CStratifiedAtmosphere& atmosphere);
//! Solves the ordinary differential equations (ODEs) of ray propagation in a stratified medium using the Euler method
/**
* @param[in] r Current wavefront position.
* @param[in] s Current slowness vector.
* @param[in] dt Time step for integration
* @param[in] atmosphere Stratified atmosphere
* @param[out] rNew New wavefront position.
* @param[out] sNew New slowness vector.
*
* Literature for solver:
* W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical recipes in C: The Art of Scientific Computing. Cambridge University Press, 2nd edition, 1992.
......@@ -57,7 +69,7 @@ namespace ITAPropagationPathSim
* The differential equations are taken from the slowness vector approach described in:
* A. D. Pierce. Acoustics: An Introduction to Its Physical Principles and Applications, volume 20. McGraw-Hill New York, 1981.
*/
void ITA_PROPAGATION_PATH_SIM_API Euler(VistaVector3D& r, VistaVector3D& s, const double& dt, const ITAGeo::CStratifiedAtmosphere& atmosphere);
void ITA_PROPAGATION_PATH_SIM_API Euler(const VistaVector3D& r, const VistaVector3D& s, const double& dt, const ITAGeo::CStratifiedAtmosphere& atmosphere, VistaVector3D& rNew, VistaVector3D& sNew);
//! Solves the ordinary differential equations (ODEs) of ray propagation in a stratified medium using the classical Runge-Kutta method (RK4)
......@@ -66,6 +78,18 @@ namespace ITAPropagationPathSim
* @param[in,out] s Current slowness vector. Will be updated.
* @param[in] dt Time step for integration
* @param[in] atmosphere Stratified atmosphere
*/
void ITA_PROPAGATION_PATH_SIM_API RungeKutta(VistaVector3D& r, VistaVector3D& s, const double& dt, const ITAGeo::CStratifiedAtmosphere& atmosphere);
//! Solves the ordinary differential equations (ODEs) of ray propagation in a stratified medium using the classical Runge-Kutta method (RK4)
/**
* @param[in] r Current wavefront position.
* @param[in] s Current slowness vector.
* @param[in] dt Time step for integration
* @param[in] atmosphere Stratified atmosphere
* @param[out] rNew New wavefront position.
* @param[out] sNew New slowness vector.
*
* Literature for solver:
* W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical recipes in C: The Art of Scientific Computing. Cambridge University Press, 2nd edition, 1992.
......@@ -73,7 +97,7 @@ namespace ITAPropagationPathSim
* The differential equations are taken from the slowness vector approach described in:
* A. D. Pierce. Acoustics: An Introduction to Its Physical Principles and Applications, volume 20. McGraw-Hill New York, 1981.
*/
void ITA_PROPAGATION_PATH_SIM_API RungeKutta(VistaVector3D& r, VistaVector3D& s, const double& dt, const ITAGeo::CStratifiedAtmosphere& atmosphere);
void ITA_PROPAGATION_PATH_SIM_API RungeKutta(const VistaVector3D& r, const VistaVector3D& s, const double& dt, const ITAGeo::CStratifiedAtmosphere& atmosphere, VistaVector3D& rNew, VistaVector3D& sNew);
}
}
}
......
......@@ -65,38 +65,54 @@ void ODEDerivatives(const double& rz, const VistaVector3D& s, const ITAGeo::CStr
void ODESolver::Euler(VistaVector3D& r, VistaVector3D& s, const double& dt, const ITAGeo::CStratifiedAtmosphere& atmosphere)
{
VistaVector3D rNew, sNew;
Euler(r, s, dt, atmosphere, rNew, sNew);
r = rNew;
s = sNew;
}
void ODESolver::Euler(const VistaVector3D& r, const VistaVector3D& s, const double& dt, const ITAGeo::CStratifiedAtmosphere& atmosphere, VistaVector3D& rNew, VistaVector3D& sNew)
{
VistaVector3D dr;
double dsz = 0;
ODEDerivatives(r[Vista::Z], s, atmosphere, dr, dsz);
s[Vista::Z] = s[Vista::Z] + dt * dsz;
r = r + dt * dr;
sNew = s;
sNew[Vista::Z] = s[Vista::Z] + dt * dsz;
rNew = r + dt * dr;
}
void ODESolver::RungeKutta(VistaVector3D& r, VistaVector3D& s, const double& dt, const ITAGeo::CStratifiedAtmosphere& atmosphere)
{
VistaVector3D rNew, sNew;
RungeKutta(r, s, dt, atmosphere, rNew, sNew);
r = rNew;
s = sNew;
}
const VistaVector3D r0 = r;
const double sz0 = s[Vista::Z];
void ODESolver::RungeKutta(const VistaVector3D& r, const VistaVector3D& s, const double& dt, const ITAGeo::CStratifiedAtmosphere& atmosphere, VistaVector3D& rNew, VistaVector3D& sNew)
{
const double sz = s[Vista::Z];
sNew = s;
auto dr = new VistaVector3D[4];
auto dsz = new double[4];
ODEDerivatives(r[Vista::Z], s, atmosphere, dr[0], dsz[0]);
r = r0 + dt/2 * dr[0];
s[Vista::Z] = sz0 + dt/2 * dsz[0];
rNew = r + dt / 2 * dr[0];
sNew[Vista::Z] = sz + dt / 2 * dsz[0];
ODEDerivatives(r[Vista::Z], s, atmosphere, dr[1], dsz[1]);
r = r0 + dt/2 * dr[1];
s[Vista::Z] = sz0 + dt/2 * dsz[1];
rNew = r + dt / 2 * dr[1];
sNew[Vista::Z] = sz + dt / 2 * dsz[1];
ODEDerivatives(r[Vista::Z], s, atmosphere, dr[2], dsz[2]);
r = r0 + dt * dr[2];
s[Vista::Z] = sz0 + dt * dsz[2];
rNew = r + dt * dr[2];
sNew[Vista::Z] = sz + dt * dsz[2];
ODEDerivatives(r[Vista::Z], s, atmosphere, dr[3], dsz[3]);
s[Vista::Z] = sz0 + dt * (dsz[0] + 2*dsz[1] + 2*dsz[2] + dsz[3]) / 6;
r = r0 + dt * (dr[0] + 2*dr[1] + 2*dr[2] + dr[3]) / 6;
sNew[Vista::Z] = sz + dt * (dsz[0] + 2 * dsz[1] + 2 * dsz[2] + dsz[3]) / 6;
rNew = r + dt * (dr[0] + 2 * dr[1] + 2 * dr[2] + dr[3]) / 6;
}
......@@ -9,24 +9,24 @@
using namespace ITAPropagationPathSim::AtmosphericRayTracing;
void Simulation::CAdaptiveSolver::Process(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere)
void Simulation::CAdaptiveSolver::Process(const VistaVector3D& r, const VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere, VistaVector3D& rNew, VistaVector3D& sNew)
{
if (rSettings.AdaptiveIntegration.bActive)
AdaptiveIntegrationStep(r, s, atmosphere);
AdaptiveIntegrationStep(r, s, atmosphere, rNew, sNew);
else
IntegrationStep(r, s, atmosphere);
IntegrationStep(r, s, atmosphere, rNew, sNew);
}
void Simulation::CAdaptiveSolver::IntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere) const
void Simulation::CAdaptiveSolver::IntegrationStep(const VistaVector3D& r, const VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere, VistaVector3D& rNew, VistaVector3D& sNew) const
{
switch (rSettings.solverMethod)
{
case SolverMethod::EULER:
ODESolver::Euler(r, s, dCurrentDt, atmosphere);
ODESolver::Euler(r, s, dCurrentDt, atmosphere, rNew, sNew);
break;
case SolverMethod::RUNGE_KUTTA:
ODESolver::RungeKutta(r, s, dCurrentDt, atmosphere);
ODESolver::RungeKutta(r, s, dCurrentDt, atmosphere, rNew, sNew);
break;
default:
ITA_EXCEPT_INVALID_PARAMETER("Unknown solver method.");
......@@ -34,18 +34,14 @@ void Simulation::CAdaptiveSolver::IntegrationStep(VistaVector3D& r, VistaVector3
}
void Simulation::CAdaptiveSolver::AdaptiveIntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere)
void Simulation::CAdaptiveSolver::AdaptiveIntegrationStep(const VistaVector3D& r, const VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere, VistaVector3D& rNew, VistaVector3D& sNew)
{
if (bIncreaseDtInNextStep)
TryIncreaseTimeStep();
VistaVector3D rNew;
VistaVector3D sNew;
while (1)
{
rNew = VistaVector3D(r);
sNew = VistaVector3D(s);
IntegrationStep(rNew, sNew, atmosphere);
IntegrationStep(r, s, atmosphere, rNew, sNew);
const double rz = abs(r[Vista::Z]);
const double rzNew = abs(rNew[Vista::Z]);
......@@ -61,8 +57,6 @@ void Simulation::CAdaptiveSolver::AdaptiveIntegrationStep(VistaVector3D& r, Vist
if (!TryDecreaseTimeStep())
break;
}
r = rNew;
s = sNew;
UpdateTimeFramePortion();
}
......
......@@ -47,17 +47,26 @@ namespace ITAPropagationPathSim
public:
double GetCurrentStepSize() const { return dCurrentDt; }
void Process(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere);
//! Performs a single integration step based on the given simulation settings.
void Process(const VistaVector3D& r, const VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere, VistaVector3D& rNew, VistaVector3D& sNew);
private:
void IntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere) const;
void AdaptiveIntegrationStep(VistaVector3D& r, VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere);
//! Performs a single integration step without adapting the step size
void IntegrationStep(const VistaVector3D& r, const VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere, VistaVector3D& rNew, VistaVector3D& sNew) const;
//! Performs a single integration step and adapting the step size if necessary
void AdaptiveIntegrationStep(const VistaVector3D& r, const VistaVector3D& s, const ITAGeo::CStratifiedAtmosphere& atmosphere, VistaVector3D& rNew, VistaVector3D& sNew);
bool TryIncreaseTimeStep();
bool TryDecreaseTimeStep();
//! Returns an estimation of the integration error which is used for the adaptation of the step size
double EstimateIntegrationError(const double& rz1, const double& rz2, const VistaVector3D& n2, const ITAGeo::CStratifiedAtmosphere& atmosphere) const;
//! Increases (doubles) time step if requirements are matched and returns true on success
bool TryIncreaseTimeStep();
//! Decreases (halves) time step if requirements are matched and returns true on success
bool TryDecreaseTimeStep();
//! Updates iTimeFramePortion after an adaptive integration step
void UpdateTimeFramePortion();
//! Returns true if increasing (doubling) the step size still allows to stay within the original resolution
bool IncreasedDtFitsInRestOfTimeFrame() const;
};
}
......
......@@ -43,13 +43,14 @@ class CEngine::CWorker
VistaVector3D s = ODESolver::NormalToSlowness(n, rz, atmosphere);
VistaVector3D rNew, sNew;
while ( !rAbortCriterion.AbortRequested(rRay) )
{
mSolver.Process(r, s, atmosphere);
mSolver.Process(r, s, atmosphere, rNew, sNew);
n = ODESolver::SlownessToNormal(s);
const double dt = mSolver.GetCurrentStepSize();
//if (ReflectionOccured(1, 1))
//if (ReflectionOccured(r[Vista::Z], rNew[Vista::Z]))
//{
// VistaVector3D rGround;
// double dtGround;
......@@ -58,6 +59,8 @@ class CEngine::CWorker
//}
r = rNew;
s = sNew;
time += dt;
rRay.Append(r, n, time);
}
......
Markdown is supported
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