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

    size_t 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_TEST_DIR
                "/ocean/test/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_TEST_DIR
                "/ocean/test/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_CASE( scattering_strength_test ) {
    cout << "=== boundary_test: scattering_strength_test ===" << endl;

    const wposition1 pos ;
    const seq_linear freq( 100.0, 0.0, 1 ) ;
    vector<double> amplitude( freq.size() ) ;
    vector<double> phase( freq.size() ) ;
    scattering_model* s = new scattering_lambert() ;
    const double de_scattered = M_PI / 4.0 ;

    const char* csv_name = USML_TEST_DIR "/ocean/test/scattering_lambert_test.csv" ;
    std::ofstream os(csv_name) ;
    cout << "writing tables to " << csv_name << endl ;
    os << "de_incident,de_scattered,amp" << endl ;
    for(int i=0; i<90; ++i) {
        double de_incident = i * M_PI / 180.0 ;
        s->scattering( pos, freq, de_incident, de_scattered, 0.0, 0.0, &amplitude ) ;
        os << de_incident << ","
           << de_scattered << ","
           << amplitude(0) << endl ;
    }
    delete s ;
}

BOOST_AUTO_TEST_CASE( ocean_volume_test ) {
    cout << "=== boundary_test: ocean_volume_test ===" << endl;

    ocean_model ocean1(
                new boundary_flat(),
                new boundary_flat(2000.0),
                new profile_linear() ) ;
    ocean1.add_volume( new volume_flat( 1000.0, 10.0, -30.0 ) ) ;

    wposition1 location(0.0,0.0) ;
    double depth, thickness ;
    ocean1.volume(0).depth( location, &depth, &thickness ) ;
    BOOST_CHECK_CLOSE(depth, wposition::earth_radius-1000.0, 1e-6);
    BOOST_CHECK_CLOSE(thickness, 10.0, 1e-6);

    seq_linear frequencies(10.0,10.0,3) ;
    boost::numeric::ublas::vector<double> amplitude( frequencies.size() ) ;
    ocean1.volume(0).scattering( location, frequencies,
                0.0, 0.0, 0.0, 0.0, &amplitude ) ;
    BOOST_CHECK_CLOSE(amplitude(2), 1E-3, 1e-6);
}


BOOST_AUTO_TEST_SUITE_END()

Generated on 4 May 2015 for USML by  doxygen 1.6.1