From c39baac3a86a97d8eba8d7137ea83b21ef55e79a Mon Sep 17 00:00:00 2001 From: Malin Ejdbo Date: Mon, 21 Jun 2021 13:43:31 +0200 Subject: [PATCH] Add functions to convert a cartesian coordinate to a ra dec distance --- include/openspace/util/coordinateconversion.h | 18 ++++ src/scene/scene.cpp | 11 ++- src/scene/scene_lua.inl | 23 ++++- src/util/coordinateconversion.cpp | 97 ++++++++++++++++++- 4 files changed, 140 insertions(+), 9 deletions(-) diff --git a/include/openspace/util/coordinateconversion.h b/include/openspace/util/coordinateconversion.h index b0e8ba0f46..d33a46c2ae 100644 --- a/include/openspace/util/coordinateconversion.h +++ b/include/openspace/util/coordinateconversion.h @@ -49,6 +49,24 @@ glm::dvec3 icrsToGalacticCartesian(double ra, double dec, double distance); */ glm::dvec2 icrsToDecimalDegrees(const std::string& ra, const std::string& dec); +/** + * Converts from galactic cartesian coordinates to ICRS coordinates in decimal degrees + * and distance. + * \param x X coordinate + * \param y Y coordinate + * \param z Z coordinate + * \return A vector with the ra and dec decimal degrees and distance. + */ +glm::dvec3 galacticCartesianToIcrs(double x, double y, double z); + +/** + * Converts from ICRS (decimal degrees) coordinates to ICRS (hms and dms) coordinates. + * \param ra Right ascension, given in decimal degrees + * \param dec Declination, given in decimal degrees + * \return A pair with the ra and dec strings in hms and dms format. + */ +std::pair decimalDegreesToIcrs(double ra, double dec); + } // namespace openspace #endif // __OPENSPACE_CORE___COORDINATECONVERSION___H__ diff --git a/src/scene/scene.cpp b/src/scene/scene.cpp index 75506d0acb..c8e23f9573 100644 --- a/src/scene/scene.cpp +++ b/src/scene/scene.cpp @@ -733,11 +733,18 @@ scripting::LuaLibrary Scene::luaLibrary() { "Returns the world rotation matrix of the scene graph node with the given string as identifier" }, { - "convertRaDec", - &luascriptfunctions::convertRaDec, + "convertFromRaDec", + &luascriptfunctions::convertFromRaDec, {}, "string, string, double", "Returns the cartesian world position of a ra dec coordinate with distance" + }, + { + "convertToRaDec", + &luascriptfunctions::convertToRaDec, + {}, + "double, double, double", + "Returns the ra, dec strings and distance for a given cartesian world coordinate" } } }; diff --git a/src/scene/scene_lua.inl b/src/scene/scene_lua.inl index 036bdf5b9a..bdaf0c1f7c 100644 --- a/src/scene/scene_lua.inl +++ b/src/scene/scene_lua.inl @@ -950,8 +950,8 @@ int worldRotation(lua_State* L) { return 1; } -int convertRaDec(lua_State* L) { - ghoul::lua::checkArgumentsAndThrow(L, 3, "lua::convertRaDec"); +int convertFromRaDec(lua_State* L) { + ghoul::lua::checkArgumentsAndThrow(L, 3, "lua::convertFromRaDec"); std::string ra = ghoul::lua::value(L, 1); std::string dec = ghoul::lua::value(L, 2); @@ -967,4 +967,23 @@ int convertRaDec(lua_State* L) { return 1; } +int convertToRaDec(lua_State* L) { + ghoul::lua::checkArgumentsAndThrow(L, 3, "lua::convertToRaDec"); + + double x = ghoul::lua::value(L, 1); + double y = ghoul::lua::value(L, 2); + double z = ghoul::lua::value(L, 3); + lua_settop(L, 0); + + glm::dvec3 degrees = galacticCartesianToIcrs(x, y, z); + std::pair raDecPair = decimalDegreesToIcrs(degrees.x, degrees.y); + + ghoul::lua::push(L, raDecPair.first); // Ra + ghoul::lua::push(L, raDecPair.second); // Dec + ghoul::lua::push(L, degrees.z); // Distance + + ghoul_assert(lua_gettop(L) == 3, "Incorrect number of items left on stack"); + return 3; +} + } // namespace openspace::luascriptfunctions diff --git a/src/util/coordinateconversion.cpp b/src/util/coordinateconversion.cpp index 2d3ce2db52..ea64bf374b 100644 --- a/src/util/coordinateconversion.cpp +++ b/src/util/coordinateconversion.cpp @@ -174,7 +174,7 @@ namespace openspace { // into Galactic coordinates (l, b) glm::dvec3 icrsToGalacticCartesian(double ra, double dec, double distance) { // Reference: - // https://www.atnf.csiro.au/people/Tobias.Westmeier/tools_coords.php, + // https://www.atnf.csiro.au/people/Tobias.Westmeier/tools_coords.php // (Ra, Dec) -> (a, d) double a = glm::radians(ra); @@ -202,9 +202,11 @@ glm::dvec3 icrsToGalacticCartesian(double ra, double dec, double distance) { return distance * rGalactic; } -// ra format "XXhYYmZZs" or "XhYmZs" -// dec format "XXdYYmZZs", "-XXdYYmZZs", "XdYmZs" or "-XdYmZs" +// Ra format "XXhYYmZZs" or "XhYmZs" +// Dec format "XXdYYmZZs", "-XXdYYmZZs", "XdYmZs" or "-XdYmZs" glm::dvec2 icrsToDecimalDegrees(const std::string& ra, const std::string& dec) { + // Reference: + // https://math.stackexchange.com/questions/15323/how-do-i-calculate-the-cartesian-coordinates-of-stars if (ra.size() > 9 || ra.size() < 6 || dec.size() > 10 || dec.size() < 6) { throw(ghoul::lua::LuaRuntimeException(fmt::format( "Ra '{}' or Dec '{}' format is incorrect. Correct format is: Ra 'XXhYYmZZs' " @@ -224,8 +226,8 @@ glm::dvec2 icrsToDecimalDegrees(const std::string& ra, const std::string& dec) { if (!isRaDecValid(ra_hours, ra_minutes, ra_seconds, dec_degrees, dec_minutes, dec_seconds)) { - LWARNING(fmt::format("Ra '{}' or Dec '{}' is outside the allowed range", - ra, dec) + LWARNING(fmt::format("Ra '{}' or Dec '{}' is outside the allowed range, " + "result may be incorrect", ra, dec) ); } @@ -242,4 +244,89 @@ glm::dvec2 icrsToDecimalDegrees(const std::string& ra, const std::string& dec) { return glm::dvec2(ra_deg, dec_deg); } +// Convert Galactic coordinates (x, y, z) or (l, b) into +// Equatorial coordinates ICRS right ascension and declination (a, d) plus distance +glm::dvec3 galacticCartesianToIcrs(double x, double y, double z) { + // References: + // https://www.atnf.csiro.au/people/Tobias.Westmeier/tools_coords.php, + // https://en.wikipedia.org/wiki/Celestial_coordinate_system + + // Normalize + double distance = sqrt(x * x + y * y + z * z); + double n_x = x / distance; + double n_y = y / distance; + double n_z = z / distance; + + // Convert from cartesian + // (x, y, z) -> (l, b) + double l = atan2(n_y, n_x); + double b = asin(n_z); + + // J2000 Galactic reference frame + constexpr double a0 = glm::radians(192.8595); // Equatorial coordinates of the Galactic north pole + constexpr double d0 = glm::radians(27.1284); + constexpr double l0 = glm::radians(122.9320); // Galactic longitude of the equatorial north pole + + // Convert to equatorial reference frame + double a = atan2( + cos(b) * sin(l0 - l), + sin(b) * cos(d0) - cos(b) * sin(d0) * cos(l0 - l) + ) + a0; + double d = asin(sin(b) * sin(d0) + cos(b) * cos(d0) * cos(l0 - l)); + + glm::dvec3 rEquatorial = glm::dvec3(a, d, distance); + + return rEquatorial; +} + +// Return a pair with two formatted strings from the decimal degrees ra and dec +std::pair decimalDegreesToIcrs(double ra, double dec) { + // References: + // https://www.rapidtables.com/convert/number/degrees-to-degrees-minutes-seconds.html, + // https://math.stackexchange.com/questions/15323/how-do-i-calculate-the-cartesian-coordinates-of-stars + + // Radians to degrees + double ra_deg = glm::degrees(ra); + double dec_deg = glm::degrees(dec); + + // Check input + if (ra_deg < 0 || ra_deg > 360 || dec_deg < -90 || dec_deg > 90) { + LWARNING(fmt::format("Given Ra '{}' or Dec '{}' is outside the allowed range, " + "result may be incorrect", ra, dec) + ); + } + + // Calculate Ra + int ra_hours = std::trunc(ra_deg) / 15.0; + double ra_minutes_full = (ra_deg - ra_hours * 15.0) * 60.0 / 15.0; + int ra_minutes = std::trunc(ra_minutes_full); + double ra_seconds = (ra_minutes_full - ra_minutes) * 60.0; + + // Calculate Dec + int dec_degrees = std::trunc(dec_deg); + double dec_minutes_full = (abs(dec_deg) - abs(dec_degrees)) * 60.0; + int dec_minutes = std::trunc(dec_minutes_full); + double dec_seconds = (dec_minutes_full - dec_minutes) * 60.0; + + // Construct strings + std::pair result; + result.first = std::to_string(ra_hours) + 'h' + + std::to_string(ra_minutes) + 'm' + + std::to_string(ra_seconds) + 's'; + + result.second = std::to_string(dec_degrees) + 'd' + + std::to_string(dec_minutes) + 'm' + + std::to_string(dec_seconds) + 's'; + + // Check results + if (!isRaDecValid(ra_hours, ra_minutes, ra_seconds, dec_degrees, dec_minutes, + dec_seconds)) + { + LWARNING(fmt::format("Resulting Ra '{}' or Dec '{}' is outside the allowed range, " + "result may be incorrect", result.first, result.second) + ); + } + return result; +} + } // namespace openspace