[ceph.git] / ceph / src / boost / libs / geometry / test / strategies / vincenty.cpp

// Boost.Geometry (aka GGL, Generic Geometry Library)
// Unit Test

// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.

// This file was modified by Oracle on 2014-2021.
// Modifications copyright (c) 2014-2021 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.

// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)


#include <geometry_test_common.hpp>

#include <boost/concept_check.hpp>

#include <boost/geometry/algorithms/assign.hpp>
#include <boost/geometry/algorithms/distance.hpp>
#include <boost/geometry/formulas/vincenty_inverse.hpp>
#include <boost/geometry/formulas/vincenty_direct.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/srs/spheroid.hpp>
#include <boost/geometry/strategies/concepts/distance_concept.hpp>
#include <boost/geometry/strategies/geographic/distance_vincenty.hpp>
#include <boost/geometry/strategies/geographic/side_vincenty.hpp>

#include <test_common/test_point.hpp>

template <typename T>
void normalize_deg(T & deg)
{
    while ( deg > T(180) )
        deg -= T(360);
    while ( deg <= T(-180) )
        deg += T(360);
}

template <typename T>
T difference_deg(T const& a1, T const& a2)
{
    T d = a1 - a2;
    normalize_deg(d);
    return d;
}

template <typename T>
void check_deg(std::string const& name, T const& a1, T const& a2, T const& percent, T const& error)
{
    T diff = bg::math::abs(difference_deg(a1, a2));
    
    if ( bg::math::equals(a1, T(0)) || bg::math::equals(a2, T(0)) )
    {
        if ( diff > error )
        {
            BOOST_ERROR(name << " - the difference {" << diff << "} between {" << a1 << "} and {" << a2 << "} exceeds {" << error << "}");
        }
    }
    else
    {
        T greater = (std::max)(bg::math::abs(a1), bg::math::abs(a2));

        if ( diff > greater * percent / T(100) )
        {
            BOOST_ERROR(name << " the difference {" << diff << "} between {" << a1 << "} and {" << a2 << "} exceeds {" << percent << "}%");
        }
    }
}

double azimuth(double deg, double min, double sec)
{
    min = fabs(min);
    sec = fabs(sec);

    if ( deg < 0 )
    {
        min = -min;
        sec = -sec;
    }

    return deg + min/60.0 + sec/3600.0;
}

double azimuth(double deg, double min)
{
    return azimuth(deg, min, 0.0);
}

template <typename P>
bool non_precise_ct()
{
    using ct = typename bg::coordinate_type<P>::type;
    return std::is_integral<ct>::value || std::is_floating_point<ct>::value;
}

template <typename P1, typename P2, typename Spheroid>
void test_vincenty(double lon1, double lat1, double lon2, double lat2,
                   double expected_distance,
                   double expected_azimuth_12,
                   double /*expected_azimuth_21*/,
                   Spheroid const& spheroid)
{
    typedef typename bg::promote_floating_point
        <
            typename bg::select_calculation_type<P1, P2, void>::type
        >::type calc_t;

    calc_t tolerance = non_precise_ct<P1>() || non_precise_ct<P2>() ?
                       5.0 : 0.001;
    calc_t error = non_precise_ct<P1>() || non_precise_ct<P2>() ?
                   1e-5 : 1e-12;

    // formula
    {
        double const d2r = bg::math::d2r<double>();
        double const r2d = bg::math::r2d<double>();

        typedef bg::formula::vincenty_inverse<calc_t, true, true> inverse_formula;
        typename inverse_formula::result_type
            result_i = inverse_formula::apply(lon1 * d2r,
                                             lat1 * d2r,
                                             lon2 * d2r,
                                             lat2 * d2r,
                                             spheroid);
        calc_t dist = result_i.distance;
        calc_t az12 = result_i.azimuth;
        //calc_t az21 = vi.azimuth21();

        calc_t az12_deg = az12 * r2d;
        //calc_t az21_deg = az21 * r2d;
        
        BOOST_CHECK_CLOSE(dist, calc_t(expected_distance), tolerance);
        check_deg("az12_deg", az12_deg, calc_t(expected_azimuth_12), tolerance, error);
        //check_deg("az21_deg", az21_deg, calc_t(expected_azimuth_21), tolerance, error);

        typedef bg::formula::vincenty_direct<calc_t> direct_formula;
        typename direct_formula::result_type
            result_d = direct_formula::apply(lon1 * d2r,
                                             lat1 * d2r,
                                             dist,
                                             az12,
                                             spheroid);
        calc_t direct_lon2 = result_d.lon2;
        calc_t direct_lat2 = result_d.lat2;
        //calc_t direct_az21 = vd.azimuth21();

        calc_t direct_lon2_deg = direct_lon2 * r2d;
        calc_t direct_lat2_deg = direct_lat2 * r2d;
        //calc_t direct_az21_deg = direct_az21 * r2d;
        
        check_deg("direct_lon2_deg", direct_lon2_deg, calc_t(lon2), tolerance, error);
        check_deg("direct_lat2_deg", direct_lat2_deg, calc_t(lat2), tolerance, error);
        //check_deg("direct_az21_deg", direct_az21_deg, az21_deg, tolerance, error);
    }

    // distance strategies
    {
        typedef bg::strategy::distance::vincenty<Spheroid> vincenty_type;
        typedef bg::strategy::distance::geographic<bg::strategy::vincenty, Spheroid> geographic_type;

        BOOST_CONCEPT_ASSERT(
            (
                bg::concepts::PointDistanceStrategy<vincenty_type, P1, P2>)
            );

        vincenty_type vincenty(spheroid);
        geographic_type geographic(spheroid);
        typedef typename bg::strategy::distance::services::return_type<vincenty_type, P1, P2>::type return_type;

        P1 p1;
        P2 p2;

        bg::assign_values(p1, lon1, lat1);
        bg::assign_values(p2, lon2, lat2);
        
        BOOST_CHECK_CLOSE(vincenty.apply(p1, p2), return_type(expected_distance), tolerance);
        BOOST_CHECK_CLOSE(geographic.apply(p1, p2), return_type(expected_distance), tolerance);
        BOOST_CHECK_CLOSE(bg::distance(p1, p2, vincenty), return_type(expected_distance), tolerance);
    }
}

template <typename P1, typename P2>
void test_vincenty(double lon1, double lat1, double lon2, double lat2,
                   double expected_distance,
                   double expected_azimuth_12,
                   double expected_azimuth_21)
{
    test_vincenty<P1, P2>(lon1, lat1, lon2, lat2,
                          expected_distance, expected_azimuth_12, expected_azimuth_21,
                          bg::srs::spheroid<double>());
}

template <typename PS, typename P>
void test_side(double lon1, double lat1,
               double lon2, double lat2,
               double lon, double lat,
               int expected_side)
{
    // Set radius type, but for integer coordinates we want to have floating point radius type
    typedef typename bg::promote_floating_point
        <
            typename bg::coordinate_type<PS>::type
        >::type rtype;

    typedef bg::srs::spheroid<rtype> stype;

    typedef bg::strategy::side::vincenty<stype> strategy_type;
    typedef bg::strategy::side::geographic<bg::strategy::vincenty, stype> strategy2_type;

    strategy_type strategy;
    strategy2_type strategy2;

    PS p1, p2;
    P p;

    bg::assign_values(p1, lon1, lat1);
    bg::assign_values(p2, lon2, lat2);
    bg::assign_values(p, lon, lat);

    int side = strategy.apply(p1, p2, p);
    int side2 = strategy2.apply(p1, p2, p);

    BOOST_CHECK_EQUAL(side, expected_side);
    BOOST_CHECK_EQUAL(side2, expected_side);
}

template <typename P1, typename P2>
void test_all()
{
    // See:
    //  - http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp
    //  - http://www.ga.gov.au/geodesy/datums/vincenty_direct.jsp
    // Values in the comments below was calculated using the above pages
    // in some cases distances may be different, previously used values was left

    // use km
    double gda_a = 6378.1370;
    double gda_f = 1.0 / 298.25722210;
    double gda_b = gda_a * ( 1.0 - gda_f );
    bg::srs::spheroid<double> gda_spheroid(gda_a, gda_b);

    // Test fractional coordinates only for non-integral types
    if ( BOOST_GEOMETRY_CONDITION(
            ! std::is_integral<typename bg::coordinate_type<P1>::type>::value
         && ! std::is_integral<typename bg::coordinate_type<P2>::type>::value ) )
    {
        // Flinders Peak -> Buninyong
        test_vincenty<P1, P2>(azimuth(144,25,29.52440), azimuth(-37,57,3.72030),
                              azimuth(143,55,35.38390), azimuth(-37,39,10.15610),
                              54.972271, azimuth(306,52,5.37), azimuth(127,10,25.07),
                              gda_spheroid);
    }

    // Lodz -> Trondheim
    test_vincenty<P1, P2>(azimuth(19,28), azimuth(51,47),
                          azimuth(10,21), azimuth(63,23),
                          1399.032724, azimuth(340,54,25.14), azimuth(153,10,0.19),
                          gda_spheroid);
    // London -> New York
    test_vincenty<P1, P2>(azimuth(0,7,39), azimuth(51,30,26),
                          azimuth(-74,0,21), azimuth(40,42,46),
                          5602.044851, azimuth(288,31,36.82), azimuth(51,10,33.43),
                          gda_spheroid);

    // Shanghai -> San Francisco
    test_vincenty<P1, P2>(azimuth(121,30), azimuth(31,12),
                          azimuth(-122,25), azimuth(37,47),
                          9899.698550, azimuth(45,12,44.76), azimuth(309,50,20.88),
                          gda_spheroid);

    test_vincenty<P1, P2>(0, 0, 0, 50, 5540.847042, 0, 180, gda_spheroid); // N
    test_vincenty<P1, P2>(0, 0, 0, -50, 5540.847042, 180, 0, gda_spheroid); // S
    test_vincenty<P1, P2>(0, 0, 50, 0, 5565.974540, 90, -90, gda_spheroid); // E
    test_vincenty<P1, P2>(0, 0, -50, 0, 5565.974540, -90, 90, gda_spheroid); // W
    
    test_vincenty<P1, P2>(0, 0, 50, 50, 7284.879297, azimuth(32,51,55.87), azimuth(237,24,50.12), gda_spheroid); // NE
    
    // The original distance values, azimuths calculated using the web form mentioned above
    // Using default spheroid units (meters)
    test_vincenty<P1, P2>(0, 89, 1, 80, 1005153.5769, azimuth(178,53,23.85), azimuth(359,53,18.35)); // sub-polar
    test_vincenty<P1, P2>(4, 52, 4, 52, 0.0, 0, 0); // no point difference
    test_vincenty<P1, P2>(4, 52, 3, 40, 1336039.890, azimuth(183,41,29.08), azimuth(2,58,5.13)); // normal case

    test_side<P1, P2>(0, 0, 0, 1, 0, 2, 0);
    test_side<P1, P2>(0, 0, 0, 1, 0, -2, 0);
    test_side<P1, P2>(10, 0, 10, 1, 10, 2, 0);
    test_side<P1, P2>(10, 0, 10, -1, 10, 2, 0);

    test_side<P1, P2>(10, 0, 10, 1, 0, 2, 1); // left
    test_side<P1, P2>(10, 0, 10, -1, 0, 2, -1); // right

    test_side<P1, P2>(-10, -10, 10, 10, 10, 0, -1); // right
    test_side<P1, P2>(-10, -10, 10, 10, -10, 0, 1); // left
    test_side<P1, P2>(170, -10, -170, 10, -170, 0, -1); // right
    test_side<P1, P2>(170, -10, -170, 10, 170, 0, 1); // left
}

template <typename P>
void test_all()
{
    test_all<P, P>();
}

int test_main(int, char* [])
{
    //test_all<float[2]>();
    //test_all<double[2]>();
    test_all<bg::model::point<double, 2, bg::cs::geographic<bg::degree> > >();
    test_all<bg::model::point<float, 2, bg::cs::geographic<bg::degree> > >();
    test_all<bg::model::point<int, 2, bg::cs::geographic<bg::degree> > >();

    return 0;
}
Commit	Line	Data
7c673cae FG	1	// Boost.Geometry (aka GGL, Generic Geometry Library)
	2	// Unit Test
	3
	4	// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
	5	// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
	6	// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
	7
1e59de90 TL	8	// This file was modified by Oracle on 2014-2021.
1e59de90 TL	9	// Modifications copyright (c) 2014-2021 Oracle and/or its affiliates.
7c673cae FG	10	// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
	11
	12	// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
	13	// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
	14
	15	// Use, modification and distribution is subject to the Boost Software License,
	16	// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
	17	// http://www.boost.org/LICENSE_1_0.txt)
	18
	19
	20	#include <geometry_test_common.hpp>
	21
	22	#include <boost/concept_check.hpp>
	23
11fdf7f2 TL	24	#include <boost/geometry/algorithms/assign.hpp>
11fdf7f2 TL	25	#include <boost/geometry/algorithms/distance.hpp>
b32b8144 FG	26	#include <boost/geometry/formulas/vincenty_inverse.hpp>
b32b8144 FG	27	#include <boost/geometry/formulas/vincenty_direct.hpp>
7c673cae	28	#include <boost/geometry/geometries/point.hpp>
11fdf7f2 TL	29	#include <boost/geometry/srs/spheroid.hpp>
	30	#include <boost/geometry/strategies/concepts/distance_concept.hpp>
	31	#include <boost/geometry/strategies/geographic/distance_vincenty.hpp>
	32	#include <boost/geometry/strategies/geographic/side_vincenty.hpp>
	33
7c673cae FG	34	#include <test_common/test_point.hpp>
7c673cae FG	35
7c673cae FG	36	template <typename T>
	37	void normalize_deg(T & deg)
	38	{
	39	while ( deg > T(180) )
	40	deg -= T(360);
	41	while ( deg <= T(-180) )
	42	deg += T(360);
	43	}
	44
	45	template <typename T>
	46	T difference_deg(T const& a1, T const& a2)
	47	{
	48	T d = a1 - a2;
	49	normalize_deg(d);
	50	return d;
	51	}
	52
	53	template <typename T>
	54	void check_deg(std::string const& name, T const& a1, T const& a2, T const& percent, T const& error)
	55	{
	56	T diff = bg::math::abs(difference_deg(a1, a2));
	57
	58	if ( bg::math::equals(a1, T(0)) \|\| bg::math::equals(a2, T(0)) )
	59	{
	60	if ( diff > error )
	61	{
	62	BOOST_ERROR(name << " - the difference {" << diff << "} between {" << a1 << "} and {" << a2 << "} exceeds {" << error << "}");
	63	}
	64	}
	65	else
	66	{
	67	T greater = (std::max)(bg::math::abs(a1), bg::math::abs(a2));
	68
	69	if ( diff > greater * percent / T(100) )
	70	{
	71	BOOST_ERROR(name << " the difference {" << diff << "} between {" << a1 << "} and {" << a2 << "} exceeds {" << percent << "}%");
	72	}
	73	}
	74	}
	75
	76	double azimuth(double deg, double min, double sec)
	77	{
	78	min = fabs(min);
	79	sec = fabs(sec);
	80
	81	if ( deg < 0 )
	82	{
	83	min = -min;
	84	sec = -sec;
	85	}
	86
	87	return deg + min/60.0 + sec/3600.0;
	88	}
	89
	90	double azimuth(double deg, double min)
	91	{
	92	return azimuth(deg, min, 0.0);
	93	}
	94
	95	template <typename P>
	96	bool non_precise_ct()
	97	{
1e59de90 TL	98	using ct = typename bg::coordinate_type<P>::type;
1e59de90 TL	99	return std::is_integral<ct>::value \|\| std::is_floating_point<ct>::value;
7c673cae FG	100	}
	101
	102	template <typename P1, typename P2, typename Spheroid>
	103	void test_vincenty(double lon1, double lat1, double lon2, double lat2,
	104	double expected_distance,
	105	double expected_azimuth_12,
	106	double /expected_azimuth_21/,
	107	Spheroid const& spheroid)
	108	{
	109	typedef typename bg::promote_floating_point
	110	<
	111	typename bg::select_calculation_type<P1, P2, void>::type
	112	>::type calc_t;
	113
	114	calc_t tolerance = non_precise_ct<P1>() \|\| non_precise_ct<P2>() ?
	115	5.0 : 0.001;
	116	calc_t error = non_precise_ct<P1>() \|\| non_precise_ct<P2>() ?
	117	1e-5 : 1e-12;
	118
	119	// formula
	120	{
	121	double const d2r = bg::math::d2r<double>();
	122	double const r2d = bg::math::r2d<double>();
	123
b32b8144	124	typedef bg::formula::vincenty_inverse<calc_t, true, true> inverse_formula;
7c673cae FG	125	typename inverse_formula::result_type
	126	result_i = inverse_formula::apply(lon1 * d2r,
	127	lat1 * d2r,
	128	lon2 * d2r,
	129	lat2 * d2r,
	130	spheroid);
	131	calc_t dist = result_i.distance;
	132	calc_t az12 = result_i.azimuth;
	133	//calc_t az21 = vi.azimuth21();
	134
	135	calc_t az12_deg = az12 * r2d;
	136	//calc_t az21_deg = az21 * r2d;
	137
	138	BOOST_CHECK_CLOSE(dist, calc_t(expected_distance), tolerance);
	139	check_deg("az12_deg", az12_deg, calc_t(expected_azimuth_12), tolerance, error);
	140	//check_deg("az21_deg", az21_deg, calc_t(expected_azimuth_21), tolerance, error);
	141
b32b8144	142	typedef bg::formula::vincenty_direct<calc_t> direct_formula;
7c673cae FG	143	typename direct_formula::result_type
	144	result_d = direct_formula::apply(lon1 * d2r,
	145	lat1 * d2r,
	146	dist,
	147	az12,
	148	spheroid);
	149	calc_t direct_lon2 = result_d.lon2;
	150	calc_t direct_lat2 = result_d.lat2;
	151	//calc_t direct_az21 = vd.azimuth21();
	152
	153	calc_t direct_lon2_deg = direct_lon2 * r2d;
	154	calc_t direct_lat2_deg = direct_lat2 * r2d;
	155	//calc_t direct_az21_deg = direct_az21 * r2d;
	156
	157	check_deg("direct_lon2_deg", direct_lon2_deg, calc_t(lon2), tolerance, error);
	158	check_deg("direct_lat2_deg", direct_lat2_deg, calc_t(lat2), tolerance, error);
	159	//check_deg("direct_az21_deg", direct_az21_deg, az21_deg, tolerance, error);
	160	}
	161
11fdf7f2	162	// distance strategies
7c673cae FG	163	{
7c673cae FG	164	typedef bg::strategy::distance::vincenty<Spheroid> vincenty_type;
11fdf7f2	165	typedef bg::strategy::distance::geographic<bg::strategy::vincenty, Spheroid> geographic_type;
7c673cae FG	166
	167	BOOST_CONCEPT_ASSERT(
	168	(
	169	bg::concepts::PointDistanceStrategy<vincenty_type, P1, P2>)
	170	);
	171
	172	vincenty_type vincenty(spheroid);
11fdf7f2	173	geographic_type geographic(spheroid);
7c673cae FG	174	typedef typename bg::strategy::distance::services::return_type<vincenty_type, P1, P2>::type return_type;
	175
	176	P1 p1;
	177	P2 p2;
	178
	179	bg::assign_values(p1, lon1, lat1);
	180	bg::assign_values(p2, lon2, lat2);
	181
	182	BOOST_CHECK_CLOSE(vincenty.apply(p1, p2), return_type(expected_distance), tolerance);
11fdf7f2	183	BOOST_CHECK_CLOSE(geographic.apply(p1, p2), return_type(expected_distance), tolerance);
7c673cae FG	184	BOOST_CHECK_CLOSE(bg::distance(p1, p2, vincenty), return_type(expected_distance), tolerance);
	185	}
	186	}
	187
	188	template <typename P1, typename P2>
	189	void test_vincenty(double lon1, double lat1, double lon2, double lat2,
	190	double expected_distance,
	191	double expected_azimuth_12,
	192	double expected_azimuth_21)
	193	{
	194	test_vincenty<P1, P2>(lon1, lat1, lon2, lat2,
	195	expected_distance, expected_azimuth_12, expected_azimuth_21,
	196	bg::srs::spheroid<double>());
	197	}
	198
	199	template <typename PS, typename P>
	200	void test_side(double lon1, double lat1,
	201	double lon2, double lat2,
	202	double lon, double lat,
	203	int expected_side)
	204	{
	205	// Set radius type, but for integer coordinates we want to have floating point radius type
	206	typedef typename bg::promote_floating_point
	207	<
	208	typename bg::coordinate_type<PS>::type
	209	>::type rtype;
	210
	211	typedef bg::srs::spheroid<rtype> stype;
	212
	213	typedef bg::strategy::side::vincenty<stype> strategy_type;
11fdf7f2	214	typedef bg::strategy::side::geographic<bg::strategy::vincenty, stype> strategy2_type;
7c673cae FG	215
7c673cae FG	216	strategy_type strategy;
11fdf7f2	217	strategy2_type strategy2;
7c673cae FG	218
	219	PS p1, p2;
	220	P p;
	221
	222	bg::assign_values(p1, lon1, lat1);
	223	bg::assign_values(p2, lon2, lat2);
	224	bg::assign_values(p, lon, lat);
	225
	226	int side = strategy.apply(p1, p2, p);
11fdf7f2	227	int side2 = strategy2.apply(p1, p2, p);
7c673cae FG	228
7c673cae FG	229	BOOST_CHECK_EQUAL(side, expected_side);
11fdf7f2	230	BOOST_CHECK_EQUAL(side2, expected_side);
7c673cae FG	231	}
	232
	233	template <typename P1, typename P2>
	234	void test_all()
	235	{
	236	// See:
	237	// - http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp
	238	// - http://www.ga.gov.au/geodesy/datums/vincenty_direct.jsp
	239	// Values in the comments below was calculated using the above pages
	240	// in some cases distances may be different, previously used values was left
	241
	242	// use km
	243	double gda_a = 6378.1370;
	244	double gda_f = 1.0 / 298.25722210;
	245	double gda_b = gda_a * ( 1.0 - gda_f );
	246	bg::srs::spheroid<double> gda_spheroid(gda_a, gda_b);
	247
	248	// Test fractional coordinates only for non-integral types
	249	if ( BOOST_GEOMETRY_CONDITION(
1e59de90 TL	250	! std::is_integral<typename bg::coordinate_type<P1>::type>::value
1e59de90 TL	251	&& ! std::is_integral<typename bg::coordinate_type<P2>::type>::value ) )
7c673cae FG	252	{
	253	// Flinders Peak -> Buninyong
	254	test_vincenty<P1, P2>(azimuth(144,25,29.52440), azimuth(-37,57,3.72030),
	255	azimuth(143,55,35.38390), azimuth(-37,39,10.15610),
	256	54.972271, azimuth(306,52,5.37), azimuth(127,10,25.07),
	257	gda_spheroid);
	258	}
	259
	260	// Lodz -> Trondheim
	261	test_vincenty<P1, P2>(azimuth(19,28), azimuth(51,47),
	262	azimuth(10,21), azimuth(63,23),
	263	1399.032724, azimuth(340,54,25.14), azimuth(153,10,0.19),
	264	gda_spheroid);
	265	// London -> New York
	266	test_vincenty<P1, P2>(azimuth(0,7,39), azimuth(51,30,26),
	267	azimuth(-74,0,21), azimuth(40,42,46),
	268	5602.044851, azimuth(288,31,36.82), azimuth(51,10,33.43),
	269	gda_spheroid);
	270
	271	// Shanghai -> San Francisco
	272	test_vincenty<P1, P2>(azimuth(121,30), azimuth(31,12),
	273	azimuth(-122,25), azimuth(37,47),
	274	9899.698550, azimuth(45,12,44.76), azimuth(309,50,20.88),
	275	gda_spheroid);
	276
	277	test_vincenty<P1, P2>(0, 0, 0, 50, 5540.847042, 0, 180, gda_spheroid); // N
	278	test_vincenty<P1, P2>(0, 0, 0, -50, 5540.847042, 180, 0, gda_spheroid); // S
b32b8144	279	test_vincenty<P1, P2>(0, 0, 50, 0, 5565.974540, 90, -90, gda_spheroid); // E
7c673cae FG	280	test_vincenty<P1, P2>(0, 0, -50, 0, 5565.974540, -90, 90, gda_spheroid); // W
	281
	282	test_vincenty<P1, P2>(0, 0, 50, 50, 7284.879297, azimuth(32,51,55.87), azimuth(237,24,50.12), gda_spheroid); // NE
	283
	284	// The original distance values, azimuths calculated using the web form mentioned above
	285	// Using default spheroid units (meters)
	286	test_vincenty<P1, P2>(0, 89, 1, 80, 1005153.5769, azimuth(178,53,23.85), azimuth(359,53,18.35)); // sub-polar
	287	test_vincenty<P1, P2>(4, 52, 4, 52, 0.0, 0, 0); // no point difference
	288	test_vincenty<P1, P2>(4, 52, 3, 40, 1336039.890, azimuth(183,41,29.08), azimuth(2,58,5.13)); // normal case
	289
	290	test_side<P1, P2>(0, 0, 0, 1, 0, 2, 0);
	291	test_side<P1, P2>(0, 0, 0, 1, 0, -2, 0);
	292	test_side<P1, P2>(10, 0, 10, 1, 10, 2, 0);
	293	test_side<P1, P2>(10, 0, 10, -1, 10, 2, 0);
	294
	295	test_side<P1, P2>(10, 0, 10, 1, 0, 2, 1); // left
	296	test_side<P1, P2>(10, 0, 10, -1, 0, 2, -1); // right
	297
	298	test_side<P1, P2>(-10, -10, 10, 10, 10, 0, -1); // right
	299	test_side<P1, P2>(-10, -10, 10, 10, -10, 0, 1); // left
	300	test_side<P1, P2>(170, -10, -170, 10, -170, 0, -1); // right
	301	test_side<P1, P2>(170, -10, -170, 10, 170, 0, 1); // left
	302	}
	303
	304	template <typename P>
	305	void test_all()
	306	{
	307	test_all<P, P>();
	308	}
	309
	310	int test_main(int, char* [])
	311	{
	312	//test_all<float[2]>();
	313	//test_all<double[2]>();
	314	test_all<bg::model::point<double, 2, bg::cs::geographic<bg::degree> > >();
	315	test_all<bg::model::point<float, 2, bg::cs::geographic<bg::degree> > >();
	316	test_all<bg::model::point<int, 2, bg::cs::geographic<bg::degree> > >();
	317
7c673cae FG	318	return 0;
7c673cae FG	319	}