#include <boost/test/unit_test.hpp>
#include <usml/waveq3d/waveq3d.h>
#include <usml/netcdf/netcdf_files.h>
#include <iostream>
#include <iomanip>
#include <fstream>
BOOST_AUTO_TEST_SUITE(proploss_test)
using namespace boost::unit_test;
using namespace usml::waveq3d;
static const
double f0 = 2000 ;
{
cout << "=== proploss_test: proploss_basic ===" << endl;
const char* csvname = USML_TEST_DIR "/waveq3d/test/proploss_basic.csv";
const double c0 = 1500.0;
const double src_alt = -15000.0;
const double time_max = 8.0;
wposition::compute_earth_radius(src_lat);
attenuation_model* attn = new attenuation_constant(0.0);
profile_model* profile = new profile_linear(c0, attn);
boundary_model* surface = new boundary_flat();
boundary_model* bottom =
new boundary_flat(
bot_depth);
ocean_model ocean(surface, bottom, profile);
seq_log
freq(10.0, 10.0, 4);
cout << "frequencies: " << freq << endl;
wposition1 pos(src_lat, src_lng, src_alt);
seq_linear de(-10.7, 1.0, 10.0);
seq_linear az(-10.5, 2.0, 10.0);
wposition target(10, 1, src_lat, src_lng, src_alt);
for (unsigned n = 0; n < target.size1(); ++n)
{
target.latitude(n, 0, src_lat + 0.01 * (n + 2.0));
}
proploss loss(freq, pos, de, az, time_step, &target);
wave_queue wave( ocean, freq, pos, de, az, time_step, &target, wave_queue::CLASSIC_RAY ) ;
wave.addProplossListener(&loss);
cout << "propagate wavefronts" << endl;
while (wave.time() < time_max)
{
wave.step();
}
loss.sum_eigenrays();
std::ofstream os(csvname);
os << "target,time,intensity,phase,src de,src az,trg de,trg az,surf,bot" << endl;
os << std::setprecision(18);
cout << "writing spreadsheets to " << csvname << endl;
for (unsigned n = 0; n < target.size1(); ++n)
{
for (eigenray_list::const_iterator iter = raylist->begin();
iter != raylist->end(); ++n, ++iter)
{
const eigenray &ray = *iter;
os << n
<< "," << ray.time
<< "," << ray.intensity(0)
<< "," << ray.phase(0)
<< "," << ray.source_de
<< "," << ray.source_az
<< "," << ray.target_de
<< "," << ray.target_az
<< "," << ray.surface
<< "," << ray.bottom
<< endl;
double range = c0 * ray.time;
double pl = 20.0 *
log10(range - 2.0);
cout << "range=" << range
<< " theory=" << pl
<< " model=" << ray.intensity << endl;
for (unsigned f = 0; f < freq.size(); ++f)
{
BOOST_CHECK(fabs(ray.intensity(f) - pl) < 0.2);
}
}
}
}
{
cout << "=== proploss_test: proploss_lloyds_range ===" << endl;
const char* csvname = USML_TEST_DIR "/waveq3d/test/proploss_lloyds_range.csv";
const char* ncname = USML_TEST_DIR "/waveq3d/test/proploss_lloyds_range.nc";
const char* ncname_wave = USML_TEST_DIR "/waveq3d/test/proploss_lloyds_range_wave.nc";
const double c0 = 1500.0;
const double src_alt = -25.0;
const double trg_alt = -200.0;
const double time_max = 8.0;
wposition::compute_earth_radius(src_lat);
attenuation_model* attn = new attenuation_constant(0.0);
profile_model* profile = new profile_linear(c0, attn);
boundary_model* surface = new boundary_flat();
boundary_model* bottom =
new boundary_flat(
bot_depth);
ocean_model ocean(surface, bottom, profile);
profile->flat_earth(true);
const double wavenum = TWO_PI *
freq(0) /
c0;
wposition1 pos(src_lat, src_lng, src_alt);
seq_rayfan de ;
seq_linear az( -4.0, 1.0, 4.0 );
seq_linear range(200.0, 10.0, 10e3);
for (unsigned n = 0; n < target.size1(); ++n)
{
double degrees = src_lat + range(n) / (1852.0 * 60.0);
target.latitude(n, 0, degrees);
}
wave.addProplossListener(&loss);
cout << "propagate wavefronts" << endl;
cout << "writing wavefronts to " << ncname_wave << endl;
wave.init_netcdf(ncname_wave);
wave.save_netcdf();
while (wave.time() < time_max)
{
wave.step();
wave.save_netcdf();
}
wave.close_netcdf();
cout << "writing eigenrays to " << ncname << endl;
loss.sum_eigenrays();
loss.write_netcdf(ncname);
cout << "writing spreadsheets to " << csvname << endl;
std::ofstream os(csvname);
os << "range,model,theory,m1amp,m1time,t1amp,t1time,m2amp,m2time,t2amp,t2time" << endl;
os << std::setprecision(18);
vector<double> tl_model(range.size());
vector<double> tl_analytic(range.size());
double mean_model = 0.0;
double mean_analytic = 0.0;
const double z1 = trg_alt - src_alt;
const double z2 = trg_alt + src_alt;
for (unsigned n = 0; n < range.size(); ++n)
{
tl_model[n] = -loss.total(n, 0)->intensity(0);
const double R1 =
sqrt(range(n) * range(n) + z1 * z1);
const double R2 =
sqrt(range(n) * range(n) + z2 * z2);
complex<double> p1(0.0, wavenum * R1);
complex<double> p2(0.0, wavenum * R2);
tl_analytic[n] = 10.0 *
log10(norm(p1 + p2));
eigenray_list::const_iterator iter = loss.eigenrays(n, 0)->begin();
if (iter != loss.eigenrays(n, 0)->end())
{
os << range(n)
<< "," << tl_model[n]
<< "," << tl_analytic[n]
<< "," << -(*iter).intensity(0)
<< "," << (*iter).time
<<
"," << 10.0 *
log10(norm(p1))
++iter;
if (iter != loss.eigenrays(n, 0)->end())
{
os << "," << -(*iter).intensity(0)
<< "," << (*iter).time
<<
"," << 10.0 *
log10(norm(p2))
}
}
os << endl;
mean_model += tl_model[n];
mean_analytic += tl_analytic[n];
}
mean_model /= range.size();
mean_analytic /= range.size();
double bias = 0.0;
double dev = 0.0;
double Sxx = 0.0;
double Syy = 0.0;
double Sxy = 0.0;
for (unsigned n = 0; n < range.size(); ++n)
{
const double diff = (tl_model[n] - tl_analytic[n]);
bias += diff ;
dev += (diff * diff);
const double diff_analytic = (tl_analytic[n] - mean_analytic);
Sxx += (diff_analytic * diff_analytic);
const double diff_model = (tl_model[n] - mean_model);
Syy += (diff_model * diff_model);
Sxy += (diff_analytic * diff_model);
}
bias /= range.size() ;
dev =
sqrt( dev / range.size() - bias*bias ) ;
double detcoef = Sxy * Sxy / (Sxx * Syy) * 100.0;
cout << std::setprecision(4);
cout << "bias = " << bias << " dB"
<< " dev = " << dev << " dB"
<< " detcoef = " << detcoef << "%" << endl ;
BOOST_CHECK(
abs(bias) <= 0.5 );
BOOST_CHECK( dev <= 4.0 );
BOOST_CHECK( detcoef >= 80.0 );
}
{
cout << "=== proploss_test: proploss_lloyds_range_freq ===" << endl;
const char* csvname = USML_TEST_DIR "/waveq3d/test/proploss_lloyds_range_freq.csv";
const char* ncname = USML_TEST_DIR "/waveq3d/test/proploss_lloyds_range_freq.nc";
const char* ncname_wave = USML_TEST_DIR "/waveq3d/test/proploss_lloyds_range_freq_wave.nc";
const double c0 = 1500.0;
const double src_lat = 45.0;
const double src_lng = -45.0;
const double src_alt = -25.0;
const double trg_alt = -200.0;
const double time_max = 8.0;
wposition::compute_earth_radius(src_lat);
attenuation_model* attn = new attenuation_constant(0.0);
profile_model* profile = new profile_linear(c0, attn);
boundary_model* surface = new boundary_flat();
boundary_model* bottom =
new boundary_flat(
bot_depth);
ocean_model ocean(surface, bottom, profile);
profile->flat_earth(true);
seq_log
freq(10.0, 10.0, 4);
cout <<
"frequencies: " <<
freq << endl;
wposition1 pos(src_lat, src_lng, src_alt);
seq_rayfan de ;
seq_linear az( -4.0, 1.0, 4.0 );
seq_linear range(200.0, 10.0, 10e3);
for (unsigned n = 0; n < target.size1(); ++n)
{
double degrees = src_lat + range(n) / (1852.0 * 60.0);
target.latitude(n, 0, degrees);
}
wave.addProplossListener(&loss);
cout << "propagate wavefronts" << endl;
cout << "writing wavefronts to " << ncname_wave << endl;
wave.init_netcdf(ncname_wave);
wave.save_netcdf();
while (wave.time() < time_max)
{
wave.step();
wave.save_netcdf();
}
wave.close_netcdf();
cout << "writing eigenrays to " << ncname << endl;
loss.sum_eigenrays();
loss.write_netcdf(ncname);
cout << "writing spreadsheets to " << csvname << endl;
std::ofstream os(csvname);
os << "range,model,theory,m1amp,m1time,t1amp,t1time,m2amp,m2time,t2amp,t2time" << endl;
os << std::setprecision(18);
vector<double> tl_model(range.size());
vector<double> tl_analytic(range.size());
double mean_model = 0.0;
double mean_analytic = 0.0;
const double z1 = trg_alt - src_alt;
const double z2 = trg_alt + src_alt;
for (
unsigned f=0; f <
freq.size(); ++f)
{
const double wavenum = TWO_PI *
freq(f) /
c0;
for (unsigned n = 0; n < range.size(); ++n)
{
tl_model[n] = -loss.total(n, 0)->intensity(f);
const double R1 =
sqrt(range(n) * range(n) + z1 * z1);
const double R2 =
sqrt(range(n) * range(n) + z2 * z2);
complex<double> p1(0.0, wavenum * R1);
complex<double> p2(0.0, wavenum * R2);
tl_analytic[n] = 10.0 *
log10(norm(p1 + p2));
eigenray_list::const_iterator iter = loss.eigenrays(n, 0)->begin();
if (iter != loss.eigenrays(n, 0)->end())
{
os << range(n)
<< "," << tl_model[n]
<< "," << tl_analytic[n]
<< "," << -(*iter).intensity(f)
<< "," << (*iter).time
<<
"," << 10.0 *
log10(norm(p1))
++iter;
if (iter != loss.eigenrays(n, 0)->end())
{
os << "," << -(*iter).intensity(f)
<< "," << (*iter).time
<<
"," << 10.0 *
log10(norm(p2))
}
}
os << endl;
mean_model += tl_model[n];
mean_analytic += tl_analytic[n];
}
mean_model /= range.size();
mean_analytic /= range.size();
double bias = 0.0;
double dev = 0.0;
double Sxx = 0.0;
double Syy = 0.0;
double Sxy = 0.0;
for (unsigned n = 0; n < range.size(); ++n)
{
const double diff = (tl_model[n] - tl_analytic[n]);
bias += diff ;
dev += (diff * diff);
const double diff_analytic = (tl_analytic[n] - mean_analytic);
Sxx += (diff_analytic * diff_analytic);
const double diff_model = (tl_model[n] - mean_model);
Syy += (diff_model * diff_model);
Sxy += (diff_analytic * diff_model);
}
bias /= range.size() ;
dev =
sqrt( dev / range.size() - bias*bias ) ;
double detcoef = Sxy * Sxy / (Sxx * Syy) * 100.0;
cout << std::setprecision(4);
cout <<
" freq =" <<
freq(f) <<
"Hz: bias = " << bias <<
" dB"
<< " dev = " << dev << " dB"
<< " detcoef = " << detcoef << "%" << endl ;
BOOST_CHECK( detcoef >= 80.0 );
int iFreq = (int)
freq(f);
switch (iFreq) {
case 10:
BOOST_CHECK(
abs(bias) <= 10.0 );
BOOST_CHECK( dev <= 5.0 );
break;
case 100:
BOOST_CHECK(
abs(bias) <= 1.0 );
BOOST_CHECK( dev <= 4.0 );
break;
case 10000:
BOOST_CHECK( dev <= 5.0 );
break;
default:
BOOST_CHECK(
abs(bias) <= 0.5 );
BOOST_CHECK( dev <= 4.0 );
}
}
}
{
cout << "=== proploss_test: proploss_lloyds_depth ===" << endl;
const char* csvname = USML_TEST_DIR "/waveq3d/test/proploss_lloyds_depth.csv";
const char* ncname = USML_TEST_DIR "/waveq3d/test/proploss_lloyds_depth.nc";
const char* ncname_wave = USML_TEST_DIR "/waveq3d/test/proploss_lloyds_depth_wave.nc";
const double c0 = 1500.0;
const double src_lat = 45.0;
const double src_lng = -45.0;
const double src_alt = -25.0;
const double range = 10e3;
const double time_max = 8.0;
wposition::compute_earth_radius(src_lat);
attenuation_model* attn = new attenuation_constant(0.0);
profile_model* profile = new profile_linear(c0, attn);
boundary_model* surface = new boundary_flat();
boundary_model* bottom =
new boundary_flat(
bot_depth);
ocean_model ocean(surface, bottom, profile);
profile->flat_earth(true);
const double wavenum = TWO_PI *
freq(0) /
c0;
wposition1 pos(src_lat, src_lng, src_alt);
seq_rayfan de ;
seq_linear az( -4.0, 1.0, 4.0 );
double degrees = src_lat +
to_degrees(range / (wposition::earth_radius+src_alt));
seq_linear depth(-0.1, -0.5, -40.1);
wposition target(depth.size(), 1, degrees,
src_lng, 0.0);
for (unsigned n = 0; n < target.size1(); ++n)
{
target.altitude(n, 0, depth(n));
}
wave.addProplossListener(&loss);
cout << "propagate wavefronts" << endl;
cout << "writing wavefronts to " << ncname_wave << endl;
wave.init_netcdf(ncname_wave);
wave.save_netcdf();
while (wave.time() < time_max)
{
wave.step();
wave.save_netcdf();
}
wave.close_netcdf();
cout << "writing eigenrays to " << ncname << endl;
loss.sum_eigenrays();
loss.write_netcdf(ncname);
cout << "writing spreadsheets to " << csvname << endl;
std::ofstream os(csvname);
os << "depth,model,theory,m1amp,t1amp,m1time,t1time,time1-diff,m2amp,t2amp,m2time,t2time,time2-diff,phase"
<< endl;
os << std::setprecision(18);
vector<double> tl_model(depth.size());
vector<double> tl_analytic(depth.size());
double mean_model = 0.0;
double mean_analytic = 0.0;
for (unsigned n = 0; n < depth.size(); ++n)
{
const double z1 = depth(n) - src_alt;
const double z2 = depth(n) + src_alt;
tl_model[n] = -loss.total(n, 0)->intensity(0);
const double R1 =
sqrt(range * range + z1 * z1);
const double R2 =
sqrt(range * range + z2 * z2);
complex<double> p1(0.0, wavenum * R1);
complex<double> p2(0.0, wavenum * R2);
tl_analytic[n] = 10.0 *
log10(norm(p1 + p2));
eigenray_list::const_iterator iter = loss.eigenrays(n, 0)->begin();
if (iter != loss.eigenrays(n, 0)->end())
{
os << depth(n)
<< "," << tl_model[n]
<< "," << tl_analytic[n]
<< "," << -(*iter).intensity(0)
<<
"," << 10.0 *
log10(norm(p1))
<< "," << (*iter).time
<< "," << R1 / c0
<<
"," << (*iter).time-R1 /
c0;
++iter;
if (iter != loss.eigenrays(n, 0)->end())
{
os << "," << -(*iter).intensity(0)
<<
"," << 10.0 *
log10(norm(p2))
<< "," << (*iter).time
<< "," << R2 / c0
<<
"," << (*iter).time-R2 /
c0;
}
os << "," << (*iter).phase(0);
}
os << endl;
mean_model += tl_model[n];
mean_analytic += tl_analytic[n];
}
mean_model /= depth.size();
mean_analytic /= depth.size();
double bias = 0.0;
double dev = 0.0;
double Sxx = 0.0;
double Syy = 0.0;
double Sxy = 0.0;
for (unsigned n = 0; n < depth.size(); ++n)
{
const double diff = (tl_model[n] - tl_analytic[n]);
bias += diff ;
dev += (diff * diff);
const double diff_analytic = (tl_analytic[n] - mean_analytic);
Sxx += (diff_analytic * diff_analytic);
const double diff_model = (tl_model[n] - mean_model);
Syy += (diff_model * diff_model);
Sxy += (diff_analytic * diff_model);
}
bias /= depth.size() ;
dev =
sqrt( dev / depth.size() - bias*bias ) ;
double detcoef = Sxy * Sxy / (Sxx * Syy) * 100.0;
cout << std::setprecision(4);
cout << "bias = " << bias << " dB"
<< " dev = " << dev << " dB"
<< " detcoef = " << detcoef << "%" << endl ;
BOOST_CHECK(
abs(bias) <= 0.7 );
BOOST_CHECK( dev <= 4.0 );
BOOST_CHECK( detcoef >= 80.0 );
}
int month = 1 ;
const double lat1 = 24.0 ;
const double lat2 = 26.0 ;
const double lng1 = 56.0 ;
const double lng2 = 58.0 ;
wposition::compute_earth_radius( (lat1+lat2)/2.0 ) ;
const double src_lat = 26.0 ;
const double src_lng = 56.75 ;
const double src_alt = -27.0 ;
const double tar_range = 1000.0 ;
const double tar_alt = -50.0 ;
const double time_max = 4.0 ;
const double dt = 0.05 ;
seq_log
freq( 6500.0, 0.0, 1 ) ;
wposition1 pos( src_lat, src_lng, src_alt ) ;
seq_rayfan de( -34.0, 36.0, 21 ) ;
seq_linear az( 0.0, 15.0, 360.0 ) ;
cout << "loading temperature & salinity data from World Ocean Atlas" << endl ;
USML_DATA_DIR "/woa09/temperature_seasonal_1deg.nc",
USML_DATA_DIR "/woa09/temperature_monthly_1deg.nc",
month, lat1, lat2, lng1, lng2 ) ;
USML_DATA_DIR "/woa09/salinity_seasonal_1deg.nc",
USML_DATA_DIR "/woa09/salinity_monthly_1deg.nc",
month, lat1, lat2, lng1, lng2 ) ;
profile_model* profile = new profile_grid<double,3>(
data_grid_mackenzie::construct(temperature, salinity) ) ;
cout << "loading bathymetry from ETOPO1 database" << endl ;
USML_DATA_DIR "/bathymetry/ETOPO1_Ice_g_gmt4.grd", lat1, lat2, lng1, lng2 ) );
bottom->reflect_loss( new reflect_loss_rayleigh( reflect_loss_rayleigh::SILT ) ) ;
boundary_model* surface = new boundary_flat() ;
ocean_model ocean( surface, bottom, profile ) ;
wposition target( 1, 1, 0.0, 0.0, tar_alt ) ;
wposition1 atarget( pos, tar_range, 0.0 ) ;
target.latitude( 0, 0, atarget.latitude() ) ;
target.longitude( 0, 0, atarget.longitude() ) ;
proploss* pLoss =
new proploss(
freq, pos, de, az, dt, &target) ;
wave_queue* pWave =
new wave_queue( ocean,
freq, pos, de, az, dt, &target ) ;
pWave->addProplossListener(pLoss);
while(pWave->time() < time_max) {
pWave->step();
}
pLoss->sum_eigenrays();
cout << "\n\n=====Eigenrays --> Rayleigh::SILT=====" << endl ;
cout << "time (s)\tbounces(s,b,c)\tlaunch angle\tarrival angle\t TL\t\t phase" << endl ;
cout << std::setprecision(5) ;
for (eigenray_list::const_iterator iter = pLoss->eigenrays(0, 0)->begin(); iter != pLoss->eigenrays(0, 0)->end(); ++iter)
{
cout << (*iter).time
<< "\t\t (" << (*iter).surface << ", " << (*iter).bottom << ", " << (*iter).caustic << ")"
<< "\t " << (*iter).source_de
<< "\t\t " << (*iter).target_de
<< "\t\t " << (*iter).intensity(0)
<< "\t " << (*iter).phase(0)
<< endl;
}
bottom->reflect_loss( new reflect_loss_rayleigh( reflect_loss_rayleigh::SAND ) ) ;
delete pLoss ;
delete pWave ;
pLoss =
new proploss(
freq, pos, de, az, dt, &target) ;
pWave =
new wave_queue( ocean,
freq, pos, de, az, dt, &target ) ;
pWave->addProplossListener(pLoss) ;
while(pWave->time() < time_max) {
pWave->step();
}
pLoss->sum_eigenrays();
cout << "=====Eigenrays --> Rayleigh::SAND=====" << endl ;
cout << "time (s)\tbounces(s,b,c)\tlaunch angle\tarrival angle\t TL\t\t phase" << endl ;
cout << std::setprecision(5) ;
for (eigenray_list::const_iterator iter = pLoss->eigenrays(0, 0)->begin(); iter != pLoss->eigenrays(0, 0)->end(); ++iter)
{
cout << (*iter).time
<< "\t\t (" << (*iter).surface << ", " << (*iter).bottom << ", " << (*iter).caustic << ")"
<< "\t " << (*iter).source_de
<< "\t\t " << (*iter).target_de
<< "\t\t " << (*iter).intensity(0)
<< "\t " << (*iter).phase(0)
<< endl;
}
bottom->reflect_loss( new reflect_loss_rayleigh( reflect_loss_rayleigh::CLAY ) ) ;
delete pLoss ;
delete pWave ;
pLoss =
new proploss(
freq, pos, de, az, dt, &target) ;
pWave =
new wave_queue( ocean,
freq, pos, de, az, dt, &target ) ;
pWave->addProplossListener(pLoss) ;
while(pWave->time() < time_max) {
pWave->step();
}
pLoss->sum_eigenrays();
cout << "=====Eigenrays --> Rayleigh::CLAY=====" << endl ;
cout << "time (s)\tbounces(s,b,c)\tlaunch angle\tarrival angle\t TL\t\t phase" << endl ;
cout << std::setprecision(5) ;
for (eigenray_list::const_iterator iter = pLoss->eigenrays(0, 0)->begin(); iter != pLoss->eigenrays(0, 0)->end(); ++iter)
{
cout << (*iter).time
<< "\t\t (" << (*iter).surface << ", " << (*iter).bottom << ", " << (*iter).caustic << ")"
<< "\t " << (*iter).source_de
<< "\t\t " << (*iter).target_de
<< "\t\t " << (*iter).intensity(0)
<< "\t " << (*iter).phase(0)
<< endl;
}
bottom->reflect_loss( new reflect_loss_rayleigh( reflect_loss_rayleigh::LIMESTONE ) ) ;
delete pLoss ;
delete pWave ;
pLoss =
new proploss(
freq, pos, de, az, dt, &target) ;
pWave =
new wave_queue( ocean,
freq, pos, de, az, dt, &target ) ;
pWave->addProplossListener(pLoss) ;
while(pWave->time() < time_max) {
pWave->step();
}
pLoss->sum_eigenrays();
cout << "=====Eigenrays --> Rayleigh::LIMESTONE=====" << endl ;
cout << "time (s)\tbounces(s,b,c)\tlaunch angle\tarrival angle\t TL\t\t phase" << endl ;
cout << std::setprecision(5) ;
for (eigenray_list::const_iterator iter = pLoss->eigenrays(0, 0)->begin(); iter != pLoss->eigenrays(0, 0)->end(); ++iter)
{
cout << (*iter).time
<< "\t\t (" << (*iter).surface << ", " << (*iter).bottom << ", " << (*iter).caustic << ")"
<< "\t " << (*iter).source_de
<< "\t\t " << (*iter).target_de
<< "\t\t " << (*iter).intensity(0)
<< "\t " << (*iter).phase(0)
<< endl;
}
bottom->reflect_loss( new reflect_loss_rayleigh( reflect_loss_rayleigh::BASALT ) ) ;
delete pLoss ;
delete pWave ;
pLoss =
new proploss(
freq, pos, de, az, dt, &target) ;
pWave =
new wave_queue( ocean,
freq, pos, de, az, dt, &target ) ;
pWave->addProplossListener(pLoss) ;
while(pWave->time() < time_max) {
pWave->step();
}
pLoss->sum_eigenrays();
cout << "=====Eigenrays --> Rayleigh::BASALT=====" << endl ;
cout << "time (s)\tbounces(s,b,c)\tlaunch angle\tarrival angle\t TL\t\t phase" << endl ;
cout << std::setprecision(5) ;
for (eigenray_list::const_iterator iter = pLoss->eigenrays(0, 0)->begin(); iter != pLoss->eigenrays(0, 0)->end(); ++iter)
{
cout << (*iter).time
<< "\t\t (" << (*iter).surface << ", " << (*iter).bottom << ", " << (*iter).caustic << ")"
<< "\t " << (*iter).source_de
<< "\t\t " << (*iter).target_de
<< "\t\t " << (*iter).intensity(0)
<< "\t " << (*iter).phase(0)
<< endl;
}
bottom->reflect_loss( new reflect_loss_rayleigh( reflect_loss_rayleigh::GRAVEL ) ) ;
delete pLoss ;
delete pWave ;
pLoss =
new proploss(
freq, pos, de, az, dt, &target) ;
pWave =
new wave_queue( ocean,
freq, pos, de, az, dt, &target ) ;
pWave->addProplossListener(pLoss) ;
while(pWave->time() < time_max) {
pWave->step();
}
pLoss->sum_eigenrays();
cout << "=====Eigenrays --> Rayleigh::GRAVEL=====" << endl ;
cout << "time (s)\tbounces(s,b,c)\tlaunch angle\tarrival angle\t TL\t\t phase" << endl ;
cout << std::setprecision(5) ;
for (eigenray_list::const_iterator iter = pLoss->eigenrays(0, 0)->begin(); iter != pLoss->eigenrays(0, 0)->end(); ++iter)
{
cout << (*iter).time
<< "\t\t (" << (*iter).surface << ", " << (*iter).bottom << ", " << (*iter).caustic << ")"
<< "\t " << (*iter).source_de
<< "\t\t " << (*iter).target_de
<< "\t\t " << (*iter).intensity(0)
<< "\t " << (*iter).phase(0)
<< endl;
}
bottom->reflect_loss( new reflect_loss_rayleigh( reflect_loss_rayleigh::MORAINE ) ) ;
delete pLoss ;
delete pWave ;
pLoss =
new proploss(
freq, pos, de, az, dt, &target) ;
pWave =
new wave_queue( ocean,
freq, pos, de, az, dt, &target ) ;
pWave->addProplossListener(pLoss) ;
while(pWave->time() < time_max) {
pWave->step();
}
pLoss->sum_eigenrays();
cout << "=====Eigenrays --> Rayleigh::MORAINE=====" << endl ;
cout << "time (s)\tbounces(s,b,c)\tlaunch angle\tarrival angle\t TL\t\t phase" << endl ;
cout << std::setprecision(5) ;
for (eigenray_list::const_iterator iter = pLoss->eigenrays(0, 0)->begin(); iter != pLoss->eigenrays(0, 0)->end(); ++iter)
{
cout << (*iter).time
<< "\t\t (" << (*iter).surface << ", " << (*iter).bottom << ", " << (*iter).caustic << ")"
<< "\t " << (*iter).source_de
<< "\t\t " << (*iter).target_de
<< "\t\t " << (*iter).intensity(0)
<< "\t " << (*iter).phase(0)
<< endl;
}
bottom->reflect_loss( new reflect_loss_rayleigh( reflect_loss_rayleigh::CHALK ) ) ;
delete pLoss ;
delete pWave ;
pLoss =
new proploss(
freq, pos, de, az, dt, &target) ;
pWave =
new wave_queue( ocean,
freq, pos, de, az, dt, &target ) ;
pWave->addProplossListener(pLoss) ;
while(pWave->time() < time_max) {
pWave->step();
}
pLoss->sum_eigenrays();
cout << "=====Eigenrays --> Rayleigh::CHALK=====" << endl ;
cout << "time (s)\tbounces(s,b,c)\tlaunch angle\tarrival angle\t TL\t\t phase" << endl ;
cout << std::setprecision(5) ;
for (eigenray_list::const_iterator iter = pLoss->eigenrays(0, 0)->begin(); iter != pLoss->eigenrays(0, 0)->end(); ++iter)
{
cout << (*iter).time
<< "\t\t (" << (*iter).surface << ", " << (*iter).bottom << ", " << (*iter).caustic << ")"
<< "\t " << (*iter).source_de
<< "\t\t " << (*iter).target_de
<< "\t\t " << (*iter).intensity(0)
<< "\t " << (*iter).phase(0)
<< endl;
}
bottom->reflect_loss( new reflect_loss_rayleigh( reflect_loss_rayleigh::MUD ) ) ;
delete pLoss ;
delete pWave ;
pLoss =
new proploss(
freq, pos, de, az, dt, &target) ;
pWave =
new wave_queue( ocean,
freq, pos, de, az, dt, &target ) ;
pWave->addProplossListener(pLoss) ;
while(pWave->time() < time_max) {
pWave->step();
}
pLoss->sum_eigenrays();
cout << "=====Eigenrays --> Rayleigh::MUD=====" << endl ;
cout << "time (s)\tbounces(s,b,c)\tlaunch angle\tarrival angle\t TL\t\t phase" << endl ;
cout << std::setprecision(5) ;
for (eigenray_list::const_iterator iter = pLoss->eigenrays(0, 0)->begin(); iter != pLoss->eigenrays(0, 0)->end(); ++iter)
{
cout << (*iter).time
<< "\t\t (" << (*iter).surface << ", " << (*iter).bottom << ", " << (*iter).caustic << ")"
<< "\t " << (*iter).source_de
<< "\t\t " << (*iter).target_de
<< "\t\t " << (*iter).intensity(0)
<< "\t " << (*iter).phase(0)
<< endl;
}
}
BOOST_AUTO_TEST_SUITE_END()