USML
waveq3d/test/proploss_test.cc
#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 time_step = 0.100 ;
static const double f0 = 2000 ;
static const double bot_depth = 1e5 ;
BOOST_AUTO_TEST_CASE(proploss_basic)
{
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_lat = 45.0;
const double src_lng = -45.0;
const double src_alt = -15000.0;
const double time_max = 8.0;
// initialize propagation model
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);
// build a series of targets at different ranges
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);
// propagate rays & record to log file
cout << "propagate wavefronts" << endl;
while (wave.time() < time_max)
{
wave.step();
}
loss.sum_eigenrays();
// save results to spreadsheet and compare to analytic results
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)
{
const eigenray_list *raylist = loss.eigenrays(n, 0);
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);
}
}
}
}
BOOST_AUTO_TEST_CASE(proploss_lloyds_range)
{
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_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;
// initialize propagation model
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(f0, 1.0, 1); // 2000 Hz
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 );
// build a series of targets at different ranges
seq_linear range(200.0, 10.0, 10e3); // range in meters
wposition target(range.size(), 1, src_lat, src_lng, trg_alt);
for (unsigned n = 0; n < target.size1(); ++n)
{
double degrees = src_lat + range(n) / (1852.0 * 60.0); // range in latitude
target.latitude(n, 0, degrees);
}
proploss loss(freq, pos, de, az, time_step, &target);
wave_queue wave( ocean, freq, pos, de, az, time_step, &target) ;
wave.addProplossListener(&loss);
// propagate rays & record to log file
cout << "propagate wavefronts" << endl;
cout << "writing wavefronts to " << ncname_wave << endl;
wave.init_netcdf(ncname_wave); // open a log file for wavefront data
wave.save_netcdf(); // write ray data to log file
while (wave.time() < time_max)
{
wave.step();
wave.save_netcdf(); // write ray data to log file
}
wave.close_netcdf(); // close log file for wavefront data
cout << "writing eigenrays to " << ncname << endl;
loss.sum_eigenrays();
loss.write_netcdf(ncname);
// save results to spreadsheet and compare to analytic results
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);
// compute analytic solution
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);
p1 = exp(p1) / R1;
p2 = -exp(p2) / R2;
tl_analytic[n] = 10.0 * log10(norm(p1 + p2));
// print to log file
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))
<< "," << R1 / c0;
++iter;
if (iter != loss.eigenrays(n, 0)->end())
{
os << "," << -(*iter).intensity(0)
<< "," << (*iter).time
<< "," << 10.0 * log10(norm(p2))
<< "," << R2 / c0;
}
}
os << endl;
mean_model += tl_model[n];
mean_analytic += tl_analytic[n];
}
mean_model /= range.size();
mean_analytic /= range.size();
// compute statistics of difference between curves
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 );
}
BOOST_AUTO_TEST_CASE(proploss_lloyds_range_freq)
{
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;
// initialize propagation model
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 );
// build a series of targets at different ranges
seq_linear range(200.0, 10.0, 10e3); // range in meters
wposition target(range.size(), 1, src_lat, src_lng, trg_alt);
for (unsigned n = 0; n < target.size1(); ++n)
{
double degrees = src_lat + range(n) / (1852.0 * 60.0); // range in latitude
target.latitude(n, 0, degrees);
}
proploss loss(freq, pos, de, az, time_step, &target);
wave_queue wave( ocean, freq, pos, de, az, time_step, &target) ;
wave.addProplossListener(&loss);
// propagate rays & record to log file
cout << "propagate wavefronts" << endl;
cout << "writing wavefronts to " << ncname_wave << endl;
wave.init_netcdf(ncname_wave); // open a log file for wavefront data
wave.save_netcdf(); // write ray data to log file
while (wave.time() < time_max)
{
wave.step();
wave.save_netcdf(); // write ray data to log file
}
wave.close_netcdf(); // close log file for wavefront data
cout << "writing eigenrays to " << ncname << endl;
loss.sum_eigenrays();
loss.write_netcdf(ncname);
// save results to spreadsheet and compare to analytic results
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);
// compute analytic solution
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);
p1 = exp(p1) / R1;
p2 = -exp(p2) / R2;
tl_analytic[n] = 10.0 * log10(norm(p1 + p2));
// print to log file
eigenray_list::const_iterator iter = loss.eigenrays(n, 0)->begin();
if (iter != loss.eigenrays(n, 0)->end())
{
os << range(n) // range
<< "," << tl_model[n] // model
<< "," << tl_analytic[n] // theory
<< "," << -(*iter).intensity(f) // m1amp
<< "," << (*iter).time // m1time
<< "," << 10.0 * log10(norm(p1)) // t1amp
<< "," << R1 / c0; // t1time
++iter;
if (iter != loss.eigenrays(n, 0)->end())
{
os << "," << -(*iter).intensity(f) // m2amp
<< "," << (*iter).time // m2time
<< "," << 10.0 * log10(norm(p2)) // t2amp
<< "," << R2 / c0; // t2time
}
}
os << endl;
mean_model += tl_model[n];
mean_analytic += tl_analytic[n];
}
mean_model /= range.size();
mean_analytic /= range.size();
// compute statistics of difference between curves for each freq
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 );
// Extreme low and high freq can have a greater dev
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 );
}
} // end Freq
}
BOOST_AUTO_TEST_CASE(proploss_lloyds_depth)
{
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;
// initialize propagation model
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(f0, 1.0, 1);
const double wavenum = TWO_PI * freq(0) / c0;
wposition1 pos(src_lat, src_lng, src_alt);
// seq_rayfan de( -17.0, 17.0, 363 ) ;
// seq_linear de( -5.0, 5.0, 720 );
seq_rayfan de ;
seq_linear az( -4.0, 1.0, 4.0 );
// build a series of targets at different depths
double degrees = src_lat + to_degrees(range / (wposition::earth_radius+src_alt)); // range in latitude
seq_linear depth(-0.1, -0.5, -40.1); // depth in meters
// seq_linear depth(-0.1, 1.0, 1); // depth in meters
wposition target(depth.size(), 1, degrees, src_lng, 0.0);
for (unsigned n = 0; n < target.size1(); ++n)
{
target.altitude(n, 0, depth(n));
}
// wposition target(1,1,degrees,src_lng,-25.0);
proploss loss(freq, pos, de, az, time_step, &target);
wave_queue wave( ocean, freq, pos, de, az, time_step, &target) ;
wave.addProplossListener(&loss);
// propagate rays & record to log file
cout << "propagate wavefronts" << endl;
cout << "writing wavefronts to " << ncname_wave << endl;
wave.init_netcdf(ncname_wave); // open a log file for wavefront data
wave.save_netcdf(); // write ray data to log file
while (wave.time() < time_max)
{
wave.step();
wave.save_netcdf(); // write ray data to log file
}
wave.close_netcdf(); // close log file for wavefront data
cout << "writing eigenrays to " << ncname << endl;
loss.sum_eigenrays();
loss.write_netcdf(ncname);
// save results to spreadsheet and compare to analytic results
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);
// compute analytic solution
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);
p1 = exp(p1) / R1;
p2 = -exp(p2) / R2;
tl_analytic[n] = 10.0 * log10(norm(p1 + p2));
// print to log file
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();
// compute statistics of difference between curves
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 );
}
BOOST_AUTO_TEST_CASE( bottom_type_effects ) {
int month = 1 ; // January
const double lat1 = 24.0 ; // gulf of oman
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 ) ;
// build sound velocity profile from World Ocean Atlas data
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) ) ;
// load bathymetry from ETOPO1 database
cout << "loading bathymetry from ETOPO1 database" << endl ;
boundary_model* bottom = new boundary_grid<double,2>( new usml::netcdf::netcdf_bathy(
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() ;
// combine sound speed and bathymetry into ocean model
ocean_model ocean( surface, bottom, profile ) ;
// create a target due north tar_range away
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()