From fbfb1f15c57ef7d037db2bf6eb11bb7a1911a780 Mon Sep 17 00:00:00 2001
From: Maximilian Schnabel
<7090722+MaximilianSchnabel@users.noreply.github.com>
Date: Thu, 7 Jul 2022 14:15:43 +0200
Subject: [PATCH] GeometricalAcoustics/PathSim: implement CCombinedWorker to
combine calculation of direct and reflective path
---
.../EigenraySearch/EigenrayEngine.h | 8 +-
.../EigenraySearch/EigenrayEngine.cpp | 51 +++-
.../EigenraySearch/Worker.cpp | 239 +++++++++++++++++-
.../EigenraySearch/Worker.h | 32 +++
4 files changed, 323 insertions(+), 7 deletions(-)
diff --git a/include/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/EigenrayEngine.h b/include/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/EigenrayEngine.h
index d2e4e95..63f6bef 100644
--- a/include/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/EigenrayEngine.h
+++ b/include/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/EigenrayEngine.h
@@ -51,7 +51,13 @@ namespace ITAPropagationPathSim
inline CEngine( ) { };
std::vector<std::shared_ptr<CRay>> Run( const ITAGeo::CStratifiedAtmosphere& atmosphere, const VistaVector3D& sourcePosition,
- const VistaVector3D& receiverPosition ) const;
+ const VistaVector3D& receiverPosition, bool combineReflCalc = false) const;
+
+ std::vector<std::shared_ptr<CRay>> Run( const ITAGeo::CStratifiedAtmosphere& atmosphere, const VistaVector3D& sourcePosition,
+ const VistaVector3D& receiverPosition, const VistaVector3D& initDirectDirection,
+ const VistaVector3D& initReflectionDirection,
+ double phiResolution, double thetaResolution,
+ bool combineReflCalc ) const;
};
} // namespace EigenraySearch
} // namespace AtmosphericRayTracing
diff --git a/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/EigenrayEngine.cpp b/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/EigenrayEngine.cpp
index 10f8444..4a00a78 100644
--- a/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/EigenrayEngine.cpp
+++ b/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/EigenrayEngine.cpp
@@ -15,7 +15,7 @@ using namespace ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch;
std::vector<std::shared_ptr<CRay>> CEngine::Run( const ITAGeo::CStratifiedAtmosphere& atmosphere, const VistaVector3D& sourcePosition,
- const VistaVector3D& receiverPosition ) const
+ const VistaVector3D& receiverPosition, bool combineReflCalc ) const
{
Simulation::Settings initialSettings = simulationSettings;
initialSettings.dIntegrationTimeStep *= 10;
@@ -23,10 +23,51 @@ std::vector<std::shared_ptr<CRay>> CEngine::Run( const ITAGeo::CStratifiedAtmosp
std::vector<CRayGrid> initialRayGrids = initialRayTracing.Run( atmosphere );
RayVector eigenrays;
- for( int reflectionOrder = 0; reflectionOrder < initialRayGrids.size( ); reflectionOrder++ )
+ if( combineReflCalc )
{
- CAdaptiveWorker worker( initialRayGrids[reflectionOrder], receiverPosition, simulationSettings, eigenraySettings, reflectionOrder );
- eigenrays.push_back( worker.Run( atmosphere ) );
+ CCombinedAdaptiveWorker worker( initialRayGrids[0], receiverPosition, simulationSettings, eigenraySettings, 2 );
+ eigenrays = worker.Run( atmosphere );
+ }
+ else
+ {
+ for( int reflectionOrder = 0; reflectionOrder < initialRayGrids.size( ); reflectionOrder++ )
+ {
+ CAdaptiveWorker worker( initialRayGrids[reflectionOrder], receiverPosition, simulationSettings, eigenraySettings, reflectionOrder );
+ eigenrays.push_back( worker.Run( atmosphere ) );
+ }
+ }
+ return eigenrays;
+}
+
+std::vector<std::shared_ptr<CRay>> ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CEngine::Run( const ITAGeo::CStratifiedAtmosphere& atmosphere,
+ const VistaVector3D& sourcePosition,
+ const VistaVector3D& receiverPosition,
+ const VistaVector3D& initDirectDirection,
+ const VistaVector3D& initReflectionDirection,
+ double phiResolution,
+ double thetaResolution,
+ bool combineReflCalc) const
+{
+ Simulation::Settings initialSettings = simulationSettings;
+ initialSettings.dIntegrationTimeStep *= 10;
+ CRayGrid initDirGrid = CPreinitializedRayGrid( sourcePosition, initDirectDirection, phiResolution, thetaResolution );
+
+ RayVector eigenrays;
+ if( combineReflCalc )
+ {
+ CCombinedAdaptiveWorker worker( initDirGrid, receiverPosition, simulationSettings, eigenraySettings, 2 );
+ eigenrays = worker.Run( atmosphere );
+ }
+ else
+ {
+ CRayGrid initReflGrid = CPreinitializedRayGrid( sourcePosition, initReflectionDirection, phiResolution, thetaResolution );
+ std::vector<CRayGrid> initialRayGrids = { initDirGrid, initReflGrid };
+
+ 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
+}
diff --git a/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/Worker.cpp b/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/Worker.cpp
index 5091a86..7d414cc 100644
--- a/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/Worker.cpp
+++ b/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/Worker.cpp
@@ -299,4 +299,241 @@ double EigenraySearch::CAdaptiveWorker::CalculateEnergyProportion( const double&
const double dDenominator = A * ( 1 + dRelVToC ) * std::sqrt( 1 + 2 * dRelVToC + vSquared / ( c * c ) );
return c / dDenominator;
}
-#pragma endregion
\ No newline at end of file
+#pragma endregion
+
+ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CCombinedAdaptiveWorker::CCombinedAdaptiveWorker( const CRayGrid& rayGrid,
+ const VistaVector3D& receiverPosition,
+ const Simulation::Settings& simSettings,
+ const Settings& eigenraySettings,
+ const int reflectionOrders )
+ : CWorkerBase( rayGrid.SourcePosition( ), receiverPosition, simSettings, eigenraySettings.rayTracing )
+ , m_reflectionOrders( reflectionOrders )
+ , m_rayAdaptationSettings( eigenraySettings.rayAdaptation )
+{
+ m_adaptiveRayGrids[0] = std::make_unique<CAdaptiveRayGrid>( rayGrid );
+ m_adaptations[0] = 0;
+
+ const float sourceReceiverDistance = ( VirtualReceiverPosition(0) - rayGrid.SourcePosition( ) ).GetLength( );
+ m_fReceiverRadius = (float)tan( m_rayAdaptationSettings.accuracy.maxSourceReceiverAngle * M_PI / 180.0 ) * sourceReceiverDistance;
+ m_fReceiverRadius = fmin( (float)m_rayAdaptationSettings.accuracy.maxReceiverRadius, m_fReceiverRadius );
+
+ ReceiverDataVector receiverData;
+ for( int reflOrder = 0; reflOrder < reflectionOrders ; ++reflOrder)
+ receiverData.push_back(CReceiverData( VirtualReceiverPosition( reflOrder ), reflOrder ));
+
+ m_pSimulationWatcher =
+ std::make_shared<CEigenraySearchWatcher>( receiverData, m_rayTracingAbortSettings.maxTime, m_rayTracingAbortSettings.bAbortOnReceiverDistIncrease );
+ m_simulationEngine.pExternalWatcher = m_pSimulationWatcher;
+}
+
+std::vector<EigenraySearch::RayPtr> ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CCombinedAdaptiveWorker::Run( const ITAGeo::CStratifiedAtmosphere& atmosphere )
+{
+ RunRayTracing( atmosphere, m_adaptiveRayGrids[0]->UniqueRays( ) );
+ m_pMinDistanceRays = FindMinimumDistanceRays( m_reflectionOrders );
+ bool keepWorking = true;
+ while( keepWorking )
+ {
+ SplitIfNecessary( );
+ keepWorking = false;
+ for( auto& [reflectionOrder, pGrid]: m_adaptiveRayGrids )
+ {
+ RayPtr& ray = m_pMinDistanceRays[reflectionOrder];
+ bool iterationWorks = !EigenrayAccuracyReached( ray, reflectionOrder ) && !AbortCriterionReached( reflectionOrder );
+ keepWorking |= iterationWorks;
+
+ if( !iterationWorks )
+ continue;
+
+ if( m_rayAdaptationSettings.advancedRayZooming.bActive )
+ pGrid->ZoomIntoRay(
+ m_pMinDistanceRays[reflectionOrder],
+ VirtualReceiverPosition( reflectionOrder ),
+ m_rayAdaptationSettings.advancedRayZooming.threshold
+ );
+ else
+ pGrid->ZoomIntoRay( m_pMinDistanceRays[reflectionOrder] );
+ RunRayTracing( atmosphere, pGrid->NewRaysOfLastAdaptation( ) );
+ m_adaptations[reflectionOrder]++;
+ }
+
+ m_pMinDistanceRays = FindMinimumDistanceRays( m_reflectionOrders );
+ }
+
+ PostProcessEigenray( atmosphere );
+ return m_pMinDistanceRays;
+}
+
+std::vector<EigenraySearch::RayPtr> ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CCombinedAdaptiveWorker::FindMinimumDistanceRays( const int underReflectionOrder ) const
+{
+ struct ClosestRayData
+ {
+ float dMin = FLT_MAX;
+ RayPtr pMinRay = nullptr;
+ };
+ std::vector<RayPtr> pMinRays = { };
+ std::vector<ClosestRayData> rayData(underReflectionOrder);
+
+ for( int reflOrder = 0; reflOrder < underReflectionOrder; ++reflOrder )
+ {
+ bool wasSplit = m_adaptiveRayGrids.count( reflOrder ) == 1;
+ const std::unique_ptr<CAdaptiveRayGrid>* usedGrid = nullptr;
+ if( wasSplit )
+ usedGrid = &m_adaptiveRayGrids.at(reflOrder);
+ else
+ usedGrid = &m_adaptiveRayGrids.at(0);
+
+ for( const RayPtr& pRay: (*usedGrid)->UniqueRays() )
+ {
+ const VistaVector3D& receiverPos = VirtualReceiverPosition( reflOrder );
+ const CRayReceiverData* distanceData = pRay->ReceiverDistanceData( receiverPos );
+ if( distanceData && distanceData->distance < rayData[reflOrder].dMin )
+ {
+ rayData[reflOrder].dMin = distanceData->distance;
+ rayData[reflOrder].pMinRay = pRay;
+ }
+ }
+ }
+
+ for( ClosestRayData& data: rayData)
+ pMinRays.push_back( data.pMinRay );
+
+ return pMinRays;
+}
+
+void ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CCombinedAdaptiveWorker::SplitIfNecessary( )
+{
+ for( int reflOrder = 1; reflOrder < m_reflectionOrders; ++reflOrder )
+ {
+ //TODO:: check if actual ray comparison is needed instead of pointer obj comparison
+ if( m_pMinDistanceRays[0] != m_pMinDistanceRays[reflOrder] )
+ {
+ if( m_adaptiveRayGrids.count( reflOrder ) != 1 )
+ {
+ m_adaptations[reflOrder] = 0;
+ m_adaptiveRayGrids[reflOrder] = std::make_unique<CAdaptiveRayGrid>( *m_adaptiveRayGrids[0] );
+ }
+ }
+ }
+}
+bool ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CCombinedAdaptiveWorker::EigenrayAccuracyReached( RayPtr& ray, int reflOrder )
+{
+ const CRayReceiverData* receiverData = ray->ReceiverDistanceData( VirtualReceiverPosition( reflOrder ) );
+ return ( receiverData && receiverData->distance <= m_fReceiverRadius );
+}
+
+bool ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CCombinedAdaptiveWorker::AbortCriterionReached( int reflOrder )
+{
+ return m_adaptations[reflOrder] > m_rayAdaptationSettings.abort.maxNAdaptations ||
+ m_adaptiveRayGrids[reflOrder]->MaxAngularResolution( ) < m_rayAdaptationSettings.abort.minAngleResolutionDeg;
+}
+
+void ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CCombinedAdaptiveWorker::PostProcessEigenray( const ITAGeo::CStratifiedAtmosphere& atmosphere )
+{
+ SetEigenrayStatus( );
+ for( auto& [reflOrder, pGrid]: m_adaptiveRayGrids )
+ {
+ SetEigenrayEndPoint( m_pMinDistanceRays[reflOrder], reflOrder );
+ CalculateEigenraySpreadingLoss( atmosphere, m_pMinDistanceRays[reflOrder], reflOrder );
+ }
+}
+
+
+void ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CCombinedAdaptiveWorker::SetEigenrayStatus( )
+{
+ for( int reflOrder = 0; reflOrder < m_pMinDistanceRays.size(); ++reflOrder )
+ {
+ RayPtr& ray = m_pMinDistanceRays[reflOrder];
+ bool bAccuracyReached = EigenrayAccuracyReached( ray, reflOrder );
+ if( !bAccuracyReached )
+ std::cout << "WARNING: EigenraySearch::Engine: Could not find eigenray with proper accuracy. Abort criterion reached." << std::endl;
+
+ ray->MarkAsEigenray( bAccuracyReached );
+ }
+}
+
+void ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CCombinedAdaptiveWorker::SetEigenrayEndPoint( RayPtr& ray, int reflOrder )
+{
+ const CRayReceiverData* receiverData = ray->ReceiverDistanceData( VirtualReceiverPosition(reflOrder) );
+ const int idxBeforeMin = receiverData->idxMinDist;
+ const VistaVector3D& rMin = receiverData->posMinDist;
+
+ if( ray->NumPoints( ) < idxBeforeMin )
+ ITA_EXCEPT_INVALID_PARAMETER( "Index for minimum distance out of bounds!" );
+
+ if( ray->NumPoints( ) == idxBeforeMin + 1 ) // Min dist point is endpoint
+ {
+ std::cout << "WARNING: EigenraySearch::Engine: Ray end point is closest point to receiver. This usually happens if maximum propagation time is too low."
+ << std::endl;
+ return;
+ }
+
+ // Interpolate to new point of minimum
+ const VistaVector3D& r1 = ray->at( idxBeforeMin ).position;
+ const VistaVector3D& r2 = ray->at( idxBeforeMin + 1 ).position;
+ const float pathFraction = ( rMin - r1 ).GetLength( ) / ( r2 - r1 ).GetLength( );
+ CRayElement newRayElement = ray->at( idxBeforeMin ).Interpolate( ray->at( idxBeforeMin + 1 ), pathFraction );
+
+ // Remove all points before and add new point
+ ray->resize( idxBeforeMin + 2 );
+ ray->back( ) = newRayElement;
+}
+
+void ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CCombinedAdaptiveWorker::CalculateEigenraySpreadingLoss(
+ const ITAGeo::CStratifiedAtmosphere& atmosphere, RayPtr& ray, int reflOrder )
+{
+ const double zRef = std::abs( ray->SourcePoint( )[Vista::Z] );
+ const VistaVector3D& nRef = ray->InitialDirection( );
+
+ const double& time = ray->back( ).timeStamp;
+ const double zReceiver = std::abs( ray->back( ).position[Vista::Z] );
+ const VistaVector3D& nReceiver = ray->back( ).wavefrontNormal;
+
+ const double dTheta = 180.0 / M_PI * std::acos( nRef[Vista::Z] );
+ const bool bCloseToPole = dTheta < 1.0 || dTheta > 179.0;
+ if( bCloseToPole || m_adaptiveRayGrids.at(reflOrder)->MaxDeltaTheta( ) > m_rayAdaptationSettings.accuracy.maxAngleForGeomSpreading || !m_adaptiveRayGrids.at(reflOrder)->Is2D( ) )
+ {
+ double dDeltaPhi = m_rayAdaptationSettings.accuracy.maxAngleForGeomSpreading;
+
+ // Lower phi resolution if close pole to avoid numerical cancellation
+ const bool bIsPole = std::abs( nRef[Vista::X] ) < FLT_EPSILON && std::abs( nRef[Vista::Y] ) < FLT_EPSILON;
+ if( bIsPole )
+ dDeltaPhi = 5.0;
+ else if( bCloseToPole )
+ dDeltaPhi *= 10;
+
+ m_adaptiveRayGrids[reflOrder]->ZoomIntoRay( ray, m_rayAdaptationSettings.accuracy.maxAngleForGeomSpreading, dDeltaPhi );
+ auto tmpSimulationEngine = Simulation::CEngine( std::make_shared<Simulation::CAbortAtMaxTime>( time ) );
+ tmpSimulationEngine.settings = m_simulationEngine.settings;
+ tmpSimulationEngine.Run( atmosphere, m_adaptiveRayGrids[reflOrder]->NewRaysOfLastAdaptation( ) );
+ }
+ else
+ m_adaptiveRayGrids[reflOrder]->SetToSurroundingGrid( ray );
+
+ const double surfaceReceiver = m_adaptiveRayGrids[reflOrder]->WavefrontSurface( time );
+ const double cReceiver = atmosphere.SpeedOfSound( zReceiver );
+ const VistaVector3D vReceiver = atmosphere.WindVector( zReceiver );
+ const double dEnergyProportionReceiver = CalculateEnergyProportion( surfaceReceiver, cReceiver, nReceiver, vReceiver );
+
+ if( std::isinf( dEnergyProportionReceiver ) )
+ {
+ std::cout << "WARNING: EigenraySearch::Engine: Numerical cancellation during spreading loss calculation. Spreading loss is unset (= -1)." << std::endl;
+ return;
+ }
+
+ const double surfaceRef = m_adaptiveRayGrids[reflOrder]->WavefrontSurfaceReference( );
+ const double cRef = atmosphere.SpeedOfSound( zRef );
+ const VistaVector3D vRef = atmosphere.WindVector( zRef );
+ const double dEnergyProportionRef = CalculateEnergyProportion( surfaceRef, cRef, nRef, vRef );
+
+ ray->SetSpreadingLoss( std::sqrt( dEnergyProportionReceiver / dEnergyProportionRef ) );
+}
+
+double ITAPropagationPathSim::AtmosphericRayTracing::EigenraySearch::CCombinedAdaptiveWorker::CalculateEnergyProportion( const double& A, const double& c,
+ const VistaVector3D& vecN,
+ const VistaVector3D& vecV )
+{
+ const double dRelVToC = (double)( vecN * vecV ) / c;
+ const double vSquared = vecV.GetLengthSquared( );
+ const double dDenominator = A * ( 1 + dRelVToC ) * std::sqrt( 1 + 2 * dRelVToC + vSquared / ( c * c ) );
+ return c / dDenominator;
+}
diff --git a/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/Worker.h b/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/Worker.h
index 80ef570..f8d5bd3 100644
--- a/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/Worker.h
+++ b/src/ITAPropagationPathSim/AtmosphericRayTracing/EigenraySearch/Worker.h
@@ -158,6 +158,38 @@ namespace ITAPropagationPathSim
void CalculateEigenraySpreadingLoss( const ITAGeo::CStratifiedAtmosphere& atmosphere );
double CalculateEnergyProportion( const double& A, const double& c, const VistaVector3D& vecN, const VistaVector3D& vecV );
};
+
+ class CCombinedAdaptiveWorker : public CWorkerBase
+ {
+ private:
+ const int m_reflectionOrders;
+ RayAdaptationSettings m_rayAdaptationSettings;
+ float m_fReceiverRadius;
+
+ std::vector<RayPtr> m_pMinDistanceRays;
+ //int m_nAdaptations = 0;
+
+ std::map<int, std::unique_ptr<CAdaptiveRayGrid>> m_adaptiveRayGrids;
+ std::map<int, int> m_adaptations;
+
+ public:
+ CCombinedAdaptiveWorker( const CRayGrid& rayGrid, const VistaVector3D& receiverPosition, const Simulation::Settings& simSettings,
+ const Settings& eigenraySettings, const int reflectionOrders );
+
+ std::vector<RayPtr> Run( const ITAGeo::CStratifiedAtmosphere& atmosphere );
+
+ private:
+ std::vector<RayPtr> FindMinimumDistanceRays( const int underReflectionOrder ) const;
+
+ void SplitIfNecessary( );
+ bool EigenrayAccuracyReached( RayPtr& ray, int reflOrder );
+ bool AbortCriterionReached( int reflOrder );
+ void PostProcessEigenray( const ITAGeo::CStratifiedAtmosphere& atmosphere );
+ void SetEigenrayStatus( );
+ void SetEigenrayEndPoint( RayPtr& ray, int reflOrder );
+ void CalculateEigenraySpreadingLoss( const ITAGeo::CStratifiedAtmosphere& atmosphere, RayPtr& ray, int reflOrder );
+ double CalculateEnergyProportion( const double& A, const double& c, const VistaVector3D& vecN, const VistaVector3D& vecV );
+ };
} // namespace EigenraySearch
} // namespace AtmosphericRayTracing
} // namespace ITAPropagationPathSim
--
GitLab