Conversion Framework to convert between Coordinate Systems

include/boost/astronomy/coordinate/conversion/conversion_graph.hpp

@@ -0,0 +1,125 @@
+/*=============================================================================
+Copyright 2020 Syed Ali Hasan <[email protected]>
+
+Distributed under the Boost Software License, Version 1.0. (See accompanying
+file License.txt or copy at https://www.boost.org/LICENSE_1_0.txt)
+=============================================================================*/
+
+#include <iostream>
+#include <stack>
+
+//Graph
+#include <boost/graph/graphviz.hpp>
+#include <boost/graph/properties.hpp>
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/breadth_first_search.hpp>
+#include <boost/graph/named_function_params.hpp>
+
+//Matrix
+#include <boost/numeric/ublas/io.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/astronomy/coordinate/utility/utility.hpp>
+
+//Angle
+#include <boost/units/io.hpp>
+#include <boost/units/systems/angle/degrees.hpp>
+#include <boost/units/systems/si/plane_angle.hpp>
+
+namespace bu = boost::units;
+namespace bg = boost::geometry;
+namespace bac = boost::astronomy::coordinate;
+
+using namespace boost::units;
+using namespace boost::units::si;
+using namespace boost::astronomy::coordinate;
+
+enum class COORDINATE_SYSTEM {
+  HORIZON,
+  EQUATORIAL_HA_DEC,
+  EQUATORIAL_RA_DEC,
+  ECLIPTIC,
+  GALACTIC};
+
+struct CoordinateData
+{
+  COORDINATE_SYSTEM coordinate_system;
+  std::string coordinate_name;
+};
+
+struct EdgeData
+{
+  std::string edge_label;
+  matrix<double> conv_matrix;
+};
+
+using graph_t = boost::adjacency_list<boost::vecS, boost::vecS,
+    boost::directedS,
+    CoordinateData,
+    EdgeData>;
+
+using vertex_t = boost::graph_traits<graph_t>::vertex_descriptor;
+
+template
+<
+  typename CoordinateType = double,
+  typename Angle = bu::quantity<bu::si::plane_angle, CoordinateType>
+>
+matrix<double> convert(const COORDINATE_SYSTEM src,
+                       const COORDINATE_SYSTEM dest,
+                       const Angle& phi,
+                       const Angle& st,
+                       const Angle& obliquity,
+                       coord_sys<2, bg::cs::spherical<bg::radian>, CoordinateType> source_coordinate)
+{
+  bac::column_vector<double,
+      quantity<bu::si::plane_angle, double>,
+      double>
+      col_vec(bg::get<0>(source_coordinate.get_point()) * radians, bg::get<1>(source_coordinate.get_point()) * radians);
+
+  graph_t G;
+
+  const int graph_size = 5;
+
+  vertex_t vd0 = boost::add_vertex(CoordinateData{COORDINATE_SYSTEM::HORIZON,"Horizon"}, G);
+  vertex_t vd1 = boost::add_vertex(CoordinateData{COORDINATE_SYSTEM::EQUATORIAL_HA_DEC,"Equatorial_HA_Dec"}, G);
+  vertex_t vd2 = boost::add_vertex(CoordinateData{COORDINATE_SYSTEM::EQUATORIAL_RA_DEC,"Equatorial_RA_Dec"}, G);
+  vertex_t vd3 = boost::add_vertex(CoordinateData{COORDINATE_SYSTEM::ECLIPTIC,"Ecliptic"}, G);
+  vertex_t vd4 = boost::add_vertex(CoordinateData{COORDINATE_SYSTEM::GALACTIC,"Galactic"}, G);
+
+  boost::add_edge(vd0, vd1, EdgeData{"Horizon to Equatorial HA Dec", bac::ha_dec_horizon<double, quantity<bu::degree::plane_angle>, double>(phi).get()}, G);
+  boost::add_edge(vd1, vd0, EdgeData{"Equatorial HA Dec to Horizon", bac::ha_dec_horizon<double, quantity<bu::degree::plane_angle>, double>(phi).get()}, G);
+  boost::add_edge(vd1, vd2, EdgeData{"Equatorial HA Dec to Equatorial RA Dec", bac::ha_dec_ra_dec<double, quantity<bu::degree::plane_angle>, double>(st).get()}, G);
+  boost::add_edge(vd2, vd1, EdgeData{"Equatorial RA Dec to Equatorial HA Dec", bac::ha_dec_ra_dec<double, quantity<bu::degree::plane_angle>, double>(st).get()}, G);
+  boost::add_edge(vd2, vd3, EdgeData{"Equatorial RA Dec to Ecliptic", bac::ra_dec_to_ecliptic<double, quantity<bu::degree::plane_angle>, double>(obliquity).get()}, G);
+  boost::add_edge(vd3, vd2, EdgeData{"Ecliptic to Equatorial RA Dec", bac::ecliptic_to_ra_dec<double, quantity<bu::degree::plane_angle>, double>(obliquity).get()}, G);
+  boost::add_edge(vd2, vd4, EdgeData{"Equatorial RA Dec to Galactic", bac::ra_dec_to_galactic<double>().get()}, G);
+  boost::add_edge(vd4, vd2, EdgeData{"Galactic to Equatorial RA Dec", bac::galactic_to_ra_dec<double>().get()}, G);
+
+  //Predecessor Array
+  boost::array<int, graph_size> predecessors{0};
+  predecessors[(int)src] = (int)src;
+
+  boost::breadth_first_search(G, (int)src, boost::visitor(
+      boost::make_bfs_visitor(
+          boost::record_predecessors(predecessors.begin(),
+                                     boost::on_tree_edge{}))));
+
+  //Get traversed path
+  std::vector<int> store_path;
+
+  int p = (int)dest;
+  while (p != (int)src)
+  {
+    store_path.push_back(p);
+    p = predecessors[p];
+  }
+  store_path.push_back(p);
+
+  matrix<double> ans = col_vec.get();
+
+  //Matrix Multiplication
+  for(auto it = store_path.rbegin(); it + 1 != store_path.rend(); ++it)
+    ans = prod(G[boost::edge(*it,*(it+1),G).first].conv_matrix, ans);
+
+  return ans;
+}

include/boost/astronomy/time/time_conversions.hpp

@@ -66,36 +66,49 @@ decimal_hour GST(ptime t)
 
 enum class DIRECTION {WEST, EAST};
 
+decimal_hour compute_LST(double longitude, DIRECTION direction, double gst)
+{
+  //Convert longitude to hours
+  double long_hours = longitude / 15.0;
+
+  switch(direction)
+  {
+    case DIRECTION::WEST:
+      //Multiply with direction
+      long_hours = -1 * long_hours;
+      break;
+    case DIRECTION::EAST:
+      //Multiply with direction
+      long_hours = 1 * long_hours;
+      break;
+  }
+
+  long_hours = long_hours + gst;
+
+  //Bring the result into the range 0 to 24 by adding or subtracting 24 if necessary.
+  //This is the local sidereal time (LST).
+  long_hours = long_hours - 24.0 * floor(long_hours/24.0);
+
+  return {long_hours};
+}
+
 //Local Sidereal Time (LST)
 decimal_hour LST(double longitude, DIRECTION direction, ptime t)
 {
   double gst = GST(t).get();
 
-    if(longitude == 0)
-      return {gst};
-
-    //Convert longitude to hours
-    double long_hours = longitude / 15.0;
+  if(longitude == 0)
+    return {gst};
 
-    switch(direction)
-    {
-      case DIRECTION::WEST:
-          //Multiply with direction
-          long_hours = -1 * long_hours;
-            break;
-      case DIRECTION::EAST:
-          //Multiply with direction
-          long_hours = 1 * long_hours;
-            break;
-    }
-
-    long_hours = long_hours + gst;
+  return compute_LST(longitude,direction,gst);
+}
 
-    //Bring the result into the range 0 to 24 by adding or subtracting 24 if necessary.
-    //This is the local sidereal time (LST).
-    long_hours = long_hours - 24.0 * floor(long_hours/24.0);
+decimal_hour LST(double longitude, DIRECTION direction, double gst)
+{
+  if(longitude == 0)
+    return {gst};
 
-    return {long_hours};
+  return compute_LST(longitude,direction,gst);
 }
 
 #endif //BOOST_ASTRONOMY_TIME_CONVERSIONS

test/CMakeLists.txt

@@ -7,3 +7,6 @@ add_subdirectory(units)
 add_subdirectory(time)
 
 add_subdirectory(io)
+
+add_subdirectory(conversion)
+

test/Jamfile

@@ -12,3 +12,7 @@ build-project header ;
 build-project coordinate ;
 build-project units ;
 build-project io ;
+build-project time ;
+build-project conversion ;
+
+

test/conversion/CMakeLists.txt

@@ -0,0 +1,16 @@
+foreach(_name
+        conversion_graph)
+    set(_target test_coordinate_${_name})
+
+    add_executable(${_target} "")
+    target_sources(${_target} PRIVATE ${_name}.cpp)
+    target_link_libraries(${_target}
+            PRIVATE
+            astronomy_compile_options
+            astronomy_include_directories
+            astronomy_dependencies)
+    add_test(NAME test.astro.${_name} COMMAND ${_target})
+
+    unset(_name)
+    unset(_target)
+endforeach()

test/conversion/Jamfile

@@ -0,0 +1,3 @@
+import testing ;
+
+run conversion_graph.cpp ;

test/conversion/conversion_graph.cpp

@@ -0,0 +1,72 @@
+/*=============================================================================
+Copyright 2020 Syed Ali Hasan <[email protected]>
+
+Distributed under the Boost Software License, Version 1.0. (See accompanying
+file License.txt or copy at https://www.boost.org/LICENSE_1_0.txt)
+=============================================================================*/
+
+#define BOOST_TEST_MODULE conversion_graph
+
+#include <iostream>
+
+#include <boost/units/io.hpp>
+#include <boost/units/systems/angle/degrees.hpp>
+#include <boost/units/systems/si/plane_angle.hpp>
+
+#include <boost/astronomy/coordinate/coord_sys/horizon_coord.hpp>
+#include <boost/astronomy/coordinate/coord_sys/ecliptic_coord.hpp>
+#include <boost/astronomy/coordinate/coord_sys/galactic_coord.hpp>
+#include <boost/astronomy/coordinate/coord_sys/equatorial_ra_coord.hpp>
+#include <boost/astronomy/coordinate/coord_sys/equatorial_ha_coord.hpp>
+
+#include <boost/astronomy/coordinate/conversion/conversion_graph.hpp>
+
+#include <boost/test/unit_test.hpp>
+
+using namespace boost::units;
+using namespace boost::units::si;
+using namespace boost::astronomy::coordinate;
+
+namespace bud = boost::units::degree;
+namespace bac = boost::astronomy::coordinate;
+
+BOOST_AUTO_TEST_SUITE(conversion_graph)
+
+BOOST_AUTO_TEST_CASE(ecliptic_to_horizon) {
+
+  ecliptic_coord<double, quantity<bud::plane_angle>, quantity<bud::plane_angle>>
+    ec(97.638119 * bud::degrees, -17.857969 * bud::degrees);
+
+  quantity<bud::plane_angle, double> phi = 52.175 * bud::degree;
+  quantity<bud::plane_angle, double> st = 77.337 * bud::degree;
+  quantity<bud::plane_angle, double> obliquity = 23.446 * bud::degree;
+
+  matrix<double> ans = convert(COORDINATE_SYSTEM::ECLIPTIC,COORDINATE_SYSTEM::HORIZON,phi,st,obliquity,ec);
+
+  auto theta = bac::extract_coordinates(ans).get_coordinates().first;
+  auto gama = bac::extract_coordinates(ans).get_coordinates().second;
+
+  BOOST_CHECK_CLOSE(theta.value() * 180.0 / PI, 153.491944, 0.001);
+  BOOST_CHECK_CLOSE(gama.value() * 180.0 / PI, 40.399444, 0.001);
+}
+
+BOOST_AUTO_TEST_CASE(horizon_to_ecliptic) {
+
+  horizon_coord<double, quantity<bud::plane_angle>, quantity<bud::plane_angle>>
+      hc(153.491944 * bud::degrees, 40.399444  * bud::degrees);
+
+  quantity<bud::plane_angle, double> phi = 52.175 * bud::degree;
+  quantity<bud::plane_angle, double> st = 77.337 * bud::degree;
+  quantity<bud::plane_angle, double> obliquity = 23.446 * bud::degree;
+
+  matrix<double> ans = convert(COORDINATE_SYSTEM::HORIZON,COORDINATE_SYSTEM::ECLIPTIC,phi,st,obliquity,hc);
+
+  auto theta = bac::extract_coordinates(ans).get_coordinates().first;
+  auto gama = bac::extract_coordinates(ans).get_coordinates().second;
+
+  BOOST_CHECK_CLOSE(theta.value() * 180.0 / PI, 97.638119 , 0.001);
+  BOOST_CHECK_CLOSE(gama.value() * 180.0 / PI, -17.857969, 0.001);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+