USML
ocean/test/boundary_test.cc
#include <boost/test/unit_test.hpp>
#include <usml/netcdf/netcdf_files.h>
#include <usml/ocean/ocean.h>
#include <fstream>
BOOST_AUTO_TEST_SUITE(boundary_test)
using namespace boost::unit_test;
using namespace usml::netcdf;
using namespace usml::ocean;
BOOST_AUTO_TEST_CASE( flat_boundary_test ) {
cout << "=== boundary_test: flat_boundary_test ===" << endl;
// simple values for points and depth
wposition1 points;
double depth;
wvector1 normal;
// compute profile
boundary_flat model(1000.0);
model.height(points, &depth, &normal);
cout << "depth: " << (wposition::earth_radius - depth) << endl;
cout << "normal rho: " << normal.rho() << endl;
cout << "normal theta: " << normal.theta() << endl;
cout << "normal phi: " << normal.phi() << endl;
// check the answer
BOOST_CHECK_CLOSE(wposition::earth_radius - depth, 1000.0, 1e-6);
BOOST_CHECK_CLOSE(normal.rho(), 1.0, 1e-6);
BOOST_CHECK_CLOSE(normal.theta(), 0.0, 1e-6);
BOOST_CHECK_CLOSE(normal.phi(), 0.0, 1e-6);
}
BOOST_AUTO_TEST_CASE( sloped_boundary_test ) {
cout << "=== boundary_test: sloped_boundary_test ===" << endl;
wposition::compute_earth_radius(0.0);
const double dlat = 1.0;
const double d0 = 1000.0;
const double slope = d0 / (to_radians(dlat)
* (wposition::earth_radius - d0));
const double alpha = atan(slope);
// simple values for points and depth
wposition1 reference;
wposition1 points;
points.latitude(dlat); // define field point 60 nmi north of reference
double depth;
wvector1 normal;
// compute profile
boundary_slope model(reference, d0, alpha);
model.height(points, &depth, &normal);
cout << "slope: " << slope << endl;
cout << "alpha: " << to_degrees(alpha) << endl;
cout << "depth: " << wposition::earth_radius - depth << endl;
cout << "normal rho: " << normal.rho() << endl;
cout << "normal theta: " << normal.theta() << endl;
cout << "normal phi: " << normal.phi() << endl;
// check the answer
BOOST_CHECK(abs(wposition::earth_radius - depth) < 0.1);
BOOST_CHECK_CLOSE(normal.theta(), sin(alpha), 1e-6);
BOOST_CHECK_CLOSE(normal.phi(), 0.0, 1e-6);
BOOST_CHECK_CLOSE(normal.rho(), sqrt(1.0 - sin(alpha) * sin(alpha)), 1e-6);
}
BOOST_AUTO_TEST_CASE( etopo_boundary_test ) {
cout << "=== boundary_test: etopo_boundary_test ===" << endl;
boundary_grid<double, 2> model( new netcdf_bathy(
USML_DATA_DIR "/bathymetry/ETOPO1_Ice_g_gmt4.grd",
36.0, 36.2, 15.85, 16.0, wposition::earth_radius));
// simple values for points and depth
wposition1 points;
points.latitude(36.000447);
points.longitude(15.890411);
double depth;
wvector1 normal;
// compute bathymetry
model.height(points, &depth, &normal);
// check the answer
cout << "Location: lat=" << points.latitude() << " long="
<< points.longitude() << endl << "World Coords: theta="
<< points.theta() << " phi=" << points.phi() << endl
<< "Depth: " << wposition::earth_radius - depth << endl
<< "Normal: theta=" << normal.theta() << " phi="
<< normal.phi() << endl << "Slope: theta=" << to_degrees(
-asin(normal.theta())) << " phi="
<< to_degrees(-asin(normal.phi())) << " deg" << endl;
#ifdef __FAST_MATH__
const double depth_accuracy = 0.005 ;
const double normal_accuracy = 2e-4 ;
#else
const double depth_accuracy = 5e-4 ;
const double normal_accuracy = 2e-4 ;
#endif
BOOST_CHECK_CLOSE(wposition::earth_radius - depth, 3671.1557116601616, depth_accuracy );
BOOST_CHECK( abs(normal.theta()) < normal_accuracy );
BOOST_CHECK( abs(normal.phi() - 0.012764948465248139) < normal_accuracy );
}
BOOST_AUTO_TEST_CASE( ascii_arc_test ) {
cout << "=== boundary_test: ascii_arc_test ===" << endl;
// test interpolation of the raw grid
ascii_arc_bathy* grid = new ascii_arc_bathy(
USML_DATA_DIR "/arcascii/small_crm.asc" );
//USML_TEST_DIR "/ocean/test/ascii_arc_test.asc" ) ;
BOOST_CHECK_EQUAL( grid->axis(0)->size(), 241 ); //rows
BOOST_CHECK_EQUAL( grid->axis(1)->size(), 241 ); //columns
unsigned index[2] ;
index[0]=0; index[1]=0; BOOST_CHECK_CLOSE(wposition::earth_radius - grid->data(index), 684.0, 1e-6);
index[0]=240; index[1]=0; BOOST_CHECK_CLOSE(wposition::earth_radius - grid->data(index), 622.0, 1e-6);
index[0]=0; index[1]=240; BOOST_CHECK_CLOSE(wposition::earth_radius - grid->data(index), 771.0, 1e-6);
index[0]=240; index[1]=240; BOOST_CHECK_CLOSE(wposition::earth_radius - grid->data(index), 747.0, 1e-6);
// test implementation as a boundary model
grid->interp_type(0,GRID_INTERP_PCHIP);
grid->interp_type(1,GRID_INTERP_PCHIP);
std::ofstream before("usml_ascii_arc_interp_before_boundary_grid.csv");
for(int i=0; i<grid->axis(0)->size(); ++i) {
before << to_latitude((*(grid->axis(0)))[i]) << ",";
for(int j=0; j<grid->axis(1)->size(); ++j) {
double location[2];
location[0] = (*(grid->axis(0)))[i];
location[1] = (*(grid->axis(1)))[j];
before << wposition::earth_radius - grid->interpolate(location) << ",";
}
before << endl;
}
for ( int k=0 ; k < grid->axis(1)->size() ; ++k ) {
if(k==0) {before << ",";}
before << to_degrees((*(grid->axis(1)))[k]) << ",";
if(k==grid->axis(1)->size()) {before << endl;}
}
boundary_grid<double,2> bottom(grid) ;
std::ofstream after("usml_ascii_arc_interp_after_boundary_grid.csv");
for(int i=0; i<grid->axis(0)->size(); ++i) {
after << to_latitude((*(grid->axis(0)))[i]) << ",";
for(int j=0; j<grid->axis(1)->size(); ++j) {
double depth;
wposition1 location( to_latitude((*(grid->axis(0)))[i]), to_degrees((*(grid->axis(1)))[j]) );
bottom.height(location, &depth);
after << wposition::earth_radius - depth << ",";
}
after << endl;
}
for ( int k=0 ; k < grid->axis(1)->size() ; ++k ) {
if(k==0) {after << ",";}
after << to_degrees((*(grid->axis(1)))[k]) << ",";
if(k==grid->axis(1)->size()) {after << endl;}
}
wposition1 location( 29.4361, -79.7862 ) ;
double depth ;
bottom.height( location, &depth ) ;
BOOST_CHECK_CLOSE(wposition::earth_radius - depth, 700.0, 0.3);
wposition1 location2( 29.4402, -79.8853 ) ;
bottom.height( location2, &depth ) ;
BOOST_CHECK_CLOSE(wposition::earth_radius - depth, 681.0, 0.3);
}
BOOST_AUTO_TEST_SUITE_END()