boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-02-13 00:22:10 +00:00
Code Issues Projects Releases Wiki Activity

Fixed ssf for spherical equatorial coordinate system (old version was for polar)

Browse Source 
[SVN r72238]
This commit is contained in:
Barend Gehrels
2011-05-27 23:02:58 +00:00
parent 9da5cc5e3c
commit 48d5f65669
2 changed files with 100 additions and 28 deletions
Download Patch File Download Diff File Expand all files Collapse all files

58
include/boost/geometry/strategies/spherical/ssf.hpp
View File

@@ -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 <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

70
include/boost/geometry/strategies/strategy_transform.hpp
View File

@@ -155,29 +155,48 @@ namespace detail
/// 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 @@ struct from_spherical_polar_2_to_cartesian_3
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 @@ struct from_spherical_polar_3_to_cartesian_3
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 @@ struct default_strategy<spherical_polar_tag, cartesian_tag, CoordSys1, CoordSys2
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>