FM / formulas

/usr/include/boost/geometry/formulas

Name	Size	Mode	Actions
andoyer_inverse.hpp	9909	0644	edit dl rm
area_formulas.hpp	21897	0644	edit dl rm
authalic_radius_sqr.hpp	2654	0644	edit dl rm
differential_quantities.hpp	10976	0644	edit dl rm
eccentricity_sqr.hpp	1959	0644	edit dl rm
flattening.hpp	1847	0644	edit dl rm
geographic.hpp	14251	0644	edit dl rm
gnomonic_intersection.hpp	4905	0644	edit dl rm
gnomonic_spheroid.hpp	3989	0644	edit dl rm
interpolate_point_spherical.hpp	3586	0644	edit dl rm
karney_direct.hpp	10063	0644	edit dl rm
karney_inverse.hpp	35130	0644	edit dl rm
mean_radius.hpp	1894	0644	edit dl rm
meridian_direct.hpp	5139	0644	edit dl rm
meridian_inverse.hpp	4782	0644	edit dl rm
meridian_segment.hpp	2009	0644	edit dl rm
quarter_meridian.hpp	3174	0644	edit dl rm
result_direct.hpp	879	0644	edit dl rm
result_inverse.hpp	903	0644	edit dl rm
sjoberg_intersection.hpp	43099	0644	edit dl rm
spherical.hpp	8788	0644	edit dl rm
thomas_direct.hpp	9256	0644	edit dl rm
thomas_inverse.hpp	7768	0644	edit dl rm
unit_spheroid.hpp	1199	0644	edit dl rm
vertex_latitude.hpp	4489	0644	edit dl rm
vertex_longitude.hpp	10613	0644	edit dl rm
vincenty_direct.hpp	7105	0644	edit dl rm
vincenty_inverse.hpp	8254	0644	edit dl rm

Edit: /usr/include/boost/geometry/formulas/gnomonic_spheroid.hpp (3989B)

// Boost.Geometry

// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// 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)

#ifndef BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP
#define BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP

#include <boost/geometry/core/radius.hpp>

#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>

#include <boost/geometry/formulas/andoyer_inverse.hpp>
#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/thomas_inverse.hpp>
#include <boost/geometry/formulas/vincenty_direct.hpp>
#include <boost/geometry/formulas/vincenty_inverse.hpp>

namespace boost { namespace geometry { namespace formula
{

/*!
\brief Gnomonic projection on spheroid (ellipsoid of revolution).
\author See
- Charles F.F Karney, Algorithms for geodesics, 2011
https://arxiv.org/pdf/1109.4448.pdf
*/
template <
    typename CT,
    template <typename, bool, bool, bool, bool ,bool> class Inverse,
    template <typename, bool, bool, bool, bool> class Direct
>
class gnomonic_spheroid
{
    typedef Inverse<CT, false, true, true, true, true> inverse_type;
    typedef typename inverse_type::result_type inverse_result;

typedef Direct<CT, false, false, true, true> direct_quantities_type;
    typedef Direct<CT, true, false, false, false> direct_coordinates_type;
    typedef typename direct_coordinates_type::result_type direct_result;

public:
    template <typename Spheroid>
    static inline bool forward(CT const& lon0, CT const& lat0,
                               CT const& lon, CT const& lat,
                               CT & x, CT & y,
                               Spheroid const& spheroid)
    {
        inverse_result i_res = inverse_type::apply(lon0, lat0, lon, lat, spheroid);
        CT const& m = i_res.reduced_length;
        CT const& M = i_res.geodesic_scale;

if (math::smaller_or_equals(M, CT(0)))
        {
            return false;
        }

CT rho = m / M;
        x = sin(i_res.azimuth) * rho;
        y = cos(i_res.azimuth) * rho;
        
        return true;
    }

template <typename Spheroid>
    static inline bool inverse(CT const& lon0, CT const& lat0,
                               CT const& x, CT const& y,
                               CT & lon, CT & lat,
                               Spheroid const& spheroid)
    {
        CT const a = get_radius<0>(spheroid);
        CT const ds_threshold = a * std::numeric_limits<CT>::epsilon(); // TODO: 0 for non-fundamental type
        
        CT const azimuth = atan2(x, y);
        CT const rho = math::sqrt(math::sqr(x) + math::sqr(y)); // use hypot?
        CT distance = a * atan(rho / a);

bool found = false;
        for (int i = 0 ; i < 10 ; ++i)
        {
            direct_result d_res = direct_quantities_type::apply(lon0, lat0, distance, azimuth, spheroid);
            CT const& m = d_res.reduced_length;
            CT const& M = d_res.geodesic_scale;

if (math::smaller_or_equals(M, CT(0)))
            {
                // found = false;
                return found;
            }
            
            CT const drho = m / M - rho; // rho = m / M
            CT const ds = drho * math::sqr(M); // drho/ds = 1/M^2
            distance -= ds;

// ds_threshold may be 0
            if (math::abs(ds) <= ds_threshold)
            {
                found = true;
                break;
            }
        }

if (found)
        {
            direct_result d_res = direct_coordinates_type::apply(lon0, lat0, distance, azimuth, spheroid);
            lon = d_res.lon2;
            lat = d_res.lat2;
        }

return found;
    }
};

}}} // namespace boost::geometry::formula

#endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP

// Boost.Geometry // Copyright (c) 2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // 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) #ifndef BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP #define BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP #include #include #include #include #include #include #include #include namespace boost { namespace geometry { namespace formula { /*! \brief Gnomonic projection on spheroid (ellipsoid of revolution). \author See - Charles F.F Karney, Algorithms for geodesics, 2011 https://arxiv.org/pdf/1109.4448.pdf */ template < typename CT, template class Inverse, template class Direct > class gnomonic_spheroid { typedef Inverse inverse_type; typedef typename inverse_type::result_type inverse_result; typedef Direct direct_quantities_type; typedef Direct direct_coordinates_type; typedef typename direct_coordinates_type::result_type direct_result; public: template static inline bool forward(CT const& lon0, CT const& lat0, CT const& lon, CT const& lat, CT & x, CT & y, Spheroid const& spheroid) { inverse_result i_res = inverse_type::apply(lon0, lat0, lon, lat, spheroid); CT const& m = i_res.reduced_length; CT const& M = i_res.geodesic_scale; if (math::smaller_or_equals(M, CT(0))) { return false; } CT rho = m / M; x = sin(i_res.azimuth) * rho; y = cos(i_res.azimuth) * rho; return true; } template static inline bool inverse(CT const& lon0, CT const& lat0, CT const& x, CT const& y, CT & lon, CT & lat, Spheroid const& spheroid) { CT const a = get_radius<0>(spheroid); CT const ds_threshold = a * std::numeric_limits::epsilon(); // TODO: 0 for non-fundamental type CT const azimuth = atan2(x, y); CT const rho = math::sqrt(math::sqr(x) + math::sqr(y)); // use hypot? CT distance = a * atan(rho / a); bool found = false; for (int i = 0 ; i < 10 ; ++i) { direct_result d_res = direct_quantities_type::apply(lon0, lat0, distance, azimuth, spheroid); CT const& m = d_res.reduced_length; CT const& M = d_res.geodesic_scale; if (math::smaller_or_equals(M, CT(0))) { // found = false; return found; } CT const drho = m / M - rho; // rho = m / M CT const ds = drho * math::sqr(M); // drho/ds = 1/M^2 distance -= ds; // ds_threshold may be 0 if (math::abs(ds) <= ds_threshold) { found = true; break; } } if (found) { direct_result d_res = direct_coordinates_type::apply(lon0, lat0, distance, azimuth, spheroid); lon = d_res.lon2; lat = d_res.lat2; } return found; } }; }}} // namespace boost::geometry::formula #endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP