From 48d5f65669e71c178ed5e9b797f8284d4c19da66 Mon Sep 17 00:00:00 2001 From: Barend Gehrels Date: Fri, 27 May 2011 23:02:58 +0000 Subject: [PATCH] Fixed ssf for spherical equatorial coordinate system (old version was for polar) [SVN r72238] --- .../geometry/strategies/spherical/ssf.hpp | 58 +++++++++------ .../strategies/strategy_transform.hpp | 70 +++++++++++++++++-- 2 files changed, 100 insertions(+), 28 deletions(-) diff --git a/include/boost/geometry/strategies/spherical/ssf.hpp b/include/boost/geometry/strategies/spherical/ssf.hpp index fdb41fc26..33ac8eea4 100644 --- a/include/boost/geometry/strategies/spherical/ssf.hpp +++ b/include/boost/geometry/strategies/spherical/ssf.hpp @@ -69,35 +69,35 @@ public : // Convenient shortcuts typedef coordinate_type ct; - ct const phi1 = get_as_radian<0>(p1); - ct const theta1 = get_as_radian<1>(p1); - ct const phi2 = get_as_radian<0>(p2); - ct const theta2 = get_as_radian<1>(p2); - ct const phi = get_as_radian<0>(p); - ct const theta = get_as_radian<1>(p); + ct const lambda1 = get_as_radian<0>(p1); + ct const delta1 = get_as_radian<1>(p1); + ct const lambda2 = get_as_radian<0>(p2); + ct const delta2 = get_as_radian<1>(p2); + ct const lambda = get_as_radian<0>(p); + ct const delta = get_as_radian<1>(p); // Create temporary points (vectors) on unit a sphere - ct const sin_theta1 = sin(theta1); - ct const c1x = sin_theta1 * cos(phi1); - ct const c1y = sin_theta1 * sin(phi1); - ct const c1z = cos(theta1); + ct const cos_delta1 = cos(delta1); + ct const c1x = cos_delta1 * cos(lambda1); + ct const c1y = cos_delta1 * sin(lambda1); + ct const c1z = sin(delta1); - ct const sin_theta2 = sin(theta2); - ct const c2x = sin_theta2 * cos(phi2); - ct const c2y = sin_theta2 * sin(phi2); - ct const c2z = cos(theta2); + ct const cos_delta2 = cos(delta2); + ct const c2x = cos_delta2 * cos(lambda2); + ct const c2y = cos_delta2 * sin(lambda2); + ct const c2z = sin(delta2); // (Third point is converted directly) - ct const sin_theta = sin(theta); + ct const cos_delta = cos(delta); // Apply the "Spherical Side Formula" as presented on my blog ct const dist - = (c1y * c2z - c1z * c2y) * sin_theta * cos(phi) - + (c1z * c2x - c1x * c2z) * sin_theta * sin(phi) - + (c1x * c2y - c1y * c2x) * cos(theta); - - return dist < 0 ? 1 - : dist > 0 ? -1 + = (c1y * c2z - c1z * c2y) * cos_delta * cos(lambda) + + (c1z * c2x - c1x * c2z) * cos_delta * sin(lambda) + + (c1x * c2y - c1y * c2x) * sin(delta); + + return dist > 0 ? 1 + : dist < 0 ? -1 : 0; } }; @@ -105,6 +105,22 @@ public : }} // namespace strategy::side +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +template +struct strategy_side +{ + typedef strategy::side::spherical_side_formula type; +}; + +template +struct strategy_side +{ + typedef strategy::side::spherical_side_formula type; +}; + +#endif + + }} // namespace boost::geometry diff --git a/include/boost/geometry/strategies/strategy_transform.hpp b/include/boost/geometry/strategies/strategy_transform.hpp index dd6ead146..34e19fc77 100644 --- a/include/boost/geometry/strategies/strategy_transform.hpp +++ b/include/boost/geometry/strategies/strategy_transform.hpp @@ -155,29 +155,48 @@ namespace detail /// Helper function for conversion, phi/theta are in radians template - inline void spherical_to_cartesian(T phi, T theta, R r, P& p) + inline void spherical_polar_to_cartesian(T phi, T theta, R r, P& p) { assert_dimension(); // http://en.wikipedia.org/wiki/List_of_canonical_coordinate_transformations#From_spherical_coordinates // http://www.vias.org/comp_geometry/math_coord_convert_3d.htm // https://moodle.polymtl.ca/file.php/1183/Autres_Documents/Derivation_for_Spherical_Co-ordinates.pdf + // http://en.citizendium.org/wiki/Spherical_polar_coordinates // Phi = first, theta is second, r is third, see documentation on cs::spherical // (calculations are splitted to implement ttmath) T r_sin_theta = r; + T r_cos_theta = r; r_sin_theta *= sin(theta); + r_cos_theta *= cos(theta); set<0>(p, r_sin_theta * cos(phi)); set<1>(p, r_sin_theta * sin(phi)); - - T r_cos_theta = r; - r_cos_theta *= cos(theta); - set<2>(p, r_cos_theta); } + + /// Helper function for conversion, lambda/delta (lon lat) are in radians + template + inline void spherical_equatorial_to_cartesian(T lambda, T delta, R r, P& p) + { + assert_dimension(); + + // http://mathworld.wolfram.com/GreatCircle.html + // http://www.spenvis.oma.be/help/background/coortran/coortran.html WRONG + + T r_cos_delta = r; + T r_sin_delta = r; + r_cos_delta *= cos(delta); + r_sin_delta *= sin(delta); + + set<0>(p, r_cos_delta * cos(lambda)); + set<1>(p, r_cos_delta * sin(lambda)); + set<2>(p, r_sin_delta); + } + /// Helper function for conversion template @@ -238,11 +257,23 @@ struct from_spherical_polar_2_to_cartesian_3 inline bool apply(P1 const& p1, P2& p2) const { assert_dimension(); - detail::spherical_to_cartesian(get_as_radian<0>(p1), get_as_radian<1>(p1), 1.0, p2); + detail::spherical_polar_to_cartesian(get_as_radian<0>(p1), get_as_radian<1>(p1), 1.0, p2); return true; } }; +template +struct from_spherical_equatorial_2_to_cartesian_3 +{ + inline bool apply(P1 const& p1, P2& p2) const + { + assert_dimension(); + detail::spherical_equatorial_to_cartesian(get_as_radian<0>(p1), get_as_radian<1>(p1), 1.0, p2); + return true; + } +}; + + /*! \brief Transformation strategy for 3D spherical (phi,theta,r) to 3D cartesian (x,y,z) \ingroup transform @@ -255,12 +286,25 @@ struct from_spherical_polar_3_to_cartesian_3 inline bool apply(P1 const& p1, P2& p2) const { assert_dimension(); - detail::spherical_to_cartesian( + detail::spherical_polar_to_cartesian( get_as_radian<0>(p1), get_as_radian<1>(p1), get<2>(p1), p2); return true; } }; +template +struct from_spherical_equatorial_3_to_cartesian_3 +{ + inline bool apply(P1 const& p1, P2& p2) const + { + assert_dimension(); + detail::spherical_equatorial_to_cartesian( + get_as_radian<0>(p1), get_as_radian<1>(p1), get<2>(p1), p2); + return true; + } +}; + + /*! \brief Transformation strategy for 3D cartesian (x,y,z) to 2D spherical (phi,theta) \details on Unit sphere @@ -358,6 +402,18 @@ struct default_strategy type; }; +template +struct default_strategy +{ + typedef from_spherical_equatorial_2_to_cartesian_3 type; +}; + +template +struct default_strategy +{ + typedef from_spherical_equatorial_3_to_cartesian_3 type; +}; + /// Specialization to transform from XYZ to unit sphere(phi,theta) template struct default_strategy