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>
#include <time.h>

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 (size_t 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, 0, wave_queue::CLASSIC_RAY ) ;
    wave.addEigenrayListener(&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 (size_t 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 (size_t f = 0; f < freq.size(); ++f)
            {
                BOOST_CHECK(fabs(ray.intensity(f) - pl) < 0.2);
            }
        }
    }
}

BOOST_AUTO_TEST_CASE( proploss_freq )
{
    cout << "=== proploss_test: proploss_freq ===" << endl ;
    const char* ncwave = USML_TEST_DIR "/waveq3d/test/proploss_freq_wave.nc" ;
    const char* ncproploss = USML_TEST_DIR "/waveq3d/test/proploss_freq_proploss.nc" ;
    const char* csvname = USML_TEST_DIR "/waveq3d/test/proploss_freq.csv";
    const double c0 = 1500.0 ;
    const double src_lat = 45.0 ;
    const double src_lng = -45.0 ;
    const double src_alt = -1000.0 ;
    const double trg_lat = 45.02 ;
    const double time_max = 3.5 ;

    // 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(3000.0) ;
    ocean_model ocean( surface, bottom, profile ) ;

    seq_linear freq( 10.0, 10.0, 1000 ) ;
    wposition1 pos( src_lat, src_lng, src_alt ) ;
    seq_linear de( -90.0, 1.0, 90.0 ) ;
    seq_linear az( -4.0, 1.0, 4.0 ) ;

    // build a single target

    wposition target( 1, 1, trg_lat, src_lng, src_alt );

    proploss loss(freq, pos, de, az, time_step, &target);
    wave_queue wave( ocean, freq, pos, de, az, time_step, &target) ;
    wave.addEigenrayListener(&loss) ;

    // propagate rays & record to log file

    cout << "propagate wavefront for " << time_max << endl ;
    cout << "writing wavefronts to " << ncwave << endl ;

    wave.init_netcdf( ncwave ) ;
    wave.save_netcdf() ;
    while (wave.time() < time_max)
    {
        wave.step() ;
        wave.save_netcdf() ;
    }
    wave.close_netcdf() ;

    loss.sum_eigenrays() ;
    cout << "writing proploss to " << ncproploss << endl ;
    loss.write_netcdf(ncproploss,"proploss_freq test") ;

    // save results to spreadsheet and compare to analytic results

    cout << "writing data to " << csvname << endl ;
    std::ofstream os(csvname) ;
    os << "frequency,theory_direct,model_direct,theory_surface,model_surface,theory_bottom,theory_surface" << endl ;
    os << std::setprecision(18) ;

    const eigenray_list *raylist = loss.eigenrays(0, 0);
    eigenray direct_ray, surface_ray, bottom_ray ;
    for (eigenray_list::const_iterator iter = raylist->begin();
            iter != raylist->end(); ++iter)
    {
        if( iter->bottom == 1 && iter->surface == 0 ) bottom_ray = *iter ;
        else if( iter->bottom == 0 && iter->surface == 1 ) surface_ray = *iter ;
        else direct_ray = *iter ;
    }
    for(size_t j=0; j<freq.size(); ++j) {
        os << freq(j)
           << "," << 66.95
           << "," << direct_ray.intensity(j)
           << "," << 69.52
           << "," << surface_ray.intensity(j)
           << "," << 73.21
           << "," << bottom_ray.intensity(j)
           << endl ;
    }
}

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 (size_t 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.addEigenrayListener(&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 (size_t 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 (size_t 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 (size_t 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.addEigenrayListener(&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 << "freq,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 (size_t f=0; f < freq.size(); ++f)
    {
        const double wavenum = TWO_PI * freq(f) / c0 ;

        for (size_t n = 0; n < range.size(); ++n)
        {
            os << freq(f) << "," ;
            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 (size_t 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 ;

        // extremely 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 );
                BOOST_CHECK( detcoef >= 90.0 );
                break;
            case 100:
                BOOST_CHECK( abs(bias) <= 2.0 );
                BOOST_CHECK( dev <= 4.0 );
                BOOST_CHECK( detcoef >= 90.0 );
                break;
            case 10000:
                BOOST_CHECK( dev <= 6.0 );
                BOOST_CHECK( detcoef >= 75.0 );
                break;
            default:
                BOOST_CHECK( abs(bias) <= 1.0 );
                BOOST_CHECK( dev <= 4.0 );
                BOOST_CHECK( detcoef >= 90 );
        }

    } // 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 ;
    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
    wposition target(depth.size(), 1, degrees, src_lng, 0.0);
    for (size_t n = 0; n < target.size1(); ++n)
    {
        target.altitude(n, 0, depth(n));
    }

    proploss loss(freq, pos, de, az, time_step, &target);
    wave_queue wave( ocean, freq, pos, de, az, time_step, &target) ;
    wave.addEigenrayListener(&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 (size_t 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 (size_t 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 <= 5.0 );
    BOOST_CHECK( detcoef >= 80.0 );
}

BOOST_AUTO_TEST_CASE(proploss_limits)
{
    cout << "=== proploss_test: proploss_limits ===" << endl;
    const double c0 = 1500.0;
    const double src_lat = 45.0;
    const double src_lng = -45.0;
    const double src_alt = 0.0;
    const double depth = -1000 ;

    // initialize propagation model

    profile_model* profile = new profile_linear(c0);
    boundary_model* surface = new boundary_flat();
    boundary_model* bottom = new boundary_flat(depth);
    ocean_model ocean(surface, bottom, profile);

    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);
    seq_log freq(f0, 1.0, 1);  // 2000 Hz

    // try building a source above ocean surface

    wave_queue wave1( ocean, freq, pos, de, az, time_step ) ;
    BOOST_CHECK_CLOSE( wave1.source_pos().altitude(), -0.1, 1e-6 );

    // try building a source above ocean surface

    pos.altitude( depth-10.0 ) ;
    wave_queue wave2( ocean, freq, pos, de, az, time_step ) ;
    BOOST_CHECK_CLOSE( wave2.source_pos().altitude(), depth+0.1, 1e-6 );
}



BOOST_AUTO_TEST_SUITE_END()

Generated on 4 May 2015 for USML by  doxygen 1.6.1