.. _program_listing_file_interfaces_spice.cc:

Program Listing for File spice.cc
=================================

|exhale_lsh| :ref:`Return to documentation for file <file_interfaces_spice.cc>` (``interfaces/spice.cc``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   
   #include "lupnt/interfaces/spice.h"
   
   #include <cspice/SpiceUsr.h>
   #include <cspice/SpiceZfc.h>
   #include <string.h>
   
   #include <cstdlib>
   #include <ctime>
   #include <filesystem>
   #include <fstream>
   #include <iostream>
   #include <mutex>
   #include <optional>
   #include <regex>
   #include <sstream>
   #include <stdexcept>
   
   #include "lupnt/conversions/coordinate_conversions.h"
   #include "lupnt/conversions/time_conversions.h"
   #include "lupnt/core/error.h"
   #include "lupnt/core/file.h"
   #include "lupnt/core/logger.h"
   #include "lupnt/environment/body.h"
   #include "lupnt/interfaces/spice_cheby.h"
   #include "lupnt/numerics/math_utils.h"
   
   namespace lupnt {
   
     namespace spice {
       static segment_t* cheby_s;
       static long cheby_n;
   
       bool spice_loaded = false;
   
       void CheckSpiceFailure(const std::string& context) {
         if (!failed_c()) return;
   
         SpiceChar short_msg[1841];
         SpiceChar long_msg[1841];
         getmsg_c("SHORT", sizeof(short_msg), short_msg);
         getmsg_c("LONG", sizeof(long_msg), long_msg);
         reset_c();
   
         throw std::runtime_error(fmt::format("{} failed: {} {}", context, short_msg, long_msg));
       }
   
       void SetSpiceReturnMode() {
         SpiceChar action[] = "RETURN";
         erract_c("SET", 0, action);
       }
   
       namespace {
         constexpr const char* kNaifGenericKernelsUrl
             = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels";
   
         bool IsKernelDownloadDisabled() {
           const char* val = std::getenv("LUPNT_SKIP_SPICE_KERNEL_DOWNLOAD");
           if (val == nullptr) return false;
           std::string s(val);
           return !(s.empty() || s == "0" || s == "false" || s == "FALSE" || s == "OFF" || s == "off");
         }
   
         bool RunShellCommand(const std::string& cmd) { return std::system(cmd.c_str()) == 0; }
   
         bool DownloadToFile(const std::string& url, const std::filesystem::path& dest_path) {
           if (IsKernelDownloadDisabled()) return false;
   
           std::filesystem::path tmp_path = dest_path;
           tmp_path += ".part";
           std::error_code ec;
           std::filesystem::remove(tmp_path, ec);
   
           std::string cmd
               = fmt::format("curl -fsSL --connect-timeout 10 --max-time 600 -o \"{}\" \"{}\"",
                             tmp_path.string(), url);
           bool ok = RunShellCommand(cmd) && std::filesystem::exists(tmp_path)
                     && std::filesystem::file_size(tmp_path) > 0;
           if (!ok) {
             std::filesystem::remove(tmp_path, ec);
             return false;
           }
           std::filesystem::rename(tmp_path, dest_path, ec);
           if (ec) {
             // `rename` can fail across filesystems; fall back to copy + remove.
             ec.clear();
             std::filesystem::copy_file(tmp_path, dest_path,
                                        std::filesystem::copy_options::overwrite_existing, ec);
             std::filesystem::remove(tmp_path, ec);
           }
           return !ec;
         }
   
         std::optional<std::string> FetchTextFromUrl(const std::string& url) {
           if (IsKernelDownloadDisabled()) return std::nullopt;
   
           auto tmp_path
               = std::filesystem::temp_directory_path()
                 / fmt::format("lupnt_naif_listing_{}.html", static_cast<long>(std::time(nullptr)));
           std::string cmd
               = fmt::format("curl -fsSL --connect-timeout 10 --max-time 60 -o \"{}\" \"{}\"",
                             tmp_path.string(), url);
           std::error_code ec;
           if (!RunShellCommand(cmd) || !std::filesystem::exists(tmp_path)) {
             std::filesystem::remove(tmp_path, ec);
             return std::nullopt;
           }
           std::ifstream ifs(tmp_path);
           std::stringstream ss;
           ss << ifs.rdbuf();
           ifs.close();
           std::filesystem::remove(tmp_path, ec);
           return ss.str();
         }
   
         std::optional<std::string> FindLatestFilenameInListing(const std::string& listing_url,
                                                                const std::regex& pattern) {
           auto contents = FetchTextFromUrl(listing_url);
           if (!contents.has_value()) return std::nullopt;
   
           std::string latest;
           for (auto it = std::sregex_iterator(contents->begin(), contents->end(), pattern);
                it != std::sregex_iterator(); ++it) {
             const std::string match = it->str();
             if (match > latest) latest = match;
           }
           if (latest.empty()) return std::nullopt;
           return latest;
         }
   
         std::optional<std::string> FindLatestFilenameLocally(const std::filesystem::path& dir,
                                                              const std::regex& pattern) {
           std::string latest;
           std::error_code ec;
           for (const auto& entry : std::filesystem::directory_iterator(dir, ec)) {
             if (!entry.is_regular_file()) continue;
             std::string name = entry.path().filename().string();
             if (std::regex_match(name, pattern) && name > latest) latest = name;
           }
           if (latest.empty()) return std::nullopt;
           return latest;
         }
   
         bool RefreshAndLoadKernel(const std::string& filename, const std::string& url,
                                   const std::filesystem::path& kernel_dir) {
           std::filesystem::path local_path = kernel_dir / filename;
           bool had_local_copy = std::filesystem::exists(local_path);
   
           if (DownloadToFile(url, local_path)) {
             Logger::Info(fmt::format("SPICE: downloaded latest NAIF kernel '{}'", filename), "Spice");
           } else if (had_local_copy) {
             Logger::Debug(
                 fmt::format("SPICE: could not refresh '{}' from NAIF; using cached local copy",
                             filename),
                 "Spice");
           } else {
             Logger::Warn(fmt::format("SPICE: could not download '{}' from NAIF ({}) and no local "
                                      "copy is available; this kernel will not be loaded",
                                      filename, url),
                          "Spice");
             return false;
           }
   
           furnsh_c(filename.c_str());
           CheckSpiceFailure(fmt::format("Loading {}", filename));
           return true;
         }
       }  // namespace
   
       void LoadSpiceKernel(void) {
         if (spice_loaded) return;
         SpiceInt kcount;
         ktotal_c("ALL", &kcount);
         if (kcount > 0) return;
   
         SetSpiceReturnMode();
   
         std::string orig_dir = std::filesystem::current_path().string();
         std::filesystem::path kernel_dir = GetCspiceKernelDir();
         std::filesystem::current_path(kernel_dir);
   
         furnsh_c("naif0012.tls");  // leap seconds
         CheckSpiceFailure("Loading naif0012.tls");
         furnsh_c("de440.bsp");  // planetary ephemeris
         CheckSpiceFailure("Loading de440.bsp");
         furnsh_c("pck00011.tpc");  // planetary constants
         CheckSpiceFailure("Loading pck00011.tpc");
   
         // ----------------------------------------------------------------
         // High-accuracy Earth and lunar orientation kernels
         //
         // The default text PCK ("pck00011.tpc") only provides the low-accuracy
         // IAU_EARTH / IAU_MOON body-fixed frames (errors of roughly 0.06 deg /
         // 6 km for Earth, and up to ~0.005 deg / ~155 m for the Moon).
         // High-accuracy orientation requires binary PCKs that define the
         // ITRF93 (Earth) and MOON_PA_DExxx / MOON_ME_DExxx (Moon) frames; see
         // the NAIF tutorial "'High Accuracy' Orientation and Body-fixed Frames
         // for the Moon and Earth" (April 2023).
         //
         // NAIF continually republishes these kernels as new orientation data
         // becomes available, so during initialization LuPNT looks up and
         // downloads the *current latest* files directly from the NAIF generic
         // kernels server, caching them under the local kernel directory and
         // falling back to the cached copy (or older bundled kernels) if the
         // server cannot be reached. Set LUPNT_SKIP_SPICE_KERNEL_DOWNLOAD=1 to
         // force fully-offline operation using only locally-cached kernels.
         // ----------------------------------------------------------------
         const std::string kPckUrl = std::string(kNaifGenericKernelsUrl) + "/pck/";
         const std::string kLunarFkUrl = std::string(kNaifGenericKernelsUrl) + "/fk/satellites/";
   
         // --- Earth --------------------------------------------------------
         // NAIF publishes the latest reconstructed + short-term-predicted
         // high-accuracy Earth orientation under a stable file name that never
         // changes ("earth_latest_high_prec.bpc"), making it trivial to always
         // fetch the most recent data.
         const std::string kEarthLatestPck = "earth_latest_high_prec.bpc";
         bool earth_loaded
             = RefreshAndLoadKernel(kEarthLatestPck, kPckUrl + kEarthLatestPck, kernel_dir);
         if (!earth_loaded) {
           // Fall back to whichever (older, hard-coded) Earth orientation
           // kernels may already be cached locally from a previous LuPNT
           // release, so that LuPNT keeps working fully offline.
           for (const char* fallback :
                {"earth_200101_990825_predict.bpc", "earth_000101_241014_240722.bpc"}) {
             if (std::filesystem::exists(fallback)) {
               furnsh_c(fallback);
               CheckSpiceFailure(fmt::format("Loading {}", fallback));
             }
           }
         }
   
         // --- Moon ---------------------------------------------------------
         // Unlike Earth, NAIF does not publish the latest lunar orientation
         // kernels under a stable file name: a new binary PCK + frame kernel
         // (FK) pair is released -- named after the underlying JPL DE/LE
         // ephemeris model and release date -- each time a new ephemeris is
         // produced (e.g. "moon_pa_de440_200625.bpc" + "moon_de440_250416.tf").
         // Determine the current latest pair by inspecting the NAIF directory
         // listings: file names embed zero-padded DE version numbers and
         // dates, so the lexicographically-largest match is also the most
         // recent. If the listings cannot be retrieved, fall back to the
         // newest matching kernel already cached locally.
         static const std::regex kMoonPckPattern(R"(moon_pa_de\d+_[0-9]{6}\.bpc)");
         static const std::regex kMoonFkPattern(R"(moon_de\d+_[0-9]{6}\.tf)");
   
         auto moon_pck_name = FindLatestFilenameInListing(kPckUrl, kMoonPckPattern);
         if (!moon_pck_name.has_value())
           moon_pck_name = FindLatestFilenameLocally(kernel_dir, kMoonPckPattern);
   
         auto moon_fk_name = FindLatestFilenameInListing(kLunarFkUrl, kMoonFkPattern);
         if (!moon_fk_name.has_value())
           moon_fk_name = FindLatestFilenameLocally(kernel_dir, kMoonFkPattern);
   
         // Load the binary PCK before the frame kernel: the FK connects the
         // MOON_PA_DExxx frame name to the orientation data provided by the PCK.
         if (moon_pck_name.has_value()) {
           RefreshAndLoadKernel(*moon_pck_name, kPckUrl + *moon_pck_name, kernel_dir);
         } else {
           Logger::Warn(
               "SPICE: could not determine the latest lunar orientation binary PCK from "
               "NAIF; high-accuracy lunar orientation (MOON_PA / MOON_ME) will be "
               "unavailable",
               "Spice");
         }
         if (moon_fk_name.has_value()) {
           RefreshAndLoadKernel(*moon_fk_name, kLunarFkUrl + *moon_fk_name, kernel_dir);
         } else {
           Logger::Warn(
               "SPICE: could not determine the latest lunar frame kernel (FK) from NAIF; "
               "high-accuracy lunar orientation (MOON_PA / MOON_ME) will be unavailable",
               "Spice");
         }
   
         // Make the high-accuracy MOON_PA frame the default lunar body-fixed
         // frame for the legacy (pre-N0062) SPICE APIs (LSPCN, ET2LST, ILLUM,
         // SRFXPT, SUBPT, SUBSOL, ...).
         const std::string kMoonAssocPa = "moon_assoc_pa.tf";
         RefreshAndLoadKernel(kMoonAssocPa, kLunarFkUrl + kMoonAssocPa, kernel_dir);
   
         // Mars
         if (std::filesystem::exists("mars097.bsp")) {
           furnsh_c("mars097.bsp");
           CheckSpiceFailure("Loading mars097.bsp");
         }
   
         // Load Chebyshev coefficients
         cheby_s = spk_extract("de440.bsp", &cheby_n);
         if (cheby_s == nullptr) {
           throw std::runtime_error(
               "Could not load SPK file - Please Download the SPK file. See "
               "data/ephemeris/readme.md for instructions");
         }
   
         std::filesystem::current_path(orig_dir);
         spice_loaded = true;
         return;
       }
   
       void LoadSpiceKernel(const std::filesystem::path& filepath) {
         LoadSpiceKernel();
         SetSpiceReturnMode();
         furnsh_c(filepath.string().c_str());
         CheckSpiceFailure(fmt::format("Loading {}", filepath.string()));
       }
   
       int GetNaifId(const std::string& body_name) {
         LoadSpiceKernel();
         SpiceInt code;
         SpiceBoolean found;
   #pragma omp critical
         {
           bodn2c_c(body_name.c_str(), &code, &found);
         }
         CheckSpiceFailure(fmt::format("Resolving NAIF body {}", body_name));
         if (!found) throw std::runtime_error(fmt::format("NAIF body not found: {}", body_name));
         return static_cast<int>(code);
       }
   
       std::string GetNaifName(int naif_id) {
         LoadSpiceKernel();
         SpiceChar name[256];
         SpiceBoolean found;
   #pragma omp critical
         {
           bodc2n_c(static_cast<SpiceInt>(naif_id), sizeof(name), name, &found);
         }
         CheckSpiceFailure(fmt::format("Resolving NAIF body id {}", naif_id));
         if (!found) return std::to_string(naif_id);
         return std::string(name);
       }
   
       bool HasNaifBody(const std::string& body_name) {
         LoadSpiceKernel();
         SpiceInt code;
         SpiceBoolean found;
   #pragma omp critical
         {
           bodn2c_c(body_name.c_str(), &code, &found);
         }
         CheckSpiceFailure(fmt::format("Resolving NAIF body {}", body_name));
         return static_cast<bool>(found);
       }
   
       bool HasFrame(const std::string& frame_name) {
         LoadSpiceKernel();
         SpiceInt frame_code = 0;
   #pragma omp critical
         {
           namfrm_c(frame_name.c_str(), &frame_code);
         }
         CheckSpiceFailure(fmt::format("Resolving NAIF frame {}", frame_name));
         return frame_code != 0;
       }
   
       void ExtractPckCoeffs() {
         int handle;
         SpiceInt pck_handle;
         ConstSpiceChar* pck_file = "../data/ephemeris/moon_pa_de440_200625.bpc";
         double t_tdb = 8000.0;
         int body = 301;
         double descr[5];
         char ident[40];
         int found;
         SpiceDouble record[120];
         int rsize;
         int pdeg;
         // double ra;
         // double dec;
         // double w;
         // double lambda;
         static char bref[32];
         double eulang[6];
   
         // The fundamental quantities defined by PCK orientation models are actually
         // Euler angles, not matrices. These Euler angles, which we call ``RA, DEC,
         // and W,'' are related to the transformation operator returned from
         // pxform_c by the equation rotate = [ W ]   [ Pi/2 - DEC ]   [ Pi/2 + RA ]
         //              3                1               3
         // To directly retrieve these angles, use the call:
   
         // bodeul_( &body, &t_tdb, &ra, &dec, &w, &lambda );
         // std::cout << "body: " << body << std::endl;
         // std::cout << "t_tdb: " << t_tdb << std::endl;
         // std::cout << "ra: " << ra << std::endl;
         // std::cout << "dec: " << dec << std::endl;
         // std::cout << "w: " << w << std::endl;
         // std::cout << " " << std::endl;
   
         pckeul_(&body, &t_tdb, &found, bref, eulang, (ftnlen)32);
         std::cout << "found:" << found << std::endl;
         std::cout << "phi: " << eulang[0] << std::endl;
         std::cout << "delta: " << eulang[1] << std::endl;
         std::cout << " " << std::endl;
   
         pcklof_c(pck_file, &pck_handle);  // load the PCK file
         pcksfs_(&body, &t_tdb, &handle, descr, ident, &found, (ftnlen)40);
   
         std::cout << "pck handle: :" << pck_handle << std::endl;
         std::cout << "handle: :" << handle << std::endl;
         std::cout << "descr: " << &descr << std::endl;
         std::cout << "ident: " << &ident << std::endl;
         std::cout << "found:" << found << std::endl;
   
         if (found) {
           pckr02_(&handle, descr, &t_tdb, record);
           rsize = record[1];
           pdeg = (rsize - 2) / 3 - 1;
           std::cout << "Polynomial Size:" << rsize << std::endl;
           std::cout << "Polynomial Degree:" << pdeg << std::endl;
         }
   
         // extract coefficients from the CSPICE PCK file
         //    pcksfs_(body, t_tdb, &handle, descr, ident, found, (ftnlen)40);
         //    pckr02_c(handle, target)
       }
   
       Vec3d GetBodyPosSpice(Real t_tdb, BodyId obs, BodyId target, const std::string& refFrame,
                             const std::string& abCorrection) {
         if (!spice_loaded) {
           LoadSpiceKernel();
         }
   
         std::string targ_str = std::to_string((int)target);
         std::string obs_str = std::to_string((int)obs);
   
         SpiceDouble ptarg[3];
         SpiceDouble et = t_tdb.val();
         const char* targ = strcpy(new char[targ_str.length() + 1], targ_str.c_str());
         const char* ref = strcpy(new char[refFrame.length() + 1], refFrame.c_str());
         const char* abcorr = strcpy(new char[abCorrection.length() + 1], abCorrection.c_str());
         const char* obs_spice = strcpy(new char[obs_str.length() + 1], obs_str.c_str());
         SpiceDouble lt;
   
         // void spkpos_c(ConstSpiceChar * targ, SpiceDouble t_tdb, ConstSpiceChar *
         // ref, ConstSpiceChar * abcorr,
         //               ConstSpiceChar * obs, SpiceDouble ptarg[3], SpiceDouble *
         //               lt)
         spkpos_c(targ, et, ref, abcorr, obs_spice, ptarg, &lt);
   
         Vec3d r;
         for (int i = 0; i < 3; i++) r(i) = ptarg[i];
         return r;
       }
   
       Vec6d GetBodyPosVelSpice(Real t_tdb, BodyId obs, BodyId target, const std::string& refFrame,
                                const std::string& abCorrection) {
         if (!spice_loaded) LoadSpiceKernel();
   
         SpiceDouble starg[6];
         Vec6d rv;
   
         std::string targ_str = std::to_string((int)target);
         std::string obs_str = std::to_string((int)obs);
   
         SpiceDouble et = t_tdb.val();
         const char* targ = strcpy(new char[targ_str.length() + 1], targ_str.c_str());
         const char* ref = strcpy(new char[refFrame.length() + 1], refFrame.c_str());
         const char* abcorr = strcpy(new char[abCorrection.length() + 1], abCorrection.c_str());
         const char* obs_spice = strcpy(new char[obs_str.length() + 1], obs_str.c_str());
         SpiceDouble lt;
   
         //  void spkez_c ( SpiceInt            targ,
         //                 SpiceDouble         et,
         //                 ConstSpiceChar     *ref,
         //                 ConstSpiceChar     *abcorr,
         //                 SpiceInt            obs,
         //                 SpiceDouble         starg[6],
         //                 SpiceDouble        *lt        )
   #pragma omp critical
         {
           spkezr_c(targ, et, ref, abcorr, obs_spice, starg, &lt);
         }
         for (int i = 0; i < 6; i++) rv(i) = starg[i];
         return rv;
       }
   
       GroundStationSpiceData GetGroundStationDataSpice(Real t_tdb, const std::string& station_name,
                                                        BodyId center, const std::string& refFrame,
                                                        const std::string& abCorrection) {
         int station_id = GetNaifId(station_name);
         GroundStationSpiceData data
             = GetGroundStationDataSpice(t_tdb, station_id, center, refFrame, abCorrection);
         data.name = station_name;
         return data;
       }
   
       GroundStationSpiceData GetGroundStationDataSpice(Real t_tdb, int station_id, BodyId center,
                                                        const std::string& refFrame,
                                                        const std::string& abCorrection) {
         LoadSpiceKernel();
         SetSpiceReturnMode();
   
         SpiceDouble et = t_tdb.val();
         SpiceDouble pos_km[3];
         SpiceDouble lt;
         std::string station = std::to_string(station_id);
         std::string center_id = std::to_string(static_cast<int>(center));
   
   #pragma omp critical
         {
           spkpos_c(station.c_str(), et, refFrame.c_str(), abCorrection.c_str(), center_id.c_str(),
                    pos_km, &lt);
         }
         CheckSpiceFailure(fmt::format("Fetching ground station {} from SPICE", station_id));
   
         GroundStationSpiceData data;
         data.naif_id = station_id;
         data.name = GetNaifName(station_id);
         data.body_id = center;
         data.frame = refFrame;
         for (int i = 0; i < 3; ++i) data.position_m(i) = pos_km[i] * M_KM;
   
         BodyData body_data = GetBodyData(center);
         Cart3 pos(data.position_m.cast<Real>(), body_data.fixed_frame);
         LatLonAlt lla = CartToLatLonAlt(pos, body_data.R, body_data.flattening);
         data.latitude_deg = lla(0).val();
         data.longitude_deg = lla(1).val();
         data.altitude_m = lla(2).val();
         return data;
       }
   
       Mat6d GetFrameConversionMat(Real t_tdb, const std::string& from_frame,
                                   const std::string& to_frame) {
         if (!spice_loaded) LoadSpiceKernel();
   
         SpiceDouble et_spice = (SpiceDouble)t_tdb.val();
         double xform[6][6];
         Mat6d M_rot;
   
         const char* from_frame_char = strcpy(new char[from_frame.length() + 1], from_frame.c_str());
         const char* to_frame_char = strcpy(new char[to_frame.length() + 1], to_frame.c_str());
   #pragma omp critical
         {
           sxform_c(from_frame_char, to_frame_char, et_spice, xform);
         }
   
         for (int i = 0; i < 6; i++) {
           for (int j = 0; j < 6; j++) {
             M_rot(i, j) = xform[i][j];
           }
         }
   
         return M_rot;
       }
   
       Vec3d GetPlanetOrientation(BodyId id, Real t_tdb) {
         if (!spice_loaded) LoadSpiceKernel();
   
         std::string to_frame = "J2000";
   
         char bodyname[36];
         SpiceBoolean found;
   #pragma omp critical
         {
           bodc2n_c((SpiceInt)id, 36, bodyname, &found);
         }
         if (!found) {
           throw std::runtime_error("Invalid planet ID");
         }
         std::string from_frame = "IAU_" + std::string(bodyname);
   
         SpiceDouble rotmat[3][3];
         SpiceDouble et_spice = (SpiceDouble)t_tdb.val();
         const char* from_frame_char = strcpy(new char[from_frame.length() + 1], from_frame.c_str());
         const char* to_frame_char = strcpy(new char[to_frame.length() + 1], to_frame.c_str());
   #pragma omp critical
         {
           pxform_c(from_frame_char, to_frame_char, et_spice, rotmat);
         }
         SpiceDouble psi, theta, phi;
   
   #pragma omp critical
         m2eul_c(rotmat, 3, 1, 3, &psi, &theta, &phi);
   
         Vec3d angles;
         double W = -(phi + PI);          // R_z(-W)
         double delta0 = PI / 2 - theta;  // R_x(PI/2 - delta0)
         double alpha0 = PI / 2 - psi;    // R_z(-alpha0)
         angles << alpha0, delta0, W;
         return angles;
       }
   
       Real StringToTdb(const std::string& str) {
         if (!spice_loaded) LoadSpiceKernel();
         SpiceDouble t_tdb;
         str2et_c(str.c_str(), &t_tdb);
         return t_tdb;
       }
   
       Real StringToTai(const std::string& str) {
         if (!spice_loaded) LoadSpiceKernel();
         Real t_tdb = StringToTdb(str);
         Real t_tai = spice::ConvertTime(t_tdb, Time::TDB, Time::TAI);
         return t_tai;
       }
   
       std::string TDBtoStringUTC(Real t_tdb, int prec = 3) {
         if (!spice_loaded) LoadSpiceKernel();
         SpiceDouble et = t_tdb.val();
         SpiceChar str[100];
         et2utc_c(et, "C", prec, 100, str);
         return std::string(str);
       }
   
       std::string TAItoStringUTC(Real t_tai, int prec = 3) {
         if (!spice_loaded) LoadSpiceKernel();
         Real et_tdb = spice::ConvertTime(t_tai, Time::TAI, Time::TDB);
         std::string str = TDBtoStringUTC(et_tdb, prec);
         return str;
       }
   
       Real ConvertTime(Real t, Time from, Time to) {
         if (!spice_loaded) LoadSpiceKernel();
         if (from == to) return t;
         SpiceDouble t_in = t.val();
         SpiceDouble t_out_spice;
   #pragma omp critical
         {
           t_out_spice
               = unitim_c(t_in, time_to_string.at(from).c_str(), time_to_string.at(to).c_str());
         }
         double offset = t_out_spice - t_in;  // offset in seconds
         Real t_out = t + offset;             // this is to convert to real
         return t_out;
       }
   
       MatX6 GetBodyPosVel(const VecX& t_tdb, BodyId center, BodyId target) {
         if (!spice_loaded) LoadSpiceKernel();
         MatX6 retState(t_tdb.size(), 6);
         for (int i = 0; i < t_tdb.size(); i++)
           retState.row(i) = GetBodyPosVel(t_tdb(i), center, target).transpose();
         return retState;
       }
   
       Vec6 GetBodyPosVelBase(const Real t_tdb, BodyId center, BodyId target) {
         if (!spice_loaded) LoadSpiceKernel();
         if (center == target) return Vec6::Zero();
   
         for (int i = 0; i < cheby_n; i++) {
           if (cheby_s[i].center == (int)center && cheby_s[i].target == (int)target)
             return cheby_posvel_ad(t_tdb, cheby_s[i].seg, cheby_s[i].len);
           if (cheby_s[i].center == (int)target && cheby_s[i].target == (int)center)
             return -cheby_posvel_ad(t_tdb, cheby_s[i].seg, cheby_s[i].len);
         }
         LUPNT_CHECK(false, "Chebyshev coefficients not found", "SpiceInterface");
         return Vec6::Zero();
       }
   
       Vec6 GetBodyPosVel(const Real t_tdb, BodyId center, BodyId target) {
         if (!spice_loaded) LoadSpiceKernel();
         if (center == target) return Vec6::Zero();
   
         Vec6 rv = Vec6::Zero();
         if (target == BodyId::EARTH) {
           rv += GetBodyPosVelBase(t_tdb, BodyId::EMB, BodyId::EARTH);
           target = BodyId::EMB;
         }
         if (center == BodyId::EARTH) {
           rv -= GetBodyPosVelBase(t_tdb, BodyId::EMB, BodyId::EARTH);
           center = BodyId::EMB;
         }
         if (target == BodyId::MOON) {
           rv += GetBodyPosVelBase(t_tdb, BodyId::EMB, BodyId::MOON);
           target = BodyId::EMB;
         }
         if (center == BodyId::MOON) {
           rv -= GetBodyPosVelBase(t_tdb, BodyId::EMB, BodyId::MOON);
           center = BodyId::EMB;
         }
         if (target == BodyId::MERCURY) {
           rv += GetBodyPosVelBase(t_tdb, BodyId::MERCURY_BARYCENTER, BodyId::MERCURY);
           target = BodyId::MERCURY_BARYCENTER;
         }
         if (center == BodyId::MERCURY) {
           rv -= GetBodyPosVelBase(t_tdb, BodyId::MERCURY_BARYCENTER, BodyId::MERCURY);
           center = BodyId::MERCURY_BARYCENTER;
         }
         if (target == BodyId::VENUS) {
           rv += GetBodyPosVelBase(t_tdb, BodyId::VENUS_BARYCENTER, BodyId::VENUS);
           target = BodyId::VENUS_BARYCENTER;
         }
         if (center == BodyId::VENUS) {
           rv -= GetBodyPosVelBase(t_tdb, BodyId::VENUS_BARYCENTER, BodyId::VENUS);
           center = BodyId::VENUS_BARYCENTER;
         }
         rv += GetBodyPosVelBase(t_tdb, center, target);
         return rv;
       }
   
     }  // namespace spice
   
   }  // namespace lupnt