|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r72238 - in trunk/boost/geometry/strategies: . spherical
From: barend.gehrels_at_[hidden]
Date: 2011-05-27 19:02:59
Author: barendgehrels
Date: 2011-05-27 19:02:58 EDT (Fri, 27 May 2011)
New Revision: 72238
URL: http://svn.boost.org/trac/boost/changeset/72238
Log:
Fixed ssf for spherical equatorial coordinate system (old version was for polar)
Text files modified:
trunk/boost/geometry/strategies/spherical/ssf.hpp | 60 +++++++++++++++++++++------------
trunk/boost/geometry/strategies/strategy_transform.hpp | 70 ++++++++++++++++++++++++++++++++++++----
2 files changed, 101 insertions(+), 29 deletions(-)
Modified: trunk/boost/geometry/strategies/spherical/ssf.hpp
==============================================================================
--- trunk/boost/geometry/strategies/spherical/ssf.hpp (original)
+++ trunk/boost/geometry/strategies/spherical/ssf.hpp 2011-05-27 19:02:58 EDT (Fri, 27 May 2011)
@@ -69,35 +69,35 @@
// 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 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_delta1 = cos(delta1);
+ ct const c1x = cos_delta1 * cos(lambda1);
+ ct const c1y = cos_delta1 * sin(lambda1);
+ ct const c1z = sin(delta1);
+
+ 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 @@
}} // namespace strategy::side
+#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
+template <typename CalculationType>
+struct strategy_side<spherical_polar_tag, CalculationType>
+{
+ typedef strategy::side::spherical_side_formula<CalculationType> type;
+};
+
+template <typename CalculationType>
+struct strategy_side<geographic_tag, CalculationType>
+{
+ typedef strategy::side::spherical_side_formula<CalculationType> type;
+};
+
+#endif
+
+
}} // namespace boost::geometry
Modified: trunk/boost/geometry/strategies/strategy_transform.hpp
==============================================================================
--- trunk/boost/geometry/strategies/strategy_transform.hpp (original)
+++ trunk/boost/geometry/strategies/strategy_transform.hpp 2011-05-27 19:02:58 EDT (Fri, 27 May 2011)
@@ -155,29 +155,48 @@
/// Helper function for conversion, phi/theta are in radians
template <typename P, typename T, typename R>
- 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<P, 3>();
// 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 <typename P, typename T, typename R>
+ inline void spherical_equatorial_to_cartesian(T lambda, T delta, R r, P& p)
+ {
+ assert_dimension<P, 3>();
+
+ // 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 <typename P, typename T>
@@ -238,11 +257,23 @@
inline bool apply(P1 const& p1, P2& p2) const
{
assert_dimension<P1, 2>();
- 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 <typename P1, typename P2>
+struct from_spherical_equatorial_2_to_cartesian_3
+{
+ inline bool apply(P1 const& p1, P2& p2) const
+ {
+ assert_dimension<P1, 2>();
+ 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 @@
inline bool apply(P1 const& p1, P2& p2) const
{
assert_dimension<P1, 3>();
- 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 <typename P1, typename P2>
+struct from_spherical_equatorial_3_to_cartesian_3
+{
+ inline bool apply(P1 const& p1, P2& p2) const
+ {
+ assert_dimension<P1, 3>();
+ 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 @@
typedef from_spherical_polar_3_to_cartesian_3<P1, P2> type;
};
+template <typename CoordSys1, typename CoordSys2, typename P1, typename P2>
+struct default_strategy<spherical_equatorial_tag, cartesian_tag, CoordSys1, CoordSys2, 2, 3, P1, P2>
+{
+ typedef from_spherical_equatorial_2_to_cartesian_3<P1, P2> type;
+};
+
+template <typename CoordSys1, typename CoordSys2, typename P1, typename P2>
+struct default_strategy<spherical_equatorial_tag, cartesian_tag, CoordSys1, CoordSys2, 3, 3, P1, P2>
+{
+ typedef from_spherical_equatorial_3_to_cartesian_3<P1, P2> type;
+};
+
/// Specialization to transform from XYZ to unit sphere(phi,theta)
template <typename CoordSys1, typename CoordSys2, typename P1, typename P2>
struct default_strategy<cartesian_tag, spherical_polar_tag, CoordSys1, CoordSys2, 3, 2, P1, P2>
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk