minor

xanthospap · xanthospap · commit d76102b4ec41 · 2025-12-04T15:11:43.000+02:00
diff --git a/include/dpod.hpp b/include/dpod.hpp
@@ -18,10 +18,110 @@
 #include <cstring>
 
 namespace dso {
+/** @brief Parse harmonics off from a dpod_freq_cor file and compute 
+ * (cartesian) corrections (ΔX, ΔY, ΔZ).
+ *
+ * Give a dpod20*_freq_corr.txt file, do the following things:
+ * 1. search through all of the sites included in sites_crd and get their 
+ *    frequency/harmonics info,
+ * 2. compute the total, accumulated harmonics contribution (i.e. for all of 
+ *    the included frequencies) in cartesian components (ΔX, ΔY, ΔZ)
+ * 3. store results of step 2 in an std::vector<Sinex::SiteCoordinateResults> 
+ *    and return it.
+ * 
+ * Hence, if the function is successeful, to get the 'correct' position of a 
+ * given site at the sites_crd vector, the contribution returned from the call 
+ * should be added.
+ *
+ * Note that the resulting vector (std::vector<Sinex::SiteCoordinateResults>) 
+ * will have a one-to-one correspondance with the input sites_crd vector, i.e. 
+ * the sites will be placed exactly in the same order.
+ *
+ * Why do we need a std::vector<Sinex::SiteCoordinateResults> as input? Well, 
+ * this vector will (probably) hold results/positions from a respective 
+ * dpod SINEX file. But, we must know the SOLN_ID, i.e. the solution id that 
+ * was used to compute/extrapolate the coordinates. The same solution id must 
+ * be used to identify the 'correct' harmonics to be applyied.
+ *
+ * See for example the following block: 
+ * ---------------------------------------------------------------------------
+#CODE PT __DOMES__SOLN_XYZ_COSAMP__COSSTD__SINAMP__SINSTD
+ REUB  A 97401S002  2   X   2.564   0.183   0.060   0.183
+ REUB  A 97401S002  2   Y   2.498   0.151   1.706   0.151
+ REUB  A 97401S002  2   Z  -2.041   0.099   0.593   0.099
+ ROTA  A 66007S001  2   X   0.287   0.106  -0.930   0.106
+ ROTA  A 66007S001  2   Y   1.241   0.086   1.994   0.086
+ ROTA  A 66007S001  2   Z  -1.842   0.084  -2.067   0.084
+ AJAB  A 10077S002  1   X  -0.000   0.001   0.000   0.001
+ AJAB  A 10077S002  1   Y  -0.000   0.001   0.000   0.001
+ AJAB  A 10077S002  1   Z   0.000   0.001  -0.000   0.001
+ TRIB  A 30604S002  1   X   0.884   0.161  -0.812   0.163
+ TRIB  A 30604S002  1   Y   0.409   0.255  -0.705   0.259
+ TRIB  A 30604S002  1   Z  -1.148   0.139  -0.276   0.141
+ SAKA  A 12329S001  3   X   4.723   0.294   0.807   0.293
+ SAKA  A 12329S001  3   Y  -2.538   0.345  -2.942   0.345
+ * ---------------------------------------------------------------------------
+ * has harmonics insformation w.r.t different solution id's; hence, this file 
+ * cannot be treated without prior knowledge of the information in the SINEX 
+ * file.
+ *
+ * @param[in] fn Filename of the dpod20*_freq_corr.txt file
+ * @param[in] t Epoch at which to compute the harmonic contributions. If f is 
+ *              the frequency of the harmonic, then the harmonic correction is 
+ *              computed as: A_{cos} * cos(2*π*fdoy / f)
+ *                          +A_{sin} * sin(2*π*fdoy / f)
+ * @param[in] sites_crd An std::vector<Sinex::SiteCoordinateResults> where the 
+ *            requested sites and solution id's are retrieved from. We are
+ *            going to compute harmonic corrections for each of the sites 
+ *            included in this vector, tagged with the same solution id as the 
+ *            one recorded in here. Wea re matching against (i.e. identifying 
+ *            same sites) by: 
+ *            - site_code
+ *            - point_code
+ *            - domes, and
+ *            - solution_id
+ * @return A std::vector<Sinex::SiteCoordinateResults>, in exactly the same 
+ *         order as the input sites_crd vector, but in its (x,y,z) coordinates 
+ *         the accumulated frequency corrections are stored (in [m]). In case 
+ *         a site is not included in the file or has no harmonic terms, its 
+ *         (x,y,z) coordinates will be (0,0,0).
+ */
+[[nodiscard]]
+std::vector<Sinex::SiteCoordinateResults> get_dpod_freq_corr(
+    const char *fn, const datetime<dso::nanoseconds> &t,
+    const std::vector<Sinex::SiteCoordinateResults> &sites_crd);
+
+/** @brief Parse harmonics off from a dpod_freq_cor file and compute 
+ * (cartesian) corrections (ΔX, ΔY, ΔZ).
+ *
+ * This function does exaclty the same as get_dpod_freq_corr, but instead of 
+ * returning the accumulated frequency correcions, it adds them the passed-in 
+ * sites_crd vector (inplace).
+ *
+ * @param[in] fn Filename of the dpod20*_freq_corr.txt file
+ * @param[in] t Epoch at which to compute the harmonic contributions. If f is 
+ *              the frequency of the harmonic, then the harmonic correction is 
+ *              computed as: A_{cos} * cos(2*π*fdoy / f)
+ *                          +A_{sin} * sin(2*π*fdoy / f)
+ * @param[in] sites_crd An std::vector<Sinex::SiteCoordinateResults> where the 
+ *            requested sites and solution id's are retrieved from. We are
+ *            going to compute harmonic corrections for each of the sites 
+ *            included in this vector, tagged with the same solution id as the 
+ *            one recorded in here. We are matching against (i.e. identifying 
+ *            same sites) by: 
+ *            - site_code
+ *            - point_code
+ *            - domes, and
+ *            - solution_id
+ *            At output, the vector will hold the combined input+freq_cor 
+ *            terms. I.e., the computed frequency corrections terms will be
+ *            added (per site) each of the vector's input holdings.
+ * @return An integer, where anything other than 0 denotes an error.
+ */
+[[nodiscard]]
 int apply_dpod_freq_corr(
-    const char *fn, 
-    const dso::datetime<dso::nanoseconds> &t,
-    const std::vector<Sinex::SiteCoordinateResults> &sites_crd) noexcept;
+    const char *fn, const datetime<dso::nanoseconds> &t,
+    std::vector<Sinex::SiteCoordinateResults> &sites_crd) noexcept;
 } /* namespace dso */
 
 #endif
diff --git a/src/dpod_extrapolate.cpp b/src/dpod_extrapolate.cpp
@@ -1,6 +1,9 @@
 #include "utilities.hpp"
 
-int dso::dpod_extrapolate(const dso::datetime<dso::nanoseconds> &t, const std::vector<const char *> &sites_4charid, const char *dpod_snx, const char *dpod_freq=nullptr) noexcept {
+int dso::dpod_extrapolate(const dso::datetime<dso::nanoseconds> &t,
+                          const std::vector<const char *> &sites_4charid,
+                          const char *dpod_snx,
+                          const char *dpod_freq = nullptr) noexcept {
   /* create the sinex instance */
   dso::Sinex snx(dpod_snx);
 
@@ -24,10 +27,13 @@ int dso::dpod_extrapolate(const dso::datetime<dso::nanoseconds> &t, const std::v
 
   /* append harmonics signal(s) -> crd */
   if (dpod_freq) {
-  if (dso::apply_dpod_freq_corr(dpod_freq, t, crd)) {
-    fprintf(stderr, "[ERROR] Failed applying dpod frequency corrections; file is %s (traceback: %s)\n", dpod_freq, __func__);
-    return 1;
-  }
+    if (dso::apply_dpod_freq_corr(dpod_freq, t, crd)) {
+      fprintf(stderr,
+              "[ERROR] Failed applying dpod frequency corrections; file is %s "
+              "(traceback: %s)\n",
+              dpod_freq, __func__);
+      return 1;
+    }
   }
 
   return 0;
diff --git a/src/parse_dpod_freq_corr.cpp b/src/parse_dpod_freq_corr.cpp
@@ -3,7 +3,7 @@
 #include <charconv>
 #include <cstdio>
 #include <fstream>
-#include <limits>
+#include <stdexcept>
 
 namespace {
 
@@ -89,11 +89,15 @@ constexpr const int sdata_start = 27;
  * 012345678901234567890123456789012345678901234567890123456789
  *           10        20        30        40        50
  */
-int resolve_freq_cor_data_line(const char *line, char cmp, double *data) noexcept {
+int resolve_freq_cor_data_line(const char *line, char cmp,
+                               double *data) noexcept {
   /* get component (one char) */
   cmp = line[24];
   if ((cmp != 'X' && cmp != 'Y') && (cmp != 'Z')) {
-    fprintf(stderr, "[ERROR] Failed resolving dpod/freq line [%s], invalid XYZ field! (traceback: %s)\n", line, __func__);
+    fprintf(stderr,
+            "[ERROR] Failed resolving dpod/freq line [%s], invalid XYZ field! "
+            "(traceback: %s)\n",
+            line, __func__);
     return 1;
   }
   const char *str = line + 27;
@@ -108,22 +112,25 @@ int resolve_freq_cor_data_line(const char *line, char cmp, double *data) noexcep
 }
 } /* anonymous namespace */
 
-int dso::apply_dpod_freq_corr(
-    const char *fn, 
-    const dso::datetime<dso::nanoseconds> &t,
-    const std::vector<dso::Sinex::SiteCoordinateResults> &sites_crd) noexcept {
+std::vector<dso::Sinex::SiteCoordinateResults> dso::get_dpod_freq_corr(
+    const char *fn, const dso::datetime<dso::nanoseconds> &t,
+    const std::vector<dso::Sinex::SiteCoordinateResults> &sites_crd) {
   /* open dpod freq file */
   std::ifstream fin(fn);
   if (!fin.is_open()) {
     fprintf(stderr,
             "[ERROR] Failed opening dpod freq_corr file %s (traceback: %s)\n",
             fn, __func__);
-    return 1;
+    throw std::runtime_error("[ERROR] Failed opening dpod freq_corr file\n");
   }
 
   /* copy of input coordinates */
   std::vector<dso::Sinex::SiteCoordinateResults> scpy(sites_crd);
 
+  /* zero-out (cartesian) positions */
+  for (auto it = scpy.begin(); it != scpy.end(); ++it)
+    it->x = it->y = it->z = 0e0;
+
   /* fractional day of year and phase at epoch */
   const auto ymd_ = t.as_ydoy();
   const int idoy_ = ymd_.dy().as_underlying_type();
@@ -134,41 +141,46 @@ int dso::apply_dpod_freq_corr(
   char line[SZ];
   double cfreq = 0;
   double data[4];
-  char ccmp;
+  char ccmp = 'Q';
   double omega = 0e0;
-  int nfreq=0;
+  int nfreq = 0;
 
   int error = 0;
   while (fin.getline(line, SZ) && (!error)) {
     /* skip lines starting with '#', else parse */
     if (line[0] != '#') {
       /* first check site name */
       auto it = std::find_if(
-          scpy.cbegin(), scpy.cend(), [&](const dso::Sinex::SiteCoordinateResults &site) {
+          scpy.begin(), scpy.end(),
+          [&](const dso::Sinex::SiteCoordinateResults &site) {
             return ((!std::strncmp(site.msite.site_code(), line + 1, 4) &&
-                   !std::strncmp(site.msite.point_code(), line + 6, 2)) &&
-                   (!std::strncmp(site.msite.domes(), line + 9, 9) &&
-                   !std::strncmp(site.msolnid, line + 18, 4)));
+                     !std::strncmp(site.msite.point_code(), line + 6, 2)) &&
+                    (!std::strncmp(site.msite.domes(), line + 9, 9) &&
+                     !std::strncmp(site.msolnid, line + 18, 4)));
           });
       /* the station is in the list */
       if (it != scpy.cend()) {
 
         error += resolve_freq_cor_data_line(line, ccmp, data);
         if (!error) {
-          const double valmm = data[0] * std::cos(omega) + data[2] * std::sin(omega);
+          const double valmm =
+              data[0] * std::cos(omega) + data[2] * std::sin(omega);
           switch (ccmp) {
-            case 'X':
-              it->x += valmm*1e-3;
-              break;
-            case 'Y':
-              it->y += valmm*1e-3;
-              break;
-            case 'Z':
-              it->z += valmm*1e-3;
-              break;
-            default:
-              fprintf(stderr, "[ERROR] Invalid component in dpod_freq file %s; expected one of X, Y or Z; line is %s (traceback: %s)\n", fn, line);
-              error += 1;
+          case 'X':
+            it->x += valmm * 1e-3;
+            break;
+          case 'Y':
+            it->y += valmm * 1e-3;
+            break;
+          case 'Z':
+            it->z += valmm * 1e-3;
+            break;
+          default:
+            fprintf(stderr,
+                    "[ERROR] Invalid component in dpod_freq file %s; expected "
+                    "one of X, Y or Z; line is %s (traceback: %s)\n",
+                    fn, line, __func__);
+            error += 1;
           }
         } /* no coefficients parsing error */
       } /* if station in list */
@@ -181,5 +193,45 @@ int dso::apply_dpod_freq_corr(
   } /* reading through file entries */
 
   /* we should have read at least one list of frequencies */
-  return !((error==0) && (nfreq>=1));
+  if (!((error == 0) && (nfreq >= 1))) {
+    fprintf(stderr,
+            "[ERROR] Failed applying parsing/applying frequency corrections "
+            "from dpod file %s (traceback: %s)\n",
+            fn, __func__);
+    throw std::runtime_error(
+        "[ERROR] Failed applying parsing/applying frequency corrections\n");
+  }
+  return scpy;
+}
+
+int dso::apply_dpod_freq_corr(
+    const char *fn, const dso::datetime<dso::nanoseconds> &t,
+    std::vector<dso::Sinex::SiteCoordinateResults> &sites_crd) noexcept {
+  int error = 0;
+  try {
+    const auto cor = dso::get_dpod_freq_corr(fn, t, sites_crd);
+    int idx = 0;
+    for (auto it = sites_crd.begin(); it != sites_crd.end() && (!error); ++it) {
+      /* we assume here 1-to-1 correspondance between sites_crd and cor */
+      auto itc = cor.cbegin() + idx;
+      if ((std::strcmp(itc->msite.site_code(), it->msite.site_code()) ||
+           std::strcmp(itc->msite.point_code(), it->msite.point_code())) ||
+          (std::strcmp(itc->msite.domes(), it->msite.domes()) ||
+           std::strcmp(itc->msolnid, it->msolnid))) {
+        error += 1;
+        fprintf(stderr, "[ERROR] Corrupt(?) parsing of dpod_freq_cor file; expected correction for %s and found %s (traceback: %s)\n", itc->msite.site_code(), it->msite.site_code(), __func__);
+      }
+      it->x += itc->x;
+      it->y += itc->y;
+      it->z += itc->z;
+    }
+  } catch (std::exception &e) {
+    error = 1;
+    fprintf(stderr,
+            "[ERROR] Failed applying dpod frequency corrections from file %s "
+            "(traceback: %s)\n",
+            fn, __func__);
+    fprintf(stderr, "[ERROR] Caught exception %s", e.what());
+  }
+  return error;
 }
diff --git a/test/test_get_dpod_harmonics.cpp b/test/test_get_dpod_harmonics.cpp
@@ -9,36 +9,5 @@ int main(int argc, char* argv[])
     return 2;
   }
 
-  /* place sites in a vector */
-  std::vector<const char*> sites;
-  for (int i = 3; i < argc; i++) {
-    sites.push_back(argv[i]);
-  }
-
-  /* create the sinex instance */
-  dso::Sinex snx(argv[1]);
-
-  /* a vector of SiteId's to hold (intermediate) results */
-  std::vector<dso::sinex::SiteId> siteids;
-
-  /* parse the block SITE/ID to collect info for the given sites */
-  if (snx.parse_block_site_id(sites, /*use domes=*/false, siteids)) {
-    fprintf(stderr, "ERROR. Failed matching sites in SINEX file\n");
-    return 1;
-  }
-
-  std::vector<dso::SiteRealHarmonics> harm;
-
-  if (dso::parse_dpod_freq_corr(argv[2], siteids, harm)) {
-    fprintf(stderr, "ERROR. Failed parsing harmonics from %s\n", argv[2]);
-    return 3;
-  }
-
-  for (const auto& h : harm) {
-    for (int i = 0; i < h.harmonics().num_harmonics(); i++) {
-      printf("%s %.2f %.3f %.3f\n", h.site_name(), *(h.harmonics().operator()(i) + 0), *(h.harmonics().operator()(i) + 2), *(h.harmonics().operator()(i) + 1));
-    }
-  }
-
   return 0;
 }
