ARTEngine.cpp 5.74 KB
Newer Older
Philipp Schäfer's avatar
Philipp Schäfer committed
1
#include <ITAPropagationPathSim/AtmosphericRayTracing/Simulation/Engine.h>
Philipp Schäfer's avatar
Philipp Schäfer committed
2
3

// ITA includes
Philipp Schäfer's avatar
Philipp Schäfer committed
4
#include "AdaptiveSolver.h"
Philipp Schäfer's avatar
Philipp Schäfer committed
5
#include <ITAPropagationPathSim/AtmosphericRayTracing/ODESolver/ODESolver.h>
Philipp Schäfer's avatar
Philipp Schäfer committed
6
7

// STD
8
#include <cmath>
Philipp Schäfer's avatar
Philipp Schäfer committed
9

Philipp Schäfer's avatar
Philipp Schäfer committed
10
11
12
13
14
// OMP
#ifdef _OPENMP
#include <omp.h>
#endif

Philipp Schäfer's avatar
Philipp Schäfer committed
15
using namespace ITAPropagationPathSim::AtmosphericRayTracing;
16
using namespace ITAPropagationPathSim::AtmosphericRayTracing::Simulation;
Philipp Schäfer's avatar
Philipp Schäfer committed
17

Philipp Schäfer's avatar
Philipp Schäfer committed
18
19
20
namespace ITAPropagationPathSim { namespace AtmosphericRayTracing { namespace Simulation {

class CWorker
21
22
 {
 private:
Philipp Schäfer's avatar
Philipp Schäfer committed
23
	 std::shared_ptr<CRay> pRay;
Philipp Schäfer's avatar
Philipp Schäfer committed
24
	 const IExternalWatcher& rExternalWatcher; //!< Reference to externally defined abort criterion.
Philipp Schäfer's avatar
Philipp Schäfer committed
25
	 CAdaptiveSolver solver;
26
 public:
Philipp Schäfer's avatar
Philipp Schäfer committed
27
	 inline CWorker(std::shared_ptr<CRay> ray, const Settings& simSettings, const IExternalWatcher& externalWatcher)
Philipp Schäfer's avatar
Philipp Schäfer committed
28
29
30
31
		 : pRay(ray)
		 , rExternalWatcher(externalWatcher)
		 , solver( CAdaptiveSolver(simSettings) )
	 {};
32

33
 private:
Philipp Schäfer's avatar
Philipp Schäfer committed
34
	 inline bool GroundReflectionOccured(const double& rz1, const double& rz2) const
35
36
	 {
		 return signbit(rz1) != signbit(rz2);
Philipp Schäfer's avatar
Philipp Schäfer committed
37
	 };
Philipp Schäfer's avatar
Philipp Schäfer committed
38
	 inline void InterpolateToReflectionPoint(const VistaVector3D& r1, const VistaVector3D& r2, const double& dt, VistaVector3D& rReflection, double& dtReflection)
39
40
	 {
		 const VistaVector3D dr = r2 - r1;
Philipp Schäfer's avatar
Philipp Schäfer committed
41
42
		 if (dr[Vista::Z] == 0)
			 ITA_EXCEPT_INVALID_PARAMETER("Both points are already at same altitude.");
43
		 const double portion = -r1[Vista::Z] / dr[Vista::Z];
44

45
46
		 rReflection = r1 + dr * portion;
		 dtReflection = dt * portion;
Philipp Schäfer's avatar
Philipp Schäfer committed
47
	 };
Philipp Schäfer's avatar
Philipp Schäfer committed
48
	 inline void ExtendRayByOnePeriod()
Philipp Schäfer's avatar
Philipp Schäfer committed
49
	 {
Philipp Schäfer's avatar
ART    
Philipp Schäfer committed
50
		 const std::vector<int>& iReflectionIndices = pRay->ReflectionIndices();
Philipp Schäfer's avatar
Philipp Schäfer committed
51
		 const int reflectionOrder = pRay->ReflectionOrder();
Philipp Schäfer's avatar
Philipp Schäfer committed
52
53
54
55

		 if (reflectionOrder < 2)
			 ITA_EXCEPT_INVALID_PARAMETER("Cannot extend a ray with a reflection order below 2.");

Philipp Schäfer's avatar
Philipp Schäfer committed
56
		 const int idxStartReflection = iReflectionIndices[reflectionOrder - 2];
Philipp Schäfer's avatar
Philipp Schäfer committed
57
		 const int idxEndReflection = iReflectionIndices[reflectionOrder - 1];
Philipp Schäfer's avatar
Philipp Schäfer committed
58

Philipp Schäfer's avatar
Philipp Schäfer committed
59
60
		 double tOffset = pRay->at(idxEndReflection).timeStamp - pRay->at(idxStartReflection).timeStamp;
		 VistaVector3D rXYOffset = pRay->at(idxEndReflection).position - pRay->at(idxStartReflection).position;
Philipp Schäfer's avatar
Philipp Schäfer committed
61

Philipp Schäfer's avatar
Philipp Schäfer committed
62
		 for (int idx = idxStartReflection + 1; idx < idxEndReflection; idx++)
Philipp Schäfer's avatar
Philipp Schäfer committed
63
		 {
Philipp Schäfer's avatar
Philipp Schäfer committed
64
			 const double t = pRay->at(idx).timeStamp + tOffset;
Philipp Schäfer's avatar
Philipp Schäfer committed
65
66
67
68
69
			 VistaVector3D r = pRay->at(idx).position + rXYOffset;
			 VistaVector3D n = pRay->at(idx).wavefrontNormal;
			 r[Vista::Z] = -r[Vista::Z];
			 n[Vista::Z] = -n[Vista::Z];

Philipp Schäfer's avatar
Philipp Schäfer committed
70
			 pRay->Append(r, n, t);
Philipp Schäfer's avatar
Philipp Schäfer committed
71
			 rExternalWatcher.ProcessRay(pRay);
Philipp Schäfer's avatar
Philipp Schäfer committed
72

Philipp Schäfer's avatar
Philipp Schäfer committed
73
			 if (rExternalWatcher.AbortRequested(pRay))
Philipp Schäfer's avatar
Philipp Schäfer committed
74
75
76
				 return;
		 }

Philipp Schäfer's avatar
Philipp Schäfer committed
77
		 const double t = pRay->at(idxEndReflection).timeStamp + tOffset;
Philipp Schäfer's avatar
Philipp Schäfer committed
78
79
80
81
		 VistaVector3D r = pRay->at(idxEndReflection).position + rXYOffset;
		 VistaVector3D n = pRay->at(idxEndReflection).wavefrontNormal;
		 r[Vista::Z] = -r[Vista::Z];
		 n[Vista::Z] = -n[Vista::Z];
Philipp Schäfer's avatar
Philipp Schäfer committed
82
		 pRay->AppendReflection(r, n, t);
Philipp Schäfer's avatar
Philipp Schäfer committed
83
		 rExternalWatcher.ProcessRay(pRay);
Philipp Schäfer's avatar
Philipp Schäfer committed
84
	 };
Philipp Schäfer's avatar
Philipp Schäfer committed
85
	 inline void ExtendRayPeriodically()
Philipp Schäfer's avatar
Philipp Schäfer committed
86
	 {
Philipp Schäfer's avatar
Philipp Schäfer committed
87
		 while ( !rExternalWatcher.AbortRequested(pRay) )
Philipp Schäfer's avatar
Philipp Schäfer committed
88
			 ExtendRayByOnePeriod();
Philipp Schäfer's avatar
Philipp Schäfer committed
89
	 };
Philipp Schäfer's avatar
Philipp Schäfer committed
90
	 inline void CalculateRay(const ITAGeo::CStratifiedAtmosphere& atmosphere)
91
	 {
Philipp Schäfer's avatar
Philipp Schäfer committed
92
		 VistaVector3D r = pRay->LastPoint();
93
		 double rz = r[Vista::Z];
Philipp Schäfer's avatar
Philipp Schäfer committed
94
95
		 VistaVector3D n = pRay->LastWavefrontNormal();
		 double time = pRay->LastTimeStamp();
96
97
98

		 VistaVector3D s = ODESolver::NormalToSlowness(n, rz, atmosphere);

Philipp Schäfer's avatar
Philipp Schäfer committed
99
		 VistaVector3D rNew, sNew;
Philipp Schäfer's avatar
Philipp Schäfer committed
100
		 while ( !rExternalWatcher.AbortRequested(pRay) )
101
		 {
Philipp Schäfer's avatar
Philipp Schäfer committed
102
			 solver.Process(r, s, atmosphere, rNew, sNew);
103
			 n = ODESolver::SlownessToNormal(s);
Philipp Schäfer's avatar
Philipp Schäfer committed
104
			 const double dt = solver.CurrentStepSize();
105

Philipp Schäfer's avatar
Philipp Schäfer committed
106
			 if (GroundReflectionOccured(r[Vista::Z], rNew[Vista::Z]))
Philipp Schäfer's avatar
Philipp Schäfer committed
107
			 {
Philipp Schäfer's avatar
Philipp Schäfer committed
108
109
110
111
112
113
114
115
116
117
118
				 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);
				 }
119

Philipp Schäfer's avatar
Philipp Schäfer committed
120
				 if (rExternalWatcher.AbortRequested(pRay))
Philipp Schäfer's avatar
Philipp Schäfer committed
121
122
					 return;

Philipp Schäfer's avatar
Philipp Schäfer committed
123
				 if (pRay->ReflectionOrder() >= 2) //Ray is periodic
Philipp Schäfer's avatar
Philipp Schäfer committed
124
125
126
127
128
				 {
					 ExtendRayPeriodically();
					 return;
				 }
			 }
129

Philipp Schäfer's avatar
Philipp Schäfer committed
130
131
			 r = rNew;
			 s = sNew;
Philipp Schäfer's avatar
Philipp Schäfer committed
132
			 time += dt;
Philipp Schäfer's avatar
Philipp Schäfer committed
133
			 pRay->Append(r, n, time);
Philipp Schäfer's avatar
Philipp Schäfer committed
134
			 rExternalWatcher.ProcessRay(pRay);
135
		 }
Philipp Schäfer's avatar
Philipp Schäfer committed
136
	 };
Philipp Schäfer's avatar
Philipp Schäfer committed
137
 public:
Philipp Schäfer's avatar
Philipp Schäfer committed
138
	 inline void TraceRay(const ITAGeo::CStratifiedAtmosphere& atmosphere)
Philipp Schäfer's avatar
Philipp Schäfer committed
139
140
141
	 {
		 CalculateRay(atmosphere);
		 rExternalWatcher.FinalizeRay(pRay);
Philipp Schäfer's avatar
Philipp Schäfer committed
142
	 };
143
 };
Philipp Schäfer's avatar
Philipp Schäfer committed
144
145
146
}}}


Philipp Schäfer's avatar
ART    
Philipp Schäfer committed
147
std::vector<std::shared_ptr<CRay>> CEngine::Run(const ITAGeo::CStratifiedAtmosphere& atmosphere, const VistaVector3D& m_v3SourcePosition, const std::vector<VistaVector3D>& v3RayDirections) const
148
149
{
	std::vector<std::shared_ptr<CRay>> rays;
Philipp Schäfer's avatar
Philipp Schäfer committed
150
	rays.reserve(v3RayDirections.size());
151
	for each (const VistaVector3D & v3Direction in v3RayDirections)
Philipp Schäfer's avatar
ART    
Philipp Schäfer committed
152
		rays.push_back(std::make_shared<CRay>(m_v3SourcePosition, v3Direction));
153

Philipp Schäfer's avatar
Philipp Schäfer committed
154
	TraceRays(atmosphere, rays);
155
156
	return rays;
}
Philipp Schäfer's avatar
Philipp Schäfer committed
157
void CEngine::Run(const ITAGeo::CStratifiedAtmosphere& atmosphere, const std::set<std::shared_ptr<CRay>>& rays) const
Philipp Schäfer's avatar
Philipp Schäfer committed
158
{
Philipp Schäfer's avatar
Philipp Schäfer committed
159
160
161
162
163
	TraceRays(atmosphere, std::vector<std::shared_ptr<CRay>>(rays.begin(), rays.end()));
}

void CEngine::TraceRays(const ITAGeo::CStratifiedAtmosphere& atmosphere, const std::vector<std::shared_ptr<CRay>>& rays) const
{
Philipp Schäfer's avatar
Philipp Schäfer committed
164
165
166
167
	std::shared_ptr<IExternalWatcher> simulationWatcher = pExternalWatcher;
	if (!simulationWatcher)
		ITA_EXCEPT_INVALID_PARAMETER("Simulation-Engine: Interface to external watcher is not set. Cannot run simulation.");

Philipp Schäfer's avatar
Philipp Schäfer committed
168
169
	#pragma omp parallel for schedule(static)
	for (int idx = 0; idx < rays.size(); idx++)
Philipp Schäfer's avatar
Philipp Schäfer committed
170
	{
Philipp Schäfer's avatar
Philipp Schäfer committed
171
		if (rays[idx] == nullptr)
Philipp Schäfer's avatar
Philipp Schäfer committed
172
173
			continue;

Philipp Schäfer's avatar
Philipp Schäfer committed
174
		auto worker = CWorker(rays[idx], settings, *simulationWatcher);
Philipp Schäfer's avatar
Philipp Schäfer committed
175
		worker.TraceRay(atmosphere);
Philipp Schäfer's avatar
Philipp Schäfer committed
176
	}
Philipp Schäfer's avatar
Philipp Schäfer committed
177
}