...
 
Commits (156)
......@@ -28,3 +28,4 @@ svnaccess
*.exp
*.skp
*.skb
*.json
......@@ -6,10 +6,17 @@ project( ITAPropagationPathSim )
list( APPEND CMAKE_MODULE_PATH "$ENV{VISTA_CMAKE_COMMON}" )
include( VistaCommon )
if( NOT DEFINED ITA_PROPAGATION_PATH_SIM_WITH_OPENMP )
set( ITA_PROPAGATION_PATH_SIM_WITH_OPENMP ON CACHE BOOL "Build ITA propagation path sim using OpenMP parallelization." )
endif( )
# dependencies
vista_use_package( ITABase REQUIRED FIND_DEPENDENCIES )
vista_use_package( ITAGeo REQUIRED FIND_DEPENDENCIES )
# NOTE: Do not link against ITAPropagationModels. Otherwise, a circular link might be established.
if( ITA_PROPAGATION_PATH_SIM_WITH_OPENMP )
vista_use_package( OpenMP REQUIRED )
endif( )
# includes
include_directories( "include" )
......
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_PROPAGATIONPATHSIM_ART_EIGENRAYSEARCH_ADAPTIVERAYGRID
#define IW_ITA_PROPAGATIONPATHSIM_ART_EIGENRAYSEARCH_ADAPTIVERAYGRID
#include <ITAPropagationPathSim/AtmosphericRayTracing/RayGrid.h>
//STD
#include <set>
namespace ITAPropagationPathSim
{
namespace AtmosphericRayTracing
{
namespace EigenraySearch
{
class ITA_PROPAGATION_PATH_SIM_API CAdaptiveRayGrid : public CRayGrid
{
private:
double m_dMaxDeltaTheta = -1;
double m_dMaxDeltaPhi = -1;
std::set< std::shared_ptr<CRay> > m_vpNewRaysOfLastAdaptation;
public:
inline CAdaptiveRayGrid() {};
CAdaptiveRayGrid(const CRayGrid& rayGrid);
public:
void Init();
void Reset(const CRayGrid& rayGrid);
//! Maximum angular resolution in degrees
inline const double& MaxAngularResolution() { return m_dMaxDeltaTheta < m_dMaxDeltaPhi ? m_dMaxDeltaPhi : m_dMaxDeltaTheta; };
//! Maximum elevation resolution in degrees
inline const double& MaxDeltaTheta() { return m_dMaxDeltaTheta; };
//! Maximum azimuth resolution in degrees
inline const double& MaxDeltaPhi() { return m_dMaxDeltaPhi; };
//! Returns a unique set of the rays which were adding during the last adaptation
inline const std::set< std::shared_ptr<CRay> >& NewRaysOfLastAdaptation() const { return m_vpNewRaysOfLastAdaptation; };
//! Simple adaptation method: Setting the limits to neighboring rays and shoots additional rays to double the angular resolution of the grid.
void ZoomIntoRay(const std::shared_ptr<CRay> pRay);
//! Advanced adaptation method: Further reduces the angular limits using additional information compared to simple method before doubling the angular resolution.
void ZoomIntoRay(const std::shared_ptr<CRay> pRay, const VistaVector3D& v3ReceiverPosition, const double& dThreshold);
//! Zooms into the given ray using a specific angular resolution by creating up to 8 new rays surrounding the given one.
void ZoomIntoRay(const std::shared_ptr<CRay> pRay, const double& deltaTheta, const double& deltaPhi);
private:
//! Further reduces the limits of the ray grid using additional information on given rays and the receiver position
/**
* If the original ray grid looks like the sketch below where r5 is the ray with minimum distance to the receiver,
* this will decide whether the eigenray is rather above/below and left/right of it. In best case, the number of rays is reduced from 9 to 4.
*
* r1--------r2----------r3
* | c1 / c2 \
* r4------r5-------------r6
* \ c3 \ c4 /
* r7------r8--------r9
*
* Also works for a 1D ray grid (NTheta == 1 or NPhi() == 1)
*/
void SetAdvancedRayGridLimits(const std::shared_ptr<CRay> pRay, const int idxMinDist, const VistaVector3D& receiverPosition, const double& threshold);
//! By comparing the vector from three rays to the receiver, this decides between which two rays the eigenray is located and returns their indices
std::vector<int> FindAdvancedRayGridLimits1D(const std::vector< std::shared_ptr<CRay> >& pRays, const int idxMinDist, const VistaVector3D& receiverPosition, const double& threshold) const;
void DoubleRayResolution();
void DoubleThetaResolution();
void DoublePhiResolution();
};
}
}
}
#endif //IW_ITA_PROPAGATIONPATHSIM_ART_EIGENRAYSEARCH_ADAPTIVERAYGRID
\ No newline at end of file
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_PROPAGATIONPATHSIM_ART_EIGENRAYSEARCH_ENGINE
#define IW_ITA_PROPAGATIONPATHSIM_ART_EIGENRAYSEARCH_ENGINE
#include <ITAPropagationPathSim/Definitions.h>
// Vista includes
#include <VistaBase/VistaVector3D.h>
// ITA includes
#include <ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/EigenraySettings.h>
#include <ITAPropagationPathSim/AtmosphericRayTracing/Simulation/ARTSettings.h>
#include <ITAPropagationPathSim/AtmosphericRayTracing/Rays.h>
#include <ITAGeo/Atmosphere/StratifiedAtmosphere.h>
// STD
#include <vector>
#include <memory>
namespace ITAPropagationPathSim
{
namespace AtmosphericRayTracing
{
namespace EigenraySearch {
class ITA_PROPAGATION_PATH_SIM_API CEngine
{
public:
EigenraySearch::Settings eigenraySettings;
Simulation::Settings simulationSettings;
public:
inline CEngine() {};
std::vector<std::shared_ptr<CRay>> Run(const ITAGeo::CStratifiedAtmosphere& atmosphere, const VistaVector3D& sourcePosition, const VistaVector3D& receiverPosition);
};
}
}
}
#endif //IW_ITA_PROPAGATIONPATHSIM_ART_EIGENRAYSEARCH_ENGINE
\ No newline at end of file
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_PROPAGATIONPATHSIM_ART_EIGENRAYSEARCH_SETTINGS
#define IW_ITA_PROPAGATIONPATHSIM_ART_EIGENRAYSEARCH_SETTINGS
#include <ITAPropagationPathSim/Definitions.h>
namespace ITAPropagationPathSim
{
namespace AtmosphericRayTracing
{
namespace EigenraySearch {
struct ITA_PROPAGATION_PATH_SIM_API RayAdaptationSettings {
struct {
int maxNAdaptations = 30; //!< Abort after N adaptations of the ray resolution
double minAngleResolutionDeg = 0.001; //!< Abort if angle between neighboring rays is below this value []
} abort;
struct {
double maxReceiverRadius = 1; //!< Maximum radius of receiver sphere [m]
double maxSourceReceiverAngle = 1; //!< Maximum angle between vectors from source to receiver point and receiver sphere []
double maxAngleForGeomSpreading = 0.01; //!< Angular resolution of rays which is required to calculate spreading loss []
} accuracy;
struct {
bool bActive = false; //!< Switch to enable/disable advanced ray zooming
double threshold = 0.1; //!< Threshold between 0 and 2 above which advanced ray zooming is performed (0 = always, 2 = never)
} advancedRayZooming;
};
struct ITA_PROPAGATION_PATH_SIM_API RayTracingAbortSettings {
int maxReflectionOrder = 1; //!< Maximum considered order of reflections
double maxTime = 30; //!< Maximum propagation time of rays [s]
bool bAbortOnReceiverDistIncrease = true; //!< If enabled, ray tracing will be aborted as soon as ray receiver distance increases
};
struct ITA_PROPAGATION_PATH_SIM_API Settings {
RayTracingAbortSettings rayTracing;
RayAdaptationSettings rayAdaptation;
};
}
}
}
#endif //IW_ITA_PROPAGATIONPATHSIM_ART_EIGENRAYSEARCH_SETTINGS
\ No newline at end of file
# $Id:$
set( RelativeDir "include/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch" )
set( RelativeSourceGroup "Header Files\\ITAPropagationPathSim\\AtmosphericRayTracing\\EigenraySearch" )
set( DirFiles
EigenraySettings.h
EigenrayEngine.h
AdaptiveRayGrid.h
)
set( DirFiles_SourceGroup "${RelativeSourceGroup}" )
set( LocalSourceGroupFiles )
foreach( File ${DirFiles} )
list( APPEND LocalSourceGroupFiles "${RelativeDir}/${File}" )
list( APPEND ProjectSources "${RelativeDir}/${File}" )
endforeach()
source_group( ${DirFiles_SourceGroup} FILES ${LocalSourceGroupFiles} )
\ No newline at end of file
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_PROPAGATIONPATHSIM_ART_ODESOLVER
#define IW_ITA_PROPAGATIONPATHSIM_ART_ODESOLVER
#include <ITAPropagationPathSim/Definitions.h>
// Vista includes
#include <VistaBase/VistaVector3D.h>
// ITA includes
#include <ITAGeo/Atmosphere/StratifiedAtmosphere.h>
namespace ITAPropagationPathSim
{
namespace AtmosphericRayTracing
{
namespace ODESolver
{
//! Converts a wavefront normal to a slowness vector for given altitude in stratified atmosphere
ITA_PROPAGATION_PATH_SIM_API VistaVector3D NormalToSlowness(const VistaVector3D& n, const double& rz, const ITAGeo::CStratifiedAtmosphere& atmosphere);
//! Converts a slowness vector to a wavefront normal
/**
* The wavefront normal has the same direction as the slowness vector (n = s / norm(s))
*/
ITA_PROPAGATION_PATH_SIM_API VistaVector3D SlownessToNormal(const VistaVector3D& s);
//! Solves the ordinary differential equations (ODEs) of ray propagation in a stratified medium using the Euler method
/**
* @param[in,out] r Current wavefront position. Will be updated.
* @param[in,out] s Current slowness vector. Will be updated.
* @param[in] dt Time step for integration
* @param[in] atmosphere Stratified atmosphere
*/
ITA_PROPAGATION_PATH_SIM_API void 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.
*
* 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.
*/
ITA_PROPAGATION_PATH_SIM_API void 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)
/**
* @param[in,out] r Current wavefront position. Will be updated.
* @param[in,out] s Current slowness vector. Will be updated.
* @param[in] dt Time step for integration
* @param[in] atmosphere Stratified atmosphere
*/
ITA_PROPAGATION_PATH_SIM_API void 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.
*
* 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.
*/
ITA_PROPAGATION_PATH_SIM_API void RungeKutta(const VistaVector3D& r, const VistaVector3D& s, const double& dt, const ITAGeo::CStratifiedAtmosphere& atmosphere, VistaVector3D& rNew, VistaVector3D& sNew);
}
}
}
#endif //IW_ITA_PROPAGATIONPATHSIM_ART_ODESOLVER
\ No newline at end of file
# $Id:$
set( RelativeDir "include/ITAPropagationPathSim/AtmosphericRayTracing/ODESolver" )
set( RelativeSourceGroup "Header Files\\ITAPropagationPathSim\\AtmosphericRayTracing\\ODESolver" )
set( DirFiles
ODESolver.h
)
set( DirFiles_SourceGroup "${RelativeSourceGroup}" )
set( LocalSourceGroupFiles )
foreach( File ${DirFiles} )
list( APPEND LocalSourceGroupFiles "${RelativeDir}/${File}" )
list( APPEND ProjectSources "${RelativeDir}/${File}" )
endforeach()
source_group( ${DirFiles_SourceGroup} FILES ${LocalSourceGroupFiles} )
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_RAYGRID
#define IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_RAYGRID
#include <ITAPropagationPathSim/Definitions.h>
// Vista includes
#include <VistaBase/VistaVector3D.h>
// ITA includes
#include <ITAPropagationPathSim/AtmosphericRayTracing/Rays.h>
// STD
#include <vector>
#include <set>
#include <memory>
namespace ITAPropagationPathSim
{
namespace AtmosphericRayTracing
{
//! This class represents a set of rays which is ordered in a grid (matrix) referring to the angles of initial direction (theta and phi)
/**
* This special arrangement of rays allows to find the direct neigbors of any ray in this grid.
*/
class ITA_PROPAGATION_PATH_SIM_API CRayGrid
{
protected:
typedef std::shared_ptr<CRay> RayPtr;
typedef std::vector< RayPtr > RayVector;
typedef std::vector<RayVector> RayMatrix;
VistaVector3D m_v3SourcePos; //!< Origin of all rays in this ray grid.
std::vector<double> m_vdThetaDeg; //!< Sorted vector of elevation angles in degrees.
std::vector<double> m_vdPhiDeg; //!< Sorted vector of azimuth angles in degrees.
private:
RayMatrix m_vvpRayMatrix; //!< Matrix of pointers to ray. First index refers to elevation (theta), second to azimuth (phi) angle.
RayVector m_vpRays; //!< Vector of all rays contained by vvpRayMatrix. May contain duplicates at poles (theta = 0 or 180°).
std::set<RayPtr> m_vpUniqueRays;//!< Unique set of all rays contained by vvpRayMatrix.
bool m_bCircularPhi = false; //!< Indicates whether the azimuth vector is considered to be circular (phi covers full 360 degrees) or not.
private:
CRayGrid(const RayMatrix& rayMatrix, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg, const bool circularPhi = false);
void UpdateDependentRayContainers();
protected:
inline CRayGrid(bool circularPhi = false) : m_bCircularPhi(circularPhi) {};
//! Creates new rays using the angles of initial directions and inserts them into the ray matrix
void InitRays();
public:
//! Creates a set of rays using a source position and the given elevation and azimuth angles [°] for the inital directions. An optional bool indicates whether the azimuth angle covers the full 360°.
/**
* For each combination of theta and phi, one ray is created and stored.
* Note, that both vectors with angles will be sorted before creating the rays if not specified otherwise.
*/
CRayGrid(const VistaVector3D& sourcePos, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg, const bool circularPhi = false, const bool sortAngleVectors = true);
protected:
//! Updates the ray matrix and depenent vectors/set of rays.
void SetRayMatrix(const RayMatrix& newRayMatrix);
//! Filters for certain rays using indices for theta and phi.
void FilterDirections(const std::vector<int>& thetaIdxVec, const std::vector<int>& phiIdxVec, const bool circularPhi = false);
//! Returns a vector with all rays having the same initial elevation as the given ray. Throws an exception if the ray is not part of this grid.
const RayVector& GetRaysWithSameTheta(const std::shared_ptr<CRay> pRay);
//! Returns a vector with all rays having the same initial azimuth as the given ray. Throws an exception if the ray is not part of this grid.
RayVector GetRaysWithSamePhi(const std::shared_ptr<CRay> pRay);
//---Indicing Functions---
protected:
int GetIndex(const std::shared_ptr<CRay> pRay) const;
int IndexToThetaIndex(const int idx) const;
int IndexToPhiIndex(const int idx) const;
void GetNeighboringAngleIndices(const std::shared_ptr<CRay> pRay, std::vector<int>& thetaIdxVec, std::vector<int>& phiIdxVec) const;
//---Booleans---
protected:
bool IsPoleDirection(const double& thetaDeg) const;
bool HasPoleDirection(const std::shared_ptr<CRay> pRay) const;
public:
//! Returns true if the phi angles are defined circular (covering a full circle)
inline bool IsCircular() const { return m_bCircularPhi; };
bool IsEmpty() const;
bool Is2D() const;
bool Contains(const std::shared_ptr<CRay> pRay) const;
//---GET Functions---
protected:
inline const std::vector<std::shared_ptr<CRay>>& Rays() const { return m_vpRays; };
inline const RayMatrix& Matrix() const { return m_vvpRayMatrix; };
public:
inline const std::vector<double>& ThetaDeg() const { return m_vdThetaDeg; };
inline const std::vector<double>& PhiDeg() const { return m_vdPhiDeg; };
inline int NTheta() const { return m_vdThetaDeg.size(); };
inline int NPhi() const { return m_vdPhiDeg.size(); };
inline int NRays() const { return m_vpRays.size(); };
inline const VistaVector3D& SourcePosition() const { return m_v3SourcePos; };
inline const std::set<std::shared_ptr<CRay>>& UniqueRays() const { return m_vpUniqueRays; };
inline const std::shared_ptr<CRay> At(int idxTheta, int idxPhi) const { return m_vvpRayMatrix[idxTheta][idxPhi]; };
//! Calculates and returns an approximation for the surface area of the wavefront at given time by spanning triangles between rays
double WavefrontSurface(const double& time) const;
//! Calculates the surface area for a spherical wavefront at 1m distance using the limits of theta and phi of this grid
double WavefrontSurfaceReference() const;
//! Returns a ray grid containing the rays surrounding the given ray (including this ray).
/**
* Throws an exception if the given ray is not part of this ray grid.
*/
CRayGrid SurroundingGrid(const std::shared_ptr<CRay> pRay) const;
//! Sets the boundaries of this ray grid to the rays surrounding the given ray (including this ray).
/**
* Throws an exception if the given ray is not part of this ray grid.
*/
void SetToSurroundingGrid(const std::shared_ptr<CRay> pRay);
//! Creates a ray grid with the same initial ray directions but new rays
CRayGrid CopyWithNewRays() const;
};
//! A special CRayGrid with an equiangular distribution
class ITA_PROPAGATION_PATH_SIM_API CEquiangularRayDistribution : public CRayGrid
{
public:
//! Creates an equiangular ray distribution with N rays along elevation and azimuth respectively
CEquiangularRayDistribution(const VistaVector3D& sourcePos, const int nAngles);
//! Creates an equiangular ray distribution with N1 rays along elevation and N2 rays along azimuth
CEquiangularRayDistribution(const VistaVector3D& sourcePos, const int nTheta, const int nPhi);
};
}
}
#endif //IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_RAYGRID
\ No newline at end of file
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_PROPAGATIONPATHSIM_ART_RAYS
#define IW_ITA_PROPAGATIONPATHSIM_ART_RAYS
#include <ITAPropagationPathSim/Definitions.h>
// Vista includes
#include <VistaBase/VistaVector3D.h>
// ITA includes
// STD
#include <vector>
#include <map>
namespace ITAPropagationPathSim
{
namespace AtmosphericRayTracing
{
class ITA_PROPAGATION_PATH_SIM_API CRayElement {
public:
VistaVector3D position;
VistaVector3D wavefrontNormal;
double timeStamp;
inline CRayElement() {};
inline CRayElement(const VistaVector3D& r, const VistaVector3D& n, const double& time) : position(r), wavefrontNormal(n), timeStamp(time) {};
//! Does a linear interpolation from this element to the given one using a factor [0.0 1.0] and returns the new element.
CRayElement Interpolate(const CRayElement& target, const float fFraction) const;
};
//! Compares the timestamps of given CRayElements using the < operator
inline bool operator<(const CRayElement& lhs, const CRayElement& rhs) { return lhs.timeStamp < rhs.timeStamp; };
class ITA_PROPAGATION_PATH_SIM_API CRayReceiverData {
public:
bool bDistanceUpdatedInLastIt; //!< indicates whether receiver distance was updated in last during last integration step
float distance;
int idxMinDist;
int reflectionOrder;
VistaVector3D posMinDist;
inline CRayReceiverData() : idxMinDist(-1), distance(-1), reflectionOrder(-1), posMinDist(VistaVector3D()), bDistanceUpdatedInLastIt(false) {};
inline CRayReceiverData(int iMin, float dMin, int reflOrder, const VistaVector3D& rMin = VistaVector3D()) : idxMinDist(iMin), distance(dMin), reflectionOrder(reflOrder), posMinDist(rMin), bDistanceUpdatedInLastIt(true){};
};
class ITA_PROPAGATION_PATH_SIM_API CRay : public std::vector<CRayElement>
{
private:
std::vector<int> m_viReflectionIndices;
double m_dSpreadingLoss = -1; //!< Spreadingloss at last point of ray (receiver), calculated at the end of Eigenray search.
typedef const VistaVector3D* ReceiverPositionPtr;
typedef std::map<ReceiverPositionPtr, CRayReceiverData> ReceiverDistanceMap;
ReceiverDistanceMap m_mReceiverDistanceMap;
public:
CRay(const VistaVector3D& v3SourcePos, const double& thetaDeg, const double& phiDeg);
CRay(const VistaVector3D& v3SourcePos, const VistaVector3D& v3Direction);
public:
#pragma region Get Functions
inline const std::vector<int>& ReflectionIndices() const { return m_viReflectionIndices; };
inline int NumPoints() const { return size(); };
inline const VistaVector3D& SourcePoint() const { return front().position; };
inline const VistaVector3D& InitialDirection() const { return front().wavefrontNormal; };
inline const VistaVector3D& LastPoint() const { return back().position; };
inline const VistaVector3D& LastWavefrontNormal() const { return back().wavefrontNormal; };
inline const double& LastTimeStamp() const { return back().timeStamp; };
inline int ReflectionOrder() const { return m_viReflectionIndices.size(); };
//! Returns the reflection order of the ray element with given index
//int ReflectionOrder(const int idx) const;
//! Returns the spreading loss at end point of receiver. If this has not been calculated yet, this returns -1.
inline double SpreadingLoss() const { return m_dSpreadingLoss; };
#pragma endregion
//! Appends a new element with position, wavefront normal and timestamp to the ray
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
bool SameDirection(const CRay& other) const;
bool IsReflectionIdx(const int idx) const;
//! Returns the ray element at a given time using a linear interpolation
/**
* If the time exceeds the propagation time of the ray, the last element is returned.
*/
CRayElement AtTime(const double& time) const;
//! Returns the ray element closest to a given time using a nearest neighbor approach.
const CRayElement& Nearest(const double& time) const;
private:
CRay::const_iterator IteratorAfterTime(const double& time) const;
//MINIMUM RECEIVER DISTANCE
public:
//! Returns a pointer to ray receiver distance data for the given receiver position. Returns nullptr if data does not exist.
const CRayReceiverData* ReceiverDistanceData(const VistaVector3D& receiverPos) const;
//! Updates the distance to given receiver using the last point of the ray. Returns true if distance became smaller.
/**
* To be called after each integration step.
*/
void UpdateMinimumReceiverDistance(const VistaVector3D& receiverPos);
//! Returns true if the receiver distance has been updated for at least one receiver during the last integration time step.
bool ReceiverDistanceUpdatedOnLastTimeStep();
//! Increases the accuracy of each receiver distance using a linear interpolation to adjacent ray points.
/**
* To be called after the ray tracing has been finished.
*/
void FinalizeMinimumReceiverDistances();
private:
void InterpolateToRealMinimumPosition(const VistaVector3D& receiverPos);
};
}
}
#endif //IW_ITA_PROPAGATIONPATHSIM_ART_RAYS
\ No newline at end of file
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_ENGINE
#define IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_ENGINE
#include <ITAPropagationPathSim/Definitions.h>
// Vista includes
#include <VistaBase/VistaVector3D.h>
// ITA includes
#include <ITAPropagationPathSim/AtmosphericRayTracing/Simulation/ARTSettings.h>
#include <ITAPropagationPathSim/AtmosphericRayTracing/Rays.h>
#include <ITAGeo/Atmosphere/StratifiedAtmosphere.h>
// STD
#include <vector>
#include <set>
#include <memory>
namespace ITAPropagationPathSim
{
namespace AtmosphericRayTracing
{
namespace Simulation {
class ITA_PROPAGATION_PATH_SIM_API CEngine
{
public:
std::shared_ptr<IExternalWatcher> pExternalWatcher; //!< Reference to externally defined abort criterion.
Simulation::Settings settings;
public:
inline CEngine()
{
pExternalWatcher = std::make_shared<CAbortAtMaxTime>();
};
inline CEngine(std::shared_ptr<IExternalWatcher> pWatcher) : pExternalWatcher(pWatcher) {};
public:
std::vector<std::shared_ptr<CRay>> Run(const ITAGeo::CStratifiedAtmosphere& atmosphere, const VistaVector3D& m_v3SourcePosition, const std::vector<VistaVector3D>& v3RayDirections) const;
void Run(const ITAGeo::CStratifiedAtmosphere& atmosphere, const std::set<std::shared_ptr<CRay>>& rays) const;
private:
void TraceRays(const ITAGeo::CStratifiedAtmosphere& atmosphere, const std::vector<std::shared_ptr<CRay>>& rays) const;
};
}
}
}
#endif //IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_ENGINE
\ No newline at end of file
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_SETTINGS
#define IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_SETTINGS
#include <ITAPropagationPathSim/Definitions.h>
// ITA includes
#include <ITAPropagationPathSim/AtmosphericRayTracing/Rays.h>
// STD
#include <memory>
namespace ITAPropagationPathSim
{
namespace AtmosphericRayTracing
{
namespace Simulation {
class ITA_PROPAGATION_PATH_SIM_API IExternalWatcher
{
public:
//! Returns true if the abort criterion for tracing the given ray is reached
virtual bool AbortRequested(const std::shared_ptr<CRay> pRay) const = 0;
//! Allows additional proecssing of the ray at the end of each time step.
virtual void ProcessRay(std::shared_ptr<CRay>) const = 0;
//! This is called after the tracing of a ray is finished.
virtual void FinalizeRay(std::shared_ptr<CRay>) const = 0;
};
class ITA_PROPAGATION_PATH_SIM_API CAbortAtMaxTime : public IExternalWatcher
{
private:
double dTMax;
public:
inline CAbortAtMaxTime(double tMax = 30) : dTMax(tMax) {};
inline bool AbortRequested(const std::shared_ptr<CRay> pRay) const { return pRay->LastTimeStamp() >= dTMax; };
inline virtual void ProcessRay(std::shared_ptr<CRay>) const {};
inline virtual void FinalizeRay(std::shared_ptr<CRay>) const {};
};
enum ITA_PROPAGATION_PATH_SIM_API SolverMethod
{
EULER = 0, //!< Euler method
RUNGE_KUTTA //!< Classical Runge_Kutta method (RK4)
};
struct ITA_PROPAGATION_PATH_SIM_API 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; //!< For errors above this threshold the time step is decreased.
double dUncriticalError = 0.005; //!< An error below this limit allows to increase the time step again.
unsigned int iMaxAdaptationLevel = 31; //!< Maximum times, the time step is halfed. Maximum valid value = 31.
};
struct ITA_PROPAGATION_PATH_SIM_API Settings {
SolverMethod solverMethod = SolverMethod::EULER;
double dIntegrationTimeStep = 0.1;
AdaptiveIntegrationSettings adaptiveIntegration;
};
}
}
}
#endif //IW_ITA_PROPAGATIONPATHSIM_ART_SIMULATION_SETTINGS
\ No newline at end of file
# $Id:$
set( RelativeDir "include/ITAPropagationPathSim/AtmosphericRayTracing/Simulation" )
set( RelativeSourceGroup "Header Files\\ITAPropagationPathSim\\AtmosphericRayTracing\\Simulation" )
set( DirFiles
ARTSettings.h
ARTEngine.h
)
set( DirFiles_SourceGroup "${RelativeSourceGroup}" )
set( LocalSourceGroupFiles )
foreach( File ${DirFiles} )
list( APPEND LocalSourceGroupFiles "${RelativeDir}/${File}" )
list( APPEND ProjectSources "${RelativeDir}/${File}" )
endforeach()
source_group( ${DirFiles_SourceGroup} FILES ${LocalSourceGroupFiles} )
\ No newline at end of file
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_PROPAGATIONPATHSIM_ART_UTILS_RAYTOPROPAGATIONPATH
#define IW_ITA_PROPAGATIONPATHSIM_ART_UTILS_RAYTOPROPAGATIONPATH
#include <ITAPropagationPathSim/Definitions.h>
// ITA includes
#include <ITAPropagationPathSim/AtmosphericRayTracing/Rays.h>
#include <ITAGeo/Base.h>
// STD
#include <memory>
#include <vector>
namespace ITAPropagationPathSim
{
namespace AtmosphericRayTracing
{
namespace Utils
{
ITA_PROPAGATION_PATH_SIM_API ITAGeo::CPropagationPath ToPropagationPath(const CRay& ray);
ITA_PROPAGATION_PATH_SIM_API ITAGeo::CPropagationPathList ToPropagationPath(const std::vector<std::shared_ptr<CRay>>& vpRays);
}
}
}
#endif //IW_ITA_PROPAGATIONPATHSIM_ART_UTILS_RAYTOPROPAGATIONPATH
\ No newline at end of file
# $Id:$
set( RelativeDir "include/ITAPropagationPathSim/AtmosphericRayTracing/Utils" )
set( RelativeSourceGroup "Header Files\\ITAPropagationPathSim\\AtmosphericRayTracing\\Utils" )
set( DirFiles
RayToPropagationPath.h
)
set( DirFiles_SourceGroup "${RelativeSourceGroup}" )
set( LocalSourceGroupFiles )
foreach( File ${DirFiles} )
list( APPEND LocalSourceGroupFiles "${RelativeDir}/${File}" )
list( APPEND ProjectSources "${RelativeDir}/${File}" )
endforeach()
source_group( ${DirFiles_SourceGroup} FILES ${LocalSourceGroupFiles} )
# $Id:$
set( RelativeDir "include/ITAPropagationPathSim/AtmosphericRayTracing" )
set( RelativeSourceGroup "Header Files\\ITAPropagationPathSim\\AtmosphericRayTracing" )
set( SubDirs ODESolver Simulation EigenraySearch Utils)
set( DirFiles
Rays.h
RayGrid.h
)
set( DirFiles_SourceGroup "${RelativeSourceGroup}" )
set( LocalSourceGroupFiles )
foreach( File ${DirFiles} )
list( APPEND LocalSourceGroupFiles "${RelativeDir}/${File}" )
list( APPEND ProjectSources "${RelativeDir}/${File}" )
endforeach()
source_group( ${DirFiles_SourceGroup} FILES ${LocalSourceGroupFiles} )
set( SubDirFiles "" )
foreach( Dir ${SubDirs} )
list( APPEND SubDirFiles "${RelativeDir}/${Dir}/_SourceFiles.cmake" )
endforeach()
foreach( SubDirFile ${SubDirFiles} )
include( ${SubDirFile} )
endforeach()
\ No newline at end of file
#ifndef INCLUDE_WATCHER_DIFFRACTION_LOCATOR
#define INCLUDE_WATCHER_DIFFRACTION_LOCATOR
#include <ITABase/ITAProgress.h>
// ITA includes
#include <ITAGeo/Base.h>
#include <ITAGeo/Halfedge/MeshModel.h>
#include <ITAGeo/Utils.h>
#include<ITAPropagationPathSim/CombinedModel/PropagationShapes.h>
#include <ITAPropagationPathSim/CombinedModel/PropagationShapes.h>
namespace ITAPropagationPathSim
{
......@@ -16,14 +16,14 @@ namespace ITAPropagationPathSim
using namespace ITAGeo;
using namespace std;
namespace Diffraction
namespace Diffraction
{
bool ITA_PROPAGATION_PATH_SIM_API ConstructAperturePoints(shared_ptr<const CEmitter> pEmitter, shared_ptr<const CSensor> pSensor, const int iNumberIterations, const vector<CPropagationShapeShared> pPropagationListsIn, vector<CPropagationShapeShared>& pPropagationListsOut);
bool ITA_PROPAGATION_PATH_SIM_API ConstructAperturePoints( shared_ptr< const CEmitter > pEmitter, shared_ptr< const CSensor > pSensor, const int iNumberIterations, const vector< CPropagationShapeShared > pPropagationListsIn, vector< CPropagationShapeShared >& pPropagationListsOut, ITABase::IProgressHandler* pProgressHandler = nullptr );
//!< Exclude propagation paths with an accumulated diffraction angle that is bigger than the angle threshold.
bool ITA_PROPAGATION_PATH_SIM_API AccumulatedAngleCulling(const float fAngleThreshold, shared_ptr<const CEmitter> pEmitter, const vector<CPropagationShapeShared> pPropagationTreeIn, vector<CPropagationShapeShared>& pPropagationTreeOut);
bool ITA_PROPAGATION_PATH_SIM_API AccumulatedAngleCulling(const float fAngleThreshold, shared_ptr<const CEmitter> pEmitter, shared_ptr<const CSensor> pSensor, const vector<CPropagationShapeShared> pPropagationListsIn, vector<CPropagationShapeShared>& pPropagationListsOut);
bool ITA_PROPAGATION_PATH_SIM_API AccumulatedAngleCulling( const float fAngleThreshold, shared_ptr< const CEmitter > pEmitter, const vector< CPropagationShapeShared > pPropagationTreeIn, vector< CPropagationShapeShared >& pPropagationTreeOut );
bool ITA_PROPAGATION_PATH_SIM_API AccumulatedAngleCulling( const float fAngleThreshold, shared_ptr< const CEmitter > pEmitter, shared_ptr<const CSensor> pSensor, const vector<CPropagationShapeShared> pPropagationListsIn, vector< CPropagationShapeShared >& pPropagationListsOut );
}
}
......
......@@ -2,7 +2,7 @@
set( RelativeDir "include/ITAPropagationPathSim" )
set( RelativeSourceGroup "Header Files\\ITAPropagationPathSim" )
set( SubDirs MirrorImage ConvexDiffraction CombinedModel BaseEngine UrbanEngine)
set( SubDirs MirrorImage ConvexDiffraction CombinedModel BaseEngine UrbanEngine AtmosphericRayTracing)
set( DirFiles
Definitions.h
......
#include <ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/AdaptiveRayGrid.h>
#include <ITAException.h>
#include <math.h>
using namespace ITAPropagationPathSim::AtmosphericRayTracing;
using namespace ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch;
#pragma region CONSTRUCTOR / INITIATION
EigenraySearch::CAdaptiveRayGrid::CAdaptiveRayGrid(const CRayGrid& rayGrid) : CRayGrid(rayGrid)
{
if (rayGrid.UniqueRays().size() <= 1)
ITA_EXCEPT_INVALID_PARAMETER("Ray Grid must atleast contain 2 unique rays");
Init();
}
void EigenraySearch::CAdaptiveRayGrid::Init()
{
m_vpNewRaysOfLastAdaptation.clear();
for (int idx = 1; idx < m_vdThetaDeg.size(); idx++)
{
double deltaTheta = m_vdThetaDeg[idx] - m_vdThetaDeg[idx - 1];
if (m_dMaxDeltaTheta < deltaTheta)
m_dMaxDeltaTheta = deltaTheta;
}
for (int idx = 1; idx < m_vdPhiDeg.size(); idx++)
{
double deltaPhi = m_vdPhiDeg[idx] - m_vdPhiDeg[idx - 1];
if (deltaPhi < 0)
deltaPhi += 360;
if (m_dMaxDeltaPhi < deltaPhi)
m_dMaxDeltaPhi = deltaPhi;
}
}
void CAdaptiveRayGrid::Reset(const CRayGrid& rayGrid)
{
if (rayGrid.UniqueRays().size() <= 1)
ITA_EXCEPT_INVALID_PARAMETER("Ray Grid must atleast contain 2 unique rays");
*this = rayGrid;
Init();
}
#pragma endregion
#pragma region RAY ZOOMING - public
void CAdaptiveRayGrid::ZoomIntoRay(const std::shared_ptr<CRay> pRay)
{
if(!Contains(pRay))
ITA_EXCEPT_INVALID_PARAMETER("Given ray is not part of the adaptive ray grid.");
SetToSurroundingGrid(pRay);
DoubleRayResolution();
}
void CAdaptiveRayGrid::ZoomIntoRay(const std::shared_ptr<CRay> pRay, const VistaVector3D& v3ReceiverPosition, const double& dThreshold)
{
if (!Contains(pRay))
ITA_EXCEPT_INVALID_PARAMETER("Given ray is not part of the adapted grid.");
const CRayReceiverData* receiverData = pRay->ReceiverDistanceData(v3ReceiverPosition);
if(!receiverData)
ITA_EXCEPT_INVALID_PARAMETER("Given ray has no data about given receiver.");
SetToSurroundingGrid(pRay);
SetAdvancedRayGridLimits(pRay, receiverData->idxMinDist, v3ReceiverPosition, dThreshold);
DoubleRayResolution();
}
void CAdaptiveRayGrid::ZoomIntoRay(const std::shared_ptr<CRay> pZoomRay, const double& deltaTheta, const double& deltaPhi)
{
const int idx = GetIndex(pZoomRay);
if (idx < 0)
ITA_EXCEPT_INVALID_PARAMETER("Given ray is not part of the adapted grid.");
if (deltaTheta <= 0 || deltaPhi <= 0)
ITA_EXCEPT_INVALID_PARAMETER("Delta angles must be positive numbers > 0.");
const double thetaDeg = m_vdThetaDeg[ IndexToThetaIndex(idx) ];
const double phiDeg = m_vdPhiDeg[ IndexToPhiIndex(idx) ];
int idxRayTheta = 1; int idxRayPhi = 1;
m_vdThetaDeg = { thetaDeg - deltaTheta, thetaDeg, thetaDeg + deltaTheta };
m_vdPhiDeg = { phiDeg - deltaPhi, phiDeg, phiDeg + deltaPhi };
if (m_vdThetaDeg.front() < 0)
{
m_vdThetaDeg.erase(m_vdThetaDeg.begin());
idxRayTheta = 0;
}
if (m_vdThetaDeg.back() > 180)
m_vdThetaDeg.pop_back();
RayMatrix newRayMatrix(NTheta(), RayVector(NPhi()));
m_vpNewRaysOfLastAdaptation.clear();
for (int idxTheta = 0; idxTheta < m_vdThetaDeg.size(); idxTheta++)
{
for (int idxPhi = 0; idxPhi < m_vdPhiDeg.size(); idxPhi++)
{
RayPtr pCurrentRay;
if (idxRayTheta == idxTheta && idxRayPhi == idxPhi)
pCurrentRay = pZoomRay;
else if ( IsPoleDirection(m_vdThetaDeg[idxTheta]) && (idxRayTheta==idxTheta || idxPhi > 0) )
{
if (idxTheta == idxRayTheta)
pCurrentRay = pZoomRay;
else
pCurrentRay = newRayMatrix[idxTheta].front();
}
else
{
pCurrentRay = std::make_shared<CRay>(m_v3SourcePos, m_vdThetaDeg[idxTheta], m_vdPhiDeg[idxPhi]);
m_vpNewRaysOfLastAdaptation.insert(pCurrentRay);
}
newRayMatrix[idxTheta][idxPhi] = pCurrentRay;
}
}
SetRayMatrix(newRayMatrix);
}
#pragma endregion
#pragma region NEW RAY LIMITS - private
void CAdaptiveRayGrid::SetAdvancedRayGridLimits(const std::shared_ptr<CRay> pRay, const int idxMinDist, const VistaVector3D& v3ReceiverPosition, const double& dThreshold)
{
if (NPhi() > 3) // There can be more than three rays if pRay is located at a pole
return;
std::vector<int> phiIdxVec = FindAdvancedRayGridLimits1D( GetRaysWithSameTheta(pRay), idxMinDist, v3ReceiverPosition, dThreshold );
std::vector<int> thetaIdxVec = FindAdvancedRayGridLimits1D( GetRaysWithSamePhi(pRay), idxMinDist, v3ReceiverPosition, dThreshold );
FilterDirections(thetaIdxVec, phiIdxVec);
}
std::vector<int> CAdaptiveRayGrid::FindAdvancedRayGridLimits1D(const std::vector<std::shared_ptr<CRay>>& vpRayVector, const int idxMinDist, const VistaVector3D& v3ReceiverPosition, const double& dThreshold) const
{
if (vpRayVector.size() == 1)
return { 0 };
if (vpRayVector.size() == 2)
return {0, 1};
if (vpRayVector.size() != 3)
ITA_EXCEPT_INVALID_PARAMETER("Expected a vector with one to three rays.");
const std::shared_ptr<CRay>& pRayMin = vpRayVector[1];
const std::shared_ptr<CRay>& pRay1 = vpRayVector[0];
const std::shared_ptr<CRay>& pRay2 = vpRayVector[2];
const double tMin = pRayMin->at(idxMinDist).timeStamp;
VistaVector3D v1 = (v3ReceiverPosition - pRay1->AtTime(tMin).position).GetNormalized();
VistaVector3D v2 = (v3ReceiverPosition - pRay2->AtTime(tMin).position).GetNormalized();
VistaVector3D vMin = (v3ReceiverPosition - pRayMin->at(idxMinDist).position).GetNormalized();
const float vProd1 = v1.Dot(vMin);
const float vProd2 = v2.Dot(vMin);
if (std::abs(vProd1 - vProd2) < dThreshold)
return {0, 1, 2};
if (vProd1 < vProd2)
return {0, 1};
return {1, 2};
}
#pragma endregion
#pragma region DOUBLING RESOLUTION - private
void CAdaptiveRayGrid::DoubleRayResolution()
{
if (NRays() <= 1)
return;
DoubleThetaResolution();
DoublePhiResolution();
RayMatrix newRayMatrix( NTheta(), RayVector( NPhi() ) );
m_vpNewRaysOfLastAdaptation.clear();
RayVector::const_iterator iteratorOldRays = Rays().cbegin();
for (int idxTheta = 0; idxTheta < m_vdThetaDeg.size(); idxTheta++)
{
const bool isNewTheta = (idxTheta % 2) == 1;
for (int idxPhi = 0; idxPhi < m_vdPhiDeg.size(); idxPhi++)
{
const bool isNewPhi = (idxPhi % 2) == 1;
RayPtr pRay;
if (IsPoleDirection(m_vdThetaDeg[idxTheta]))
pRay = Matrix()[idxTheta/2].front(); //Int-Division to convert to theta index of old matrix
else if (isNewTheta || isNewPhi)
{
pRay = std::make_shared<CRay>(m_v3SourcePos, m_vdThetaDeg[idxTheta], m_vdPhiDeg[idxPhi]);
m_vpNewRaysOfLastAdaptation.insert(pRay);
}
else
pRay = *iteratorOldRays;
newRayMatrix[idxTheta][idxPhi] = pRay;
if (!isNewTheta && !isNewPhi)
iteratorOldRays++;
}
}
SetRayMatrix(newRayMatrix);
}
void CAdaptiveRayGrid::DoubleThetaResolution()
{
if (m_vdThetaDeg.size() < 2)
return;
std::vector<double> newAngleVector;
newAngleVector.reserve(2 * m_vdThetaDeg.size() - 1);
newAngleVector.push_back(m_vdThetaDeg.front());
for (int idxAngle = 1; idxAngle < m_vdThetaDeg.size(); idxAngle++)
{
const double& angle1 = m_vdThetaDeg[idxAngle - 1];
const double& angle2 = m_vdThetaDeg[idxAngle];
newAngleVector.push_back((angle1 + angle2) / 2);
newAngleVector.push_back(angle2);
}
m_vdThetaDeg = newAngleVector;
m_dMaxDeltaTheta /= 2;
}
void CAdaptiveRayGrid::DoublePhiResolution()
{
if (m_vdPhiDeg.size() < 2)
return;
std::vector<double> newAngleVector;
newAngleVector.reserve(2 * m_vdPhiDeg.size() - 1 + IsCircular());
newAngleVector.push_back(m_vdPhiDeg.front());
for (int idxAngle = 1; idxAngle < m_vdPhiDeg.size(); idxAngle++)
{
const double& angle1 = m_vdPhiDeg[idxAngle - 1];
const double& angle2 = m_vdPhiDeg[idxAngle];
double newAngle = (angle1 + angle2) / 2;
if (angle2 < angle1)
newAngle = std::fmod(newAngle + 180, 360);
newAngleVector.push_back(newAngle);
newAngleVector.push_back(angle2);
}
if (IsCircular())
{
double newAngle = (m_vdPhiDeg.front() + m_vdPhiDeg.back()) / 2;
if (m_vdPhiDeg.front() < m_vdPhiDeg.back())
newAngle = std::fmod(newAngle + 180, 360);
newAngleVector.push_back(newAngle);
}
m_vdPhiDeg = newAngleVector;
m_dMaxDeltaPhi /= 2;
}
#pragma endregion
\ No newline at end of file
#include <ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/EigenrayEngine.h>
// ITA includes
#include "Worker.h"
#include <ITAPropagationPathSim/AtmosphericRayTracing/RayGrid.h>
// Vista includes
//#include <VistaInterProcComm/Concurrency/VistaThread.h>
// STD
using namespace ITAPropagationPathSim::AtmosphericRayTracing;
using namespace ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch;
std::vector<std::shared_ptr<CRay>> CEngine::Run(const ITAGeo::CStratifiedAtmosphere& atmosphere, const VistaVector3D& sourcePosition, const VistaVector3D& receiverPosition)
{
Simulation::Settings initialSettings = simulationSettings;
initialSettings.dIntegrationTimeStep *= 10;
CInitialWorker initialRayTracing(sourcePosition, receiverPosition, initialSettings, eigenraySettings.rayTracing);
std::vector<CRayGrid> initialRayGrids = initialRayTracing.Run(atmosphere);
RayVector eigenrays;
for (int reflectionOrder = 0; reflectionOrder < initialRayGrids.size(); reflectionOrder++)
{
CAdaptiveWorker worker(initialRayGrids[reflectionOrder], receiverPosition, simulationSettings, eigenraySettings, reflectionOrder);
eigenrays.push_back( worker.Run(atmosphere) );
}
return eigenrays;
}
\ No newline at end of file
/*
* ----------------------------------------------------------------
*
* ITA geometrical acoustics
* (c) Copyright Institute of Technical Acoustics (ITA)
* RWTH Aachen University, Germany, 2015-2019
*
* ----------------------------------------------------------------
* ____ __________ _______
* // / //__ ___/ // _ |
* // / // / // /_| |
* // / // / // ___ |
* //__/ //__/ //__/ |__|
*
* ----------------------------------------------------------------
*
*/
#ifndef IW_ITA_PROPAGATIONPATHSIM_ART_EIGENRAYSEARCH_WORKER
#define IW_ITA_PROPAGATIONPATHSIM_ART_EIGENRAYSEARCH_WORKER
#include <ITAPropagationPathSim/Definitions.h>
// Vista includes
#include <VistaBase/VistaVector3D.h>
// ITA includes
#include <ITAPropagationPathSim/AtmosphericRayTracing/Simulation/ARTEngine.h>
#include <ITAPropagationPathSim/AtmosphericRayTracing/Simulation/ARTSettings.h>
#include <ITAPropagationPathSim/AtmosphericRayTracing/Rays.h>
#include <ITAPropagationPathSim/AtmosphericRayTracing/RayGrid.h>
#include <ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/EigenraySettings.h>
#include <ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/AdaptiveRayGrid.h>
#include <ITAGeo/Atmosphere/StratifiedAtmosphere.h>
// STD
#include <vector>
#include <set>
#include <memory>
namespace ITAPropagationPathSim
{
namespace AtmosphericRayTracing
{
namespace EigenraySearch {
typedef std::shared_ptr<CRay> RayPtr;
typedef std::vector< RayPtr > RayVector;
class CReceiverData {
public:
const VistaVector3D& v3Position;
const int iReflectionOrder;
inline CReceiverData(const VistaVector3D& pos, const int reflOrder): v3Position(pos), iReflectionOrder(reflOrder) {};
inline CReceiverData& operator=(const CReceiverData& other)
{
if(this != &other)
*this = CReceiverData(other.v3Position, other.iReflectionOrder);
return *this;
};
};
typedef std::vector<CReceiverData> ReceiverDataVector;
class CEigenraySearchWatcher : public Simulation::IExternalWatcher
{
private:
ReceiverDataVector m_vReceiverData;
double m_dTMax;
int m_iMaxReflectionOrderRTAbort;
bool m_bAbortOnReceiverDistIncrease;
public:
CEigenraySearchWatcher(const ReceiverDataVector& vReceiverData, const double& tMax, const bool bAbortOnReceiverDistIncrease);
public:
void InitRayReceiverDistances(const std::set<RayPtr>& rays) const;
virtual bool AbortRequested(const std::shared_ptr<CRay> pRay) const;
//! Interface function to Simulation::Engine: Updates the minimum receiver distance of the ray after each processing step
virtual void ProcessRay(std::shared_ptr<CRay>) const;
//! Interface function to Simulation::Engine: Called after tracing a ray is finished.
virtual void FinalizeRay(std::shared_ptr<CRay>) const;
};
//! Basis for a worker which tracks the distance between rays and a receiver to find eigenrays
class CWorkerBase
{
protected:
std::shared_ptr<CEigenraySearchWatcher> m_pSimulationWatcher;
Simulation::CEngine m_simulationEngine;
RayTracingAbortSettings m_rayTracingAbortSettings;
VistaVector3D m_v3SourcePosition;
private:
VistaVector3D m_v3ReceiverPosition;
VistaVector3D m_v3MirroredReceiverPosition;
public:
CWorkerBase(const VistaVector3D& sourcePosition, const VistaVector3D& receiverPosition, const Simulation::Settings& simSettings, const RayTracingAbortSettings& abortSettings);
protected:
inline const VistaVector3D& ReceiverPosition() const { return m_v3ReceiverPosition; };
inline const VistaVector3D& MirroredReceiverPosition() const { return m_v3MirroredReceiverPosition; };
const VistaVector3D& VirtualReceiverPosition(const int reflectionOrder) const;
void RunRayTracing(const ITAGeo::CStratifiedAtmosphere& atmosphere, const std::set<RayPtr>& rays) const;
RayPtr FindMinimumDistanceRay(const std::set<RayPtr>& rays, const int reflectionOrder) const;
};
//! Does a rough search for potential eigenrays directions for multiple reflection orders
class CInitialWorker : public CWorkerBase
{
private:
RayVector m_vpMinDistanceRays; //!< Vector containing the current ray of minimum receiver distance for each reflection order
public:
CInitialWorker(const VistaVector3D& sourcePosition, const VistaVector3D& receiverPosition, const Simulation::Settings& simSettings, const RayTracingAbortSettings& abortSettings);
std::vector<CRayGrid> Run(const ITAGeo::CStratifiedAtmosphere& atmosphere);
private:
CRayGrid InitRayGrid(const ITAGeo::CStratifiedAtmosphere& atmosphere);
void FindMinimumDistanceRays(const std::set<RayPtr>& rays);
std::vector<CRayGrid> FinalizeResult(const CRayGrid& initialRayGrid);
};
//! Based on a result of CInitialWorker, this worker adaptively reduces the directional limits while increasing the resolution to find an eigenray of a specific reflection order
class CAdaptiveWorker : public CWorkerBase
{
private:
CAdaptiveRayGrid m_adaptiveRayGrid;
const int m_iActiveReflexionOrder;