...
 
Commits (2)
  • Philipp Schäfer's avatar
    ART - RayGrid · 0a8f0e2c
    Philipp Schäfer authored
    - added bool to (de)activate sorting of angle vectors
    - fixed bug where ray angles where falsely sorted on CreateCopyWithNewRays()
    0a8f0e2c
  • Philipp Schäfer's avatar
    ART - EigenrayEngineTest · b61066e6
    Philipp Schäfer authored
    - added test for differenz azimuths between source and receiver
    - matlab scripts now also work with new classes from ARTMatlab app
    - split matlab test into two: one only testing the cpp data, the other comparing to matlab simulation
    b61066e6
......@@ -67,9 +67,9 @@ namespace ITAPropagationPathSim
//! 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.
* 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);
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.
......
......@@ -23,11 +23,14 @@ CRayGrid::CRayGrid(const RayMatrix& rayMatrix, const std::vector<double>& thetaD
m_v3SourcePos = m_vpRays.front()->SourcePoint();
}
CRayGrid::CRayGrid(const VistaVector3D& sourcePos, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg, const bool circularPhi /*= false*/)
CRayGrid::CRayGrid(const VistaVector3D& sourcePos, const std::vector<double>& thetaDeg, const std::vector<double>& phiDeg, const bool circularPhi /*= false*/, const bool sortAngleVectors /*= true*/)
: m_v3SourcePos(sourcePos), m_vdThetaDeg(thetaDeg), m_vdPhiDeg(phiDeg), m_bCircularPhi(circularPhi)
{
std::sort(m_vdPhiDeg.begin(), m_vdPhiDeg.end());
std::sort(m_vdThetaDeg.begin(), m_vdThetaDeg.end());
if (sortAngleVectors)
{
std::sort(m_vdPhiDeg.begin(), m_vdPhiDeg.end());
std::sort(m_vdThetaDeg.begin(), m_vdThetaDeg.end());
}
if (m_vdThetaDeg.front() < 0.0 || m_vdThetaDeg.back() > 180.0)
ITA_EXCEPT_INVALID_PARAMETER("Values for elevation angle theta are out of bounds [0, 180]");
......@@ -293,7 +296,7 @@ void CRayGrid::SetToSurroundingGrid(const std::shared_ptr<CRay> pRay)
#pragma region COPY - public
CRayGrid CRayGrid::CopyWithNewRays() const
{
return CRayGrid(m_v3SourcePos, m_vdThetaDeg, m_vdPhiDeg, m_bCircularPhi);
return CRayGrid(m_v3SourcePos, m_vdThetaDeg, m_vdPhiDeg, m_bCircularPhi, false);
}
#pragma endregion
......
......@@ -41,6 +41,8 @@
#include <stdio.h>
#include <math.h>
#include <ITAFileSystemUtils.h>
#include <omp.h>
using namespace std;
......@@ -71,28 +73,39 @@ void runTest(const CStratifiedAtmosphere& atmosphere, const VistaVector3D& sourc
cout << "Starting export..." << endl;
CPropagationPathList raysAsPropPathList = ToPropagationPath(eigenrays);
JSON::Export(raysAsPropPathList, filename + ".json");
makeDirectory("output");
std::string filepath = combinePath("output", filename + ".json");
JSON::Export(raysAsPropPathList, filepath);
cout << "Finished" << endl << endl;
}
void runTestReceiverNearGround(const CStratifiedAtmosphere& atmosphere, const int& sourceAngleDeg, const string& atmosphereSuffix)
void runTestReceiverNearGround(const CStratifiedAtmosphere& atmosphere, const int sourceElevationDeg, const int sourceAzimuthnDeg, const string& atmosphereSuffix)
{
VistaVector3D receiverPosition = VistaVector3D(0, 0, 1.8);
const double sourceReceiverDistance = 2000;
const double sourceX = sourceReceiverDistance * std::sin(sourceAngleDeg * M_PI / 180.0) + receiverPosition[Vista::X];
const double sourceY = +receiverPosition[Vista::Y];
const double sourceZ = sourceReceiverDistance * std::cos(sourceAngleDeg * M_PI / 180.0) + receiverPosition[Vista::Z];
const double sourceX = sourceReceiverDistance * std::sin(sourceElevationDeg * M_PI / 180.0) * std::cos(sourceAzimuthnDeg * M_PI / 180.0) + receiverPosition[Vista::X];
const double sourceY = sourceReceiverDistance * std::sin(sourceElevationDeg * M_PI / 180.0) * std::sin(sourceAzimuthnDeg * M_PI / 180.0) + receiverPosition[Vista::Y];
const double sourceZ = sourceReceiverDistance * std::cos(sourceElevationDeg * M_PI / 180.0) + receiverPosition[Vista::Z];
const VistaVector3D sourcePosition = VistaVector3D(sourceX, sourceY, sourceZ);
const string filename = "ReceiverNearGround_" + atmosphereSuffix +"_SourceAngle" + to_string(sourceAngleDeg);
const string filename = "ReceiverNearGround_" + atmosphereSuffix +"_SourceElevation" + to_string(sourceElevationDeg) + "_SourceAzimuth" + to_string(sourceAzimuthnDeg);
runTest(atmosphere, sourcePosition, receiverPosition, filename);
}
void TestReceiverNearGround(const CStratifiedAtmosphere& atmosphere, const string& atmosphereSuffix)
void TestReceiverNearGroundSourceElevation(const CStratifiedAtmosphere& atmosphere, const string& atmosphereSuffix)
{
const int phi = 0;
for (int sourceElevation = -90; sourceElevation <= 90; sourceElevation += 45)
runTestReceiverNearGround(atmosphere, sourceElevation, phi, atmosphereSuffix);
}
void TestReceiverNearGroundSourceAzimuth(const CStratifiedAtmosphere& atmosphere, const string& atmosphereSuffix)
{
for (int angle = -90; angle <= 90; angle += 45)
runTestReceiverNearGround(atmosphere, angle, atmosphereSuffix);
const int theta = 30;
for (int sourceAzimuth = 0; sourceAzimuth < 360; sourceAzimuth += 10)
runTestReceiverNearGround(atmosphere, theta, sourceAzimuth, atmosphereSuffix);
}
CStratifiedAtmosphere GetHomogeneousAtmosphere()
......@@ -111,16 +124,29 @@ CStratifiedAtmosphere GetInhomogeneousAtmosphere()
return CStratifiedAtmosphere(tempProfile, windProfile, humidProfile);
}
void TestSourceReceiverAzimuthAbove324Deg()
{
const CStratifiedAtmosphere atmosphere = GetInhomogeneousAtmosphere();
const VistaVector3D receiverPosition = VistaVector3D(0, 0, 1.8);
const VistaVector3D sourcePosition = VistaVector3D(-1000, 200, 600);
const string filename = "SourceReceiverAzimuthAbove324Deg";
runTest(atmosphere, sourcePosition, receiverPosition, filename);
}
int main(int iNumInArgs, char* pcInArgs[])
{
//Disable multi-threading for debugging purposes
omp_set_num_threads(1);
CStratifiedAtmosphere homAtmosphere = GetHomogeneousAtmosphere();
CStratifiedAtmosphere inhomAtmosphere = GetInhomogeneousAtmosphere();
const CStratifiedAtmosphere homAtmosphere = GetHomogeneousAtmosphere();
const CStratifiedAtmosphere inhomAtmosphere = GetInhomogeneousAtmosphere();
TestReceiverNearGroundSourceElevation(homAtmosphere, "Homogeneous");
TestReceiverNearGroundSourceElevation(inhomAtmosphere, "Inhomogeneous");
TestReceiverNearGroundSourceAzimuth(inhomAtmosphere, "Inhomogeneous");
TestReceiverNearGround(homAtmosphere, "Homogeneous");
TestReceiverNearGround(inhomAtmosphere, "Inhomogeneous");
TestSourceReceiverAzimuthAbove324Deg();
}
function EigenrayEngineTest()
%% Define Atmospheres
jsonFolder = fullfile( fileparts(mfilename('fullpath')), 'output' );
inhomAtmos = itaAtmosphere;
%% Special case azimuth > 324
cppFilename = 'SourceReceiverAzimuthAbove324Deg.json';
checkCppEigenrays( fullfile(jsonFolder, cppFilename) );
%Settings
%(These are the standard settings, itaAtmosphere is already initialized with these values)
inhomAtmos.windProfile = WindProfile.Log; %Enum: Zero, Const, Log
inhomAtmos.tempProfile = TempProfile.ISA; %Enum: Const, ISA
inhomAtmos.z0 = 0.1; %Surface Roughness for Log Wind Profile [m]
inhomAtmos.kappa = 0.4; %Karman constant for Log Wind Profile []
inhomAtmos.vStar = 0.6; %Friction velocity for Log Wind Profile [m/s]
inhomAtmos.vDirConst = [1 0 0]; %Normal in wind direction []
inhomAtmos.R = 1.4; %Air Gas Constant, used to calculate c
inhomAtmos.gamma = 402.8/1.4; %Ratio of Specific Heats, used to calculate c
homAtmos = inhomAtmos;
homAtmos.windProfile = WindProfile.Zero;
homAtmos.tempProfile = TempProfile.Const;
homAtmos.TConst = 293.15;
homAtmos.constStaticPressure = 101325;
%% Define Ray Tracer
art = itaART;
art.source = [0 0 50]; %Can also be a single itaCoordinates
art.receiver = [0 0 1.8]; %Can also be a single itaCoordinates
art.dt = 0.1; %Integration variable [s]
art.tMax = 15; %Maximum time of tracing [s]
art.bAdaptiveDt = true;
art.maxErrorV = 0.015;
art.uncriticalErrorV = 0.005;
%art: maximum level for dt adaptation = 32
art.bCopyPeriodicRefl = true;
art.bUseSlowness = true;
art.bUseWaitbar = false;
%% Eigenray search properties
art.bLimitEigenraySearchDirections = false; %If true, directions in which an eigenray is looked for will be limited
art.maxReceiverRadius = 1; %Maximum receiver radius used to calculate receiverRadius
art.maxSourceReceiverAngle = 1; %Maximum allowed angle between direct connection and eigenray used to calculate receiverRadius
art.maxEigenraySearchIterations = 30; %Abort criterion for eigenray search: maximum number of adaptations.
art.minEigenraySearchAngle = 0.001; %Abort criterion for eigenray search: minimum angle between rays.
art.rayZoomThresh = 1.0; %Threshold used for improved ray zooming. If this is undershot, the algorithm does not make a decision.
art.angleForGeomSpreading = 0.01; %Angle between rays used to calculate the geometrical spreading
art.maxReflOrder = 1;
%% Test Matlab vs Cpp
% art.atmosphere = homAtmos;
% runTestReceiverNearGround(art, 'Homogeneous');
% art.atmosphere = inhomAtmos;
% runTestReceiverNearGround(art, 'Inhomogeneous');
%% Check C++ data only
checkResultsReceiverNearGround('Homogeneous');
checkResultsReceiverNearGround('Inhomogeneous');
checkResultsReceiverNearGroundSourceElevation(jsonFolder, 'Homogeneous');
checkResultsReceiverNearGroundSourceElevation(jsonFolder, 'Inhomogeneous');
checkResultsReceiverNearGroundSourceAzimuth(jsonFolder, 'Inhomogeneous');
function runTestReceiverNearGround(art, atmosphereSuffix)
sourceReceiverDistance = 2000;
for sourceRecAngle = -90:45:90
art.source.x = sourceReceiverDistance*sind(sourceRecAngle) + art.receiver.x;
art.source.y = art.receiver.y;
art.source.z = sourceReceiverDistance*cosd(sourceRecAngle) + art.receiver.z;
cppFilename = ['ReceiverNearGround_' atmosphereSuffix '_SourceAngle' num2str(sourceRecAngle) '.json'];
runMatlabVsCppEigenrayTest(art, cppFilename, fileparts(cppFilename), true);
function checkResultsReceiverNearGroundSourceElevation(jsonFolder, atmosphereSuffix)
sourceAzimuth = 0;
for sourceElevation = -90:45:90
cppFilename = ['ReceiverNearGround_' atmosphereSuffix '_SourceElevation' num2str(sourceElevation) '_SourceAzimuth' num2str(sourceAzimuth) '.json'];
checkCppEigenrays( fullfile(jsonFolder, cppFilename) );
end
function checkResultsReceiverNearGround(atmosphereSuffix)
for sourceRecAngle = -90:45:90
cppFilename = ['ReceiverNearGround_' atmosphereSuffix '_SourceAngle' num2str(sourceRecAngle) '.json'];
checkCppEigenrays(cppFilename, fileparts(cppFilename));
function checkResultsReceiverNearGroundSourceAzimuth(jsonFolder, atmosphereSuffix)
sourceElevation = 30;
for sourceAzimuth = 0:10:350
cppFilename = ['ReceiverNearGround_' atmosphereSuffix '_SourceElevation' num2str(sourceElevation) '_SourceAzimuth' num2str(sourceAzimuth) '.json'];
checkCppEigenrays( fullfile(jsonFolder, cppFilename) );
end
function EigenrayEngineTest_CppVsMatlab()
%% Define Atmospheres
inhomAtmos = itaAtmosphere;
%Settings
%(These are the standard settings, itaAtmosphere is already initialized with these values)
inhomAtmos.windProfile = WindProfile.Log; %Enum: Zero, Const, Log
inhomAtmos.tempProfile = TempProfile.ISA; %Enum: Const, ISA
inhomAtmos.z0 = 0.1; %Surface Roughness for Log Wind Profile [m]
inhomAtmos.kappa = 0.4; %Karman constant for Log Wind Profile []
inhomAtmos.vStar = 0.6; %Friction velocity for Log Wind Profile [m/s]
inhomAtmos.vDirConst = [1 0 0]; %Normal in wind direction []
inhomAtmos.R = 1.4; %Air Gas Constant, used to calculate c
inhomAtmos.gamma = 402.8/1.4; %Ratio of Specific Heats, used to calculate c
homAtmos = inhomAtmos;
homAtmos.windProfile = WindProfile.Zero;
homAtmos.tempProfile = TempProfile.Const;
homAtmos.TConst = 293.15;
homAtmos.constStaticPressure = 101325;
%% Define Ray Tracer
art = itaART;
art.source = [0 0 50]; %Can also be a single itaCoordinates
art.receiver = [0 0 1.8]; %Can also be a single itaCoordinates
art.dt = 0.1; %Integration variable [s]
art.tMax = 15; %Maximum time of tracing [s]
art.bAdaptiveDt = true;
art.maxErrorV = 0.015;
art.uncriticalErrorV = 0.005;
%art: maximum level for dt adaptation = 32
art.bCopyPeriodicRefl = true;
art.bUseSlowness = true;
art.bUseWaitbar = false;
%% Eigenray search properties
art.bLimitEigenraySearchDirections = false; %If true, directions in which an eigenray is looked for will be limited
art.maxReceiverRadius = 1; %Maximum receiver radius used to calculate receiverRadius
art.maxSourceReceiverAngle = 1; %Maximum allowed angle between direct connection and eigenray used to calculate receiverRadius
art.maxEigenraySearchIterations = 30; %Abort criterion for eigenray search: maximum number of adaptations.
art.minEigenraySearchAngle = 0.001; %Abort criterion for eigenray search: minimum angle between rays.
art.rayZoomThresh = 1.0; %Threshold used for improved ray zooming. If this is undershot, the algorithm does not make a decision.
art.angleForGeomSpreading = 0.01; %Angle between rays used to calculate the geometrical spreading
art.maxReflOrder = 1;
%% Test Matlab vs Cpp
art.atmosphere = homAtmos;
runTestReceiverNearGroundSourceElevation(art, 'Homogeneous');
art.atmosphere = inhomAtmos;
runTestReceiverNearGroundSourceElevation(art, 'Inhomogeneous');
function runTestReceiverNearGroundSourceElevation(art, atmosphereSuffix)
sourceReceiverDistance = 2000;
phi = 0;
for sourceRecElevation = -90:45:90
art.source.x = sourceReceiverDistance*sind(sourceRecElevation) + art.receiver.x;
art.source.y = art.receiver.y;
art.source.z = sourceReceiverDistance*cosd(sourceRecElevation) + art.receiver.z;
jsonFolder = fullfile( 'output', fileparts(mfilename('fullpath')) );
cppFilename = ['ReceiverNearGround_' atmosphereSuffix '_SourceElevation' num2str(sourceRecElevation) '_SourceAzimuth' num2str(phi) '.json'];
runMatlabVsCppEigenrayTest(art, fullfile(jsonFolder, cppFilename), fileparts(cppFilename), true);
end
function checkCppEigenrays(cppFile, plotTitle)
%% Load C++ Results
cppRays = itaAtmosRay.load(cppFile);
if exist('AtmosphericRay', 'class')
cppRays = AtmosphericRay.load(cppFile);
else
cppRays = itaAtmosRay.load(cppFile);
end
if nargin == 1
[~, plotTitle, ~] = fileparts(cppFile);
end
%% Plot rays
receiverPos = [0 0 1.8];
......@@ -15,7 +22,11 @@ for idx=1:nRays
receiverDistance = norm(cppRays(idx).r.cart(end,:) - receiverPos);
reflOrder = cppRays(idx).numReflections;
legendEntries{idx} = ['order: ' num2str(reflOrder) ', dist: ' num2str(receiverDistance) 'm'];
cppRays(idx).pr(f, {'linestyle', 'linewidth'}, {'--', 2.25});
if isa(cppRays, 'AtmosphericRay')
cppRays(idx).plot([], ax, 'linestyle', '--', 'linewidth', 2.25);
else
cppRays(idx).pr(f, {'linestyle', 'linewidth'}, {'--', 2.25});
end
end
plot3(cppRays(1).r0.x, cppRays(1).r0.y, cppRays(1).r0.z, ' ok', 'markerfacecolor', 'b')
......
function runMatlabVsCppEigenrayTest(art, cppFile, plotTitleBase, doPlotRays)
function runMatlabVsCppEigenrayTest(art, cppFile, plotTitle, doPlotRays)
if nargin == 3
doPlotRays = false;
......@@ -18,11 +18,11 @@ end
%% Plot rays
if doPlotRays
plotRays(art, matlabRays, cppRays)
title(plotTitleBase)
title(plotTitle)
end
%% Plot distance
plotAbsoluteAndRelativeDistance(matlabRays, cppRays, plotTitleBase)
plotAbsoluteAndRelativeDistance(matlabRays, cppRays, plotTitle)
function plotAbsoluteAndRelativeDistance(matlabRays, cppRays, titleStr)
......