First commit

CMakeLists.txt

+cmake_minimum_required(VERSION 3.22)
+project(funwitherfa)
+
+
+set(SOURCES 
+    src/main.cpp
+)
+
+
+include(FetchContent)
+
+FetchContent_Declare(
+  aocommon
+  GIT_REPOSITORY https://gitlab.com/aroffringa/aocommon.git
+  GIT_TAG master)
+  
+FetchContent_MakeAvailable(aocommon)
+
+list(APPEND CMAKE_MODULE_PATH "${aocommon_SOURCE_DIR}/CMake")
+
+find_package(Casacore COMPONENTS measures)
+
+add_executable(funwitherfa ${SOURCES})
+target_include_directories(funwitherfa PRIVATE src ${CASACORE_INCLUDE_DIRS})
+
+target_link_libraries(funwitherfa ${CASACORE_LIBRARIES})
+
+# ERFA
+find_package(PkgConfig REQUIRED)
+pkg_check_modules(ERFA REQUIRED erfa)
+
+target_link_libraries(funwitherfa ${ERFA_LIBRARIES})
+target_include_directories(funwitherfa PRIVATE ${ERFA_INCLUDE_DIRS})
\ No newline at end of file

src/erfacpp.h

+#include <erfa.h>
+#include <erfaiers.h>
+#include <erfatime.h>
+
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+
+namespace erfacpp {
+
+double I3[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+
+typedef struct vec2d {
+  double phi;
+  double theta;
+  friend std::ostream& operator<<(std::ostream& os, const struct vec2d& obj) {
+    os << "(" << std::setprecision(9) << obj.phi << "," << obj.theta << ")";
+    return os;
+  }
+} direction;
+
+typedef struct vec3d {
+  double x;
+  double y;
+  double z;
+
+  friend std::ostream& operator<<(std::ostream& os, const struct vec3d& obj) {
+    os << "(" << std::setprecision(9) << obj.x << "," << obj.y << "," << obj.z
+       << ")";
+    return os;
+  }
+} position;
+
+void J2000_to_CIRS(direction j2000, UTCTime time, direction& cirs, double* eo) {
+  TTTime ttime = time.to<TTTime>();
+  direction icrs;
+  eraAtcc13(j2000.phi, j2000.theta, 0, 0, 0, 0, ttime.djm0, ttime.djm,
+            &icrs.phi, &icrs.theta);
+  eraAtci13(icrs.phi, icrs.theta, 0, 0, 0, 0, ttime.djm0, ttime.djm, &cirs.phi,
+            &cirs.theta, eo);
+}
+
+void J2000_to_HADEC(direction j2000, UTCTime time, position observer,
+                    direction& hadec) {
+  eraASTROM parameters;
+
+  double elong, phi, hm;
+  double eo;
+  eraGc2gd(1, &observer.x, &elong, &phi, &hm);
+  double deltaut1 = 0;
+  eraApco13(time.djm0, time.djm, 0, elong, phi, hm, 0, 0, 0, 0, 0, 0,
+            &parameters, &eo);
+
+  direction itrf;
+  eraAtccq(j2000.phi, j2000.theta, 0, 0, 0, 0, &parameters, &hadec.phi,
+           &hadec.theta);
+}
+
+void J2000_to_GCRS(direction j2000, UTCTime time, direction& gcrs) {
+  double eo = 0;
+  J2000_to_CIRS(j2000, time, gcrs, &eo);
+  gcrs.phi = eraAnp(gcrs.phi - eo);
+}
+
+void polar_motion(UTCTime time, double& xp, double& yp) {
+  // FOR NOW
+  xp = 1.696847883883376e-07;
+  yp = 1.4059596752176543e-06;
+}
+
+void J2000_to_ITRF(direction j2000, UTCTime time, position& itrf) {
+  direction cirs;
+  double eo;
+  position observer{0, 0, 0};
+  UT1Time ut1 = time.to<UT1Time>();
+  const double sp = time.to<TTTime>().ToTIO();
+  double xp, yp;
+  double rpom00[3][3];
+  double c2im[3][3];
+
+  polar_motion(time, xp, yp);
+
+  eraPom00(xp, yp, sp, rpom00);
+  const double era = eraEra00(ut1.djm0, ut1.djm);
+  eraC2tcio(I3, era, rpom00, c2im);
+  J2000_to_CIRS(j2000, time, cirs, &eo);
+  itrf.x = std::cos(cirs.phi) * std::cos(cirs.theta);
+  itrf.y = std::sin(cirs.phi) * std::cos(cirs.theta);
+  itrf.z = std::sin(cirs.theta);
+
+  itrf.x = c2im[0][0] * itrf.x + c2im[0][1] * itrf.y + c2im[0][2] * itrf.z;
+  itrf.y = c2im[1][0] * itrf.x + c2im[1][1] * itrf.y + c2im[1][2] * itrf.z;
+  itrf.z = c2im[2][0] * itrf.x + c2im[2][1] * itrf.y + c2im[2][2] * itrf.z;
+
+  //[ 1,  0, 0 ]     [ cos(lon)  sin(lon) 0]
+  //[ 0, -1, 0 ]  @  [-sin(lon)  cos(lon) 0]
+  //[ 0,  0, 1 ]     [0                0  1]
+}
+
+}  // namespace erfacpp

src/erfaiers.h

+#ifndef ERFACPP_IERS
+#define ERFACPP_IERS
+
+#include <erfatime.h>
+
+#include <fstream>
+#include <iostream>
+#include <limits>
+#include <string>
+#include <vector>
+
+namespace erfacpp {
+
+class TableParser {
+ public:
+  TableParser(const std::string& path) : path_{path} {};
+
+  inline const int GetColumnIndex(const std::string& name) {
+    return col_names_[name];
+  }
+
+  void parse() {
+    std::ifstream f_stream_in = std::ifstream(path_);
+
+    const std::vector<int> skip_map = SkipMap();
+    const size_t n_columns = cols_sizes_.size();
+
+    for (std::string line; std::getline(f_stream_in, line);) {
+      const int current_line = values.size();
+      std::vector<double> current_values;
+      for (int col_idx = 0; col_idx < n_columns; col_idx++) {
+        if (use_cols_[col_idx] && line.size()) {
+          const std::string value =
+              line.substr(skip_map[col_idx], cols_sizes_[col_idx]);
+          try {
+            current_values.push_back(std::stod(value));
+          } catch (std::invalid_argument e) {
+            std::cout << "Invalid argument: '" << value << "'" << std::endl;
+            use_cols_[col_idx] = false;
+          }
+        }
+      }
+      values.push_back(current_values);
+    }
+    std::cout << "Parsed" << std::endl;
+  };
+
+ protected:
+  std::vector<std::vector<double>> values;
+  void AddColumn(const short size, const std::string& name, const bool store) {
+    if (store) col_names_[name] = cols_sizes_.size();
+    cols_sizes_.push_back(size);
+    use_cols_.push_back(store);
+  }
+
+ private:
+  std::string path_;
+  std::vector<int> cols_sizes_;
+  std::vector<bool> use_cols_;
+  std::map<std::string, int> col_names_;
+
+  const std::vector<int> SkipMap() {
+    std::vector<int> skip_map(1);
+
+    skip_map[0] = 0;
+    for (int i = 0; i < cols_sizes_.size(); i++) {
+      skip_map.push_back(cols_sizes_[i] + skip_map[i]);
+    }
+    return skip_map;
+  };
+};
+
+class PMTable {};
+
+class IERSTable : public TableParser {
+ public:
+  IERSTable(const std::string path) : TableParser(path) {
+    AddColumn(2, "year", false);
+    AddColumn(2, "month", false);
+    AddColumn(2, "day", false);
+    AddColumn(1, "space", false);
+    AddColumn(8, "mdj", false);
+    AddColumn(1, "space", false);
+    AddColumn(1, "Flag_PM", false);
+    AddColumn(1, "space", false);
+    AddColumn(9, "PM_x", true);
+    AddColumn(9, "PM_x_err", true);
+    AddColumn(1, "space", false);
+    AddColumn(9, "PM_y", true);
+    AddColumn(9, "PM_y_err", true);
+    AddColumn(2, "space", false);
+    AddColumn(1, "Flag_U1UTC", false);
+    AddColumn(10, "U1UTC", true);
+    AddColumn(10, "U1UTC_err", true);
+
+    parse();
+  };
+
+  const double GetValue(const std::string& value, UTCTime time) {
+    const int value_idx = GetColumnIndex(value);
+    const int time_idx = GetColumnIndex(value);
+
+    for (int i = values.size(); i > 0; i--) {
+      if (values[i - 1].size() == 0) {
+        continue;
+      } else {
+        std::cout << values[i - 1][time_idx] << std::endl;
+        if (time.djm < values[i - 1][time_idx]) {
+          return values[i - 1][value_idx];
+        }
+      }
+    }
+
+    return std::numeric_limits<double>::quiet_NaN();
+  }
+};
+
+class LODTable {};
+
+class ADVTable {};
+
+class APMTable {};
+}  // namespace erfacpp
+#endif  // ERFACPP_IERS
\ No newline at end of file

src/erfatime.h

+#ifndef ERFACPP_TIME
+#define ERFACPP_TIME
+#include <erfa.h>
+
+#include <iostream>
+
+namespace erfacpp {
+template <class TBase>
+struct OTime {
+  double djm0;
+  double djm;
+  template <typename T>
+  T to() {
+    throw std::runtime_error("conversion not implemented");
+  }
+
+  friend std::ostream& operator<<(std::ostream& os, const OTime& obj) {
+    os << "(" << obj.djm0 << ",\t" << obj.djm << ")";
+    return os;
+  }
+};
+
+struct UTCTime : OTime<UTCTime> {
+ public:
+  static UTCTime FromMJD(double mjd) {
+    UTCTime res;
+    res.djm0 = 2400000.5;
+    res.djm = mjd;
+    return res;
+  }
+};
+
+struct TTTime : OTime<TTTime> {
+  const double ToTIO() { return eraSp00(djm0, djm); };
+};
+struct TAITime : OTime<TAITime> {};
+struct UT1Time : OTime<UT1Time> {};
+
+template <>
+template <>
+TAITime OTime<UTCTime>::to<TAITime>() {
+  TAITime res;
+  eraUtctai(djm0, djm, &res.djm0, &res.djm);
+  return res;
+}
+
+template <>
+template <>
+TTTime OTime<TAITime>::to<TTTime>() {
+  TTTime res;
+  eraTaitt(djm0, djm, &res.djm0, &res.djm);
+  return res;
+}
+
+template <>
+template <>
+UT1Time OTime<TTTime>::to<UT1Time>() {
+  UT1Time res;
+  eraTtut1(djm0, djm, 0, &res.djm0, &res.djm);
+  return res;
+}
+
+template <>
+template <>
+TTTime OTime<UTCTime>::to<TTTime>() {
+  return this->to<TAITime>().to<TTTime>();
+}
+
+template <>
+template <>
+UT1Time OTime<UTCTime>::to<UT1Time>() {
+  return this->to<TAITime>().to<TTTime>().to<UT1Time>();
+}
+
+}  // namespace erfacpp
+
+#endif  // ERFACPP_TIME
\ No newline at end of file

src/main.cpp

+
+#include <casacore/measures/Measures/MCDirection.h>
+#include <casacore/measures/Measures/MDirection.h>
+#include <casacore/measures/Measures/MEpoch.h>
+#include <casacore/measures/Measures/MPosition.h>
+#include <casacore/measures/Measures/MeasConvert.h>
+#include <casacore/measures/Measures/MeasFrame.h>
+#include <casacore/measures/Measures/SolarPos.h>
+#include <erfacpp.h>
+
+#include <iostream>
+#include <string>
+using casacore::MDirection;
+using casacore::MeasFrame;
+using casacore::MVPosition;
+using casacore::SolarPos;
+using erfacpp::direction;
+using erfacpp::position;
+
+const direction test_dir{1.234, 1.57};
+const position test_position{826577.022720000, 461022.995082000, 5064892.814};
+const double mjd = 60395.47027855601;
+const casacore::MEpoch epoch(casacore::Quantity(mjd, casacore::Unit('d')));
+
+void j2000_to_jnat_casacore(direction d, position p, direction& d_new) {
+  casacore::MVDirection in = casacore::MVDirection(d.phi, d.theta);
+
+  casacore::MVPosition mv_position(p.x, p.y, p.z);
+  casacore::MPosition m_position(mv_position, casacore::MPosition::ITRF);
+  MeasFrame frame = casacore::MeasFrame(epoch, m_position);
+
+  MDirection dir(in, MDirection::J2000);
+  casacore::MDirection::Convert converter = casacore::MDirection::Convert(
+      dir, casacore::MDirection::Ref(casacore::MDirection::JNAT, frame));
+  const casacore::MDirection jnat = converter();
+  const casacore::MVDirection converted = jnat.getValue();
+  d_new.phi = converted.getLong();
+  d_new.theta = converted.getLat();
+}
+
+void j2000_to_jnat(direction d, position p, direction& d_new) {
+  SolarPos pos(SolarPos::STANDARD);
+  double epoch = casacore::MEpoch().getData()->getVector()[0];
+  std::cout << "Epoch at 0 " << casacore::MEpoch() << std::endl;
+  std::cout << "Epoch at 0 " << epoch << std::endl;
+  MVPosition solar_pos = pos(0);
+}
+
+void erfa_j2000_to_cirs(direction d, position p, direction& d_new) {
+  double eo;
+  erfacpp::J2000_to_CIRS(d, erfacpp::UTCTime::FromMJD(mjd), d_new, &eo);
+}
+
+void erfa_j2000_to_gcrs(direction d, position p, direction& d_new) {
+  erfacpp::J2000_to_GCRS(d, erfacpp::UTCTime::FromMJD(mjd), d_new);
+}
+
+void j2000_to_ITRF_casacore(direction d, position p, direction& d_new) {
+  casacore::MVDirection in = casacore::MVDirection(d.phi, d.theta);
+
+  casacore::MVPosition mv_position(p.x, p.y, p.z);
+  casacore::MPosition m_position(mv_position, casacore::MPosition::ITRF);
+  MeasFrame frame = casacore::MeasFrame(epoch, m_position);
+
+  MDirection dir(in, MDirection::J2000);
+  casacore::MDirection::Convert converter = casacore::MDirection::Convert(
+      dir, casacore::MDirection::Ref(casacore::MDirection::ITRF, frame));
+  const casacore::MDirection itrf = converter();
+  const casacore::MVDirection converted = itrf.getValue();
+  d_new.phi = converted.getLong();
+  d_new.theta = converted.getLat();
+  std::cout << converted.getVector() << std::endl;
+}
+
+void j2000_to_AZEL_casacore(direction d, position p, direction& d_new) {
+  casacore::MVDirection in = casacore::MVDirection(d.phi, d.theta);
+
+  casacore::MVPosition mv_position(p.x, p.y, p.z);
+  casacore::MPosition m_position(mv_position, casacore::MPosition::ITRF);
+  MeasFrame frame = casacore::MeasFrame(epoch, m_position);
+
+  MDirection dir(in, MDirection::J2000);
+  casacore::MDirection::Convert converter = casacore::MDirection::Convert(
+      dir, casacore::MDirection::Ref(casacore::MDirection::AZEL, frame));
+  const casacore::MDirection azel = converter();
+  const casacore::MVDirection converted = azel.getValue();
+  d_new.phi = converted.getLong();
+  d_new.theta = converted.getLat();
+}
+
+void j2000_to_HADEC_casacore(direction d, position p, direction& d_new) {
+  casacore::MVDirection in = casacore::MVDirection(d.phi, d.theta);
+
+  casacore::MVPosition mv_position(p.x, p.y, p.z);
+  casacore::MPosition m_position(mv_position, casacore::MPosition::ITRF);
+  MeasFrame frame = casacore::MeasFrame(epoch, m_position);
+
+  MDirection dir(in, MDirection::J2000);
+  casacore::MDirection::Convert converter = casacore::MDirection::Convert(
+      dir, casacore::MDirection::Ref(casacore::MDirection::HADEC, frame));
+  const casacore::MDirection hadec = converter();
+  const casacore::MVDirection converted = hadec.getValue();
+  d_new.phi = converted.getLong();
+  d_new.theta = converted.getLat();
+}
+
+void erfa_j2000_to_hadec(direction d, position p, direction& d_new) {
+  erfacpp::J2000_to_HADEC(d, erfacpp::UTCTime::FromMJD(mjd), p, d_new);
+}
+void erfa_j2000_to_itrf(direction d, position& p_new) {
+  erfacpp::J2000_to_ITRF(d, erfacpp::UTCTime::FromMJD(mjd), p_new);
+}
+
+int main(int argc, char* argv[]) {
+  std::cout << "Welcome at fun with ERFA.\n Your coordinate transformation "
+               "best library.\n"
+            << std::endl;
+
+  erfacpp::IERSTable table(
+      std::string("/home/mancini/git/funwitherfa/src/tables/finals2000A.all"));
+
+  std::cout << table.GetValue("mjd", erfacpp::UTCTime::FromMJD(41703))
+            << std::endl;
+
+  direction new_dir;
+  direction dir_cirs{2.8503633615820227, 1.5686321729185995};
+  direction dir_itrf{0.007436305256449476, 1.3858688716237009};
+  position pos_itrf{-0.00213085, 0.00037689, 0.99999766};
+  position new_pos;
+
+  std::cout << "J2000_to_JNAT_casacore" << std::endl;
+  j2000_to_jnat_casacore(test_dir, test_position, new_dir);
+  std::cout << "Direction before: " << test_dir << std::endl;
+  std::cout << "Direction after: " << new_dir << std::endl;
+
+  std::cout << "J2000_to_CIRS_erfa" << std::endl;
+
+  erfa_j2000_to_cirs(test_dir, test_position, new_dir);
+  std::cout << "Direction before: " << test_dir << std::endl;
+  std::cout << "Direction after: " << new_dir << std::endl;
+  std::cout << "Direction expected: " << dir_cirs << std::endl;
+
+  std::cout << "J2000_to_ITRF_casacore" << std::endl;
+  j2000_to_ITRF_casacore(test_dir, test_position, new_dir);
+  std::cout << "Direction before: " << test_dir << std::endl;
+  std::cout << "Direction after: " << new_dir << std::endl;
+
+  std::cout << "J2000_to_GCRS_erfa" << std::endl;
+  erfa_j2000_to_gcrs(test_dir, test_position, new_dir);
+  std::cout << "Direction before: " << test_dir << std::endl;
+  std::cout << "Direction after: " << new_dir << std::endl;
+
+  std::cout << "J2000_to_AZEL_casacore" << std::endl;
+  j2000_to_AZEL_casacore(test_dir, test_position, new_dir);
+  std::cout << "Direction before: " << test_dir << std::endl;
+  std::cout << "Direction after: " << new_dir << std::endl;
+
+  std::cout << "J2000_to_ITRF_erfa" << std::endl;
+  erfa_j2000_to_itrf(test_dir, new_pos);
+  std::cout << "Direction before: " << test_dir << std::endl;
+  std::cout << "Direction after: " << new_pos << std::endl;
+  std::cout << "Direction after: " << pos_itrf << std::endl;
+
+  std::cout << "J2000_to_HADEC_casacore" << std::endl;
+  j2000_to_HADEC_casacore(test_dir, test_position, new_dir);
+  std::cout << "Direction before: " << test_dir << std::endl;
+  std::cout << "Direction after: " << new_dir << std::endl;
+
+  std::cout << "J2000_to_hadec_erfa" << std::endl;
+  erfa_j2000_to_hadec(test_dir, test_position, new_dir);
+  std::cout << "Direction before: " << test_dir << std::endl;
+  std::cout << "Direction after: " << new_dir << std::endl;
+  return 0;
+}
\ No newline at end of file