Merge pull request #1671 from OpenSpace/feature/radec-conversion

Feature/radec conversion
This commit is contained in:
Malin E
2021-07-05 10:31:16 +02:00
committed by GitHub
6 changed files with 424 additions and 16 deletions

View File

@@ -26,18 +26,45 @@
#define __OPENSPACE_CORE___COORDINATECONVERSION___H__
#include <ghoul/glm.h>
#include <string>
namespace openspace {
/**
* Converts from ICRS coordinates to galactic cartesian coordinates.
* Converts from ICRS decimal degrees coordinates to galactic cartesian coordinates.
* \param ra Right ascension, given in decimal degrees
* \param dec Declination, given in decimal degrees
* \param distance The distance, or radius, to the position given in any unit.
* \return A position in galactic cartesian coordinates, given in the same unit as the
* distance parameter.
*/
glm::dvec3 icrsToGalacticCartesian(float ra, float dec, double distance);
glm::dvec3 icrsToGalacticCartesian(double ra, double dec, double distance);
/**
* Converts from ICRS (hms and dms) coordinates to decimal degrees.
* \param ra Right ascension, given as a string in format 'XhYmZs'
* \param dec Declination, given as a string in format 'XdYmZs'
* \return The decimal degrees coordinate in degrees
*/
glm::dvec2 icrsToDecimalDegrees(const std::string& ra, const std::string& dec);
/**
* Converts from galactic cartesian coordinates to ICRS decimal degrees coordinates
* and distance.
* \param x X coordinate
* \param y Y coordinate
* \param z Z coordinate
* \return A vector with the ra and dec decimal degrees in 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<std::string, std::string> decimalDegreesToIcrs(double ra, double dec);
} // namespace openspace

View File

@@ -44,6 +44,7 @@ set(HEADER_FILES
source_group("Header Files" FILES ${HEADER_FILES})
set(SOURCE_FILES
spacemodule_lua.inl
speckloader.cpp
rendering/planetgeometry.cpp
rendering/renderableconstellationbounds.cpp

View File

@@ -39,11 +39,14 @@
#include <openspace/documentation/documentation.h>
#include <openspace/rendering/renderable.h>
#include <openspace/rendering/screenspacerenderable.h>
#include <openspace/scripting/lualibrary.h>
#include <openspace/util/factorymanager.h>
#include <openspace/util/spicemanager.h>
#include <ghoul/misc/assert.h>
#include <ghoul/misc/templatefactory.h>
#include "spacemodule_lua.inl"
namespace {
constexpr openspace::properties::Property::PropertyInfo SpiceExceptionInfo = {
"ShowExceptions",
@@ -133,4 +136,31 @@ std::vector<documentation::Documentation> SpaceModule::documentations() const {
};
}
scripting::LuaLibrary SpaceModule::luaLibrary() const {
scripting::LuaLibrary res;
res.name = "space";
res.functions = {
{
"convertFromRaDec",
&space::luascriptfunctions::convertFromRaDec,
{},
"string/double, string/double, double",
"Returns the cartesian world position of a ra dec coordinate with distance. "
"If the coordinate is given as strings the format should be ra 'XhYmZs' and "
"dec 'XdYmZs'. If the coordinate is given as numbers the values should be "
"in degrees."
},
{
"convertToRaDec",
&space::luascriptfunctions::convertToRaDec,
{},
"double, double, double",
"Returns the formatted ra, dec strings and distance for a given cartesian "
"world coordinate."
}
};
return res;
}
} // namespace openspace

View File

@@ -42,6 +42,8 @@ public:
static ghoul::opengl::ProgramObjectManager ProgramObjectManager;
scripting::LuaLibrary luaLibrary() const override;
private:
void internalInitialize(const ghoul::Dictionary&) override;
void internalDeinitializeGL() override;

View File

@@ -0,0 +1,78 @@
/*****************************************************************************************
* *
* OpenSpace *
* *
* Copyright (c) 2014-2021 *
* *
* Permission is hereby granted, free of charge, to any person obtaining a copy of this *
* software and associated documentation files (the "Software"), to deal in the Software *
* without restriction, including without limitation the rights to use, copy, modify, *
* merge, publish, distribute, sublicense, and/or sell copies of the Software, and to *
* permit persons to whom the Software is furnished to do so, subject to the following *
* conditions: *
* *
* The above copyright notice and this permission notice shall be included in all copies *
* or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, *
* INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A *
* PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT *
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF *
* CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE *
* OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. *
****************************************************************************************/
#include <openspace/util/coordinateconversion.h>
namespace openspace::space::luascriptfunctions {
int convertFromRaDec(lua_State* L) {
ghoul::lua::checkArgumentsAndThrow(L, 3, "lua::convertFromRaDec");
glm::dvec2 degrees = glm::dvec2(0.0);
if (lua_type(L, 1) == LUA_TSTRING && lua_type(L, 2) == LUA_TSTRING) {
std::string ra = ghoul::lua::value<std::string>(L, 1);
std::string dec = ghoul::lua::value<std::string>(L, 2);
degrees = icrsToDecimalDegrees(ra, dec);
}
else if (lua_type(L, 1) == LUA_TNUMBER && lua_type(L, 2) == LUA_TNUMBER) {
degrees.x = ghoul::lua::value<double>(L, 1);
degrees.y = ghoul::lua::value<double>(L, 2);
}
else {
throw ghoul::lua::LuaRuntimeException("lua::convertFromRaDec: Ra and Dec have to "
"be of the same type, either String or Number"
);
}
double distance = ghoul::lua::value<double>(L, 3);
lua_settop(L, 0);
glm::dvec3 pos = icrsToGalacticCartesian(degrees.x, degrees.y, distance);
ghoul::lua::push(L, pos);
ghoul_assert(lua_gettop(L) == 1, "Incorrect number of items left on stack");
return 1;
}
int convertToRaDec(lua_State* L) {
ghoul::lua::checkArgumentsAndThrow(L, 3, "lua::convertToRaDec");
double x = ghoul::lua::value<double>(L, 1);
double y = ghoul::lua::value<double>(L, 2);
double z = ghoul::lua::value<double>(L, 3);
lua_settop(L, 0);
glm::dvec3 degrees = galacticCartesianToIcrs(x, y, z);
std::pair<std::string, std::string> 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::space::luascriptfunctions

View File

@@ -24,26 +24,296 @@
#include <openspace/util/coordinateconversion.h>
#include <ghoul/fmt.h>
#include <ghoul/logging/logmanager.h>
#include <ghoul/lua/ghoul_lua.h>
namespace {
constexpr const char* _loggerCat = "Coordinateconversion";
// 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
void parseString(const std::string& str, int& hoursOrDegrees, int& minutes,
double& seconds)
{
// Find hms or dms indicies
size_t hOrDIndex =
(str.find('h') != std::string::npos) ? str.find('h') : str.find('d');
size_t mIndex = str.find('m');
size_t sIndex = str.find('s');
if (hOrDIndex == std::string::npos || mIndex == std::string::npos ||
sIndex == std::string::npos)
{
throw ghoul::lua::LuaRuntimeException(fmt::format(
"Ra or Dec '{}' format is incorrect. Correct format is: Ra 'XhYmZs', "
"and Dec 'XdYmZs'", str));
}
// Construct the number strings
std::string sHoursOrDegrees = str.substr(0, hOrDIndex);
std::string sMinutes = str.substr(hOrDIndex + 1, mIndex - hOrDIndex - 1);
std::string sSeconds = str.substr(mIndex + 1, sIndex - mIndex - 1);
// Convert the strings to numbers
try {
// Hours or degrees must be an integer
double temp = std::stod(sHoursOrDegrees);
if (std::floor(temp) != temp) {
throw ghoul::lua::LuaRuntimeException(fmt::format(
"Ra or Dec '{}' format is incorrect. Correct format is: Ra 'XhYmZs', "
"and Dec 'XdYmZs', where X must be an integer", str));
}
hoursOrDegrees = std::stoi(sHoursOrDegrees);
// Minutes must be an integer
temp = std::stod(sMinutes);
if (std::floor(temp) != temp) {
throw ghoul::lua::LuaRuntimeException(fmt::format(
"Ra or Dec '{}' format is incorrect. Correct format is: Ra 'XhYmZs', "
"and Dec 'XdYmZs', where Y must be an integer", str));
}
minutes = std::stoi(sMinutes);
// Seconds is a double
seconds = std::stod(sSeconds);
} catch (const std::invalid_argument& ia) {
throw ghoul::lua::LuaRuntimeException(fmt::format(
"Ra or Dec '{}' format is incorrect. Correct format is: Ra 'XhYmZs', "
"and Dec 'XdYmZs'", str));
}
}
void parseRa(const std::string& ra, int& hours, int& minutes, double& seconds) {
if (ra.find('d') != std::string::npos) {
throw ghoul::lua::LuaRuntimeException(fmt::format(
"Ra '{}' format is incorrect. Correct format is: 'XhYmZs'", ra));
}
parseString(ra, hours, minutes, seconds);
}
void parseDec(const std::string& dec, int& degrees, int& minutes,
double& seconds)
{
if (dec.find('h') != std::string::npos) {
throw ghoul::lua::LuaRuntimeException(fmt::format(
"Dec '{}' format is incorrect. Correct format is: 'XdYmZs'", dec));
}
parseString(dec, degrees, minutes, seconds);
}
bool isRaDecValid(int raH, int raM, double raS, int decD,
int decM, double decS)
{
// Ra
if (raH < 0.0 || raH >= 24.0) {
LWARNING(fmt::format("Right ascension hours '{}' is outside the allowed "
"range of 0 to 24 hours (exclusive)", raH)
);
return false;
}
if (raM < 0.0 || raM >= 60.0) {
LWARNING(fmt::format("Right ascension minutes '{}' is outside the allowed "
"range of 0 to 60 minutes (exclusive)", raM)
);
return false;
}
if (raS < 0.0 || raS >= 60.0) {
LWARNING(fmt::format("Right ascension seconds '{}' is outside the allowed "
"range of 0 to 60 seconds (exclusive)", raS)
);
return false;
}
// Dec
if (decD < -90.0 || decD > 90.0) {
LWARNING(fmt::format("Declination degrees '{}' is outside the allowed range "
"of -90 to 90 degrees (inclusive)", decD)
);
return false;
}
else if ((decD == -90.0 || decD == 90.0) && (decM != 0 || decS != 0)) {
LWARNING("Total declination is outside the allowed range of -90 to 90 "
"degrees (inclusive)"
);
return false;
}
if (decM < 0.0 || decM >= 60.0) {
LWARNING(fmt::format("Declination minutes '{}' is outside the allowed range "
"of 0 to 60 minutes (exclusive)", decM)
);
return false;
}
if (decS < 0.0 || decS >= 60.0) {
LWARNING(fmt::format("Declination seconds '{}' is outside the allowed range "
"of 0 to 60 seconds (exclusive)", decS)
);
return false;
}
return true;
}
} // namespace
namespace openspace {
glm::dvec3 icrsToGalacticCartesian(float ra, float dec, double distance) {
// Convert to Galactic Coordinates from ICRS right ascension and declination
// https://gea.esac.esa.int/archive/documentation/GDR2/Data_processing/
// chap_cu3ast/sec_cu3ast_intro/ssec_cu3ast_intro_tansforms.html#SSS1
const glm::dmat3 conversionMatrix = glm::dmat3({
-0.0548755604162154, 0.4941094278755837, -0.8676661490190047, // col 0
-0.8734370902348850, -0.4448296299600112, -0.1980763734312015, // col 1
-0.4838350155487132, 0.7469822444972189, 0.4559837761750669 // col 2
});
// Convert Equatorial coordinates ICRS right ascension and declination (a, d)
// into Galactic coordinates (l, b)
// Reference:
// https://www.atnf.csiro.au/people/Tobias.Westmeier/tools_coords.php
glm::dvec3 icrsToGalacticCartesian(double ra, double dec, double distance) {
// (Ra, Dec) -> (a, d)
double a = glm::radians(ra);
double d = glm::radians(dec);
glm::dvec3 rICRS = glm::dvec3(
cos(glm::radians(ra)) * cos(glm::radians(dec)),
sin(glm::radians(ra)) * cos(glm::radians(dec)),
sin(glm::radians(dec))
// Convert to galactic reference frame
double l = L0 - atan2(
cos(d) * sin(a - A0),
sin(d) * cos(D0) - cos(d) * sin(D0) * cos(a - A0)
);
double b = asin(sin(d) * sin(D0) + cos(d) * cos(D0) * cos(a - A0));
// Convert to cartesian
glm::dvec3 rGalactic = glm::dvec3(
cos(b) * cos(l),
cos(b) * sin(l),
sin(b)
);
glm::dvec3 rGalactic = conversionMatrix * rICRS; // on the unit sphere
return distance * rGalactic;
}
// Ra format 'XhYmZs', where X and Y are positive integers and Z is a positive double
// Dec format 'XdYmZs', where X is a signed integer, Y is a positive integer and Z is a
// positive double
// Reference:
// https://math.stackexchange.com/questions/15323/how-do-i-calculate-the-cartesian-coordinates-of-stars
glm::dvec2 icrsToDecimalDegrees(const std::string& ra, const std::string& dec) {
if (ra.size() < 6 || dec.size() < 6) {
throw ghoul::lua::LuaRuntimeException(fmt::format(
"Ra '{}' or Dec '{}' format is incorrect. Correct format is: Ra 'XhYmZs', "
"and Dec 'XdYmZs'", ra, dec));
}
// Parse right ascension
int raHours, raMinutes;
double raSeconds;
parseRa(ra, raHours, raMinutes, raSeconds);
// Parse declination
int decDegrees, decMinutes;
double decSeconds;
parseDec(dec, decDegrees, decMinutes, decSeconds);
const bool isValid = isRaDecValid(raHours,
raMinutes,
raSeconds,
decDegrees,
decMinutes,
decSeconds
);
if (!isValid) {
LWARNING(fmt::format("Ra '{}' or Dec '{}' is outside the allowed range, "
"result may be incorrect", ra, dec)
);
}
// Convert from hours/degrees, minutes, seconds to decimal degrees
double sign = signbit(static_cast<float>(decDegrees)) ? -1.0 : 1.0;
double raDeg = (raHours * 15.0) +
(raMinutes * 15.0 / 60.0) +
(raSeconds * 15.0 / 3600.0);
double decDeg = (abs(decDegrees) +
(decMinutes / 60.0) +
(decSeconds / 3600.0)) * sign;
return glm::dvec2(raDeg, decDeg);
}
// Convert Galactic coordinates (x, y, z) or (l, b) into Equatorial coordinates ICRS
// right ascension and declination in decimal degrees (a, d) plus distance
// References:
// https://www.atnf.csiro.au/people/Tobias.Westmeier/tools_coords.php,
// https://en.wikipedia.org/wiki/Celestial_coordinate_system
glm::dvec3 galacticCartesianToIcrs(double x, double y, double z) {
// Normalize
double distance = sqrt(x*x + y*y + z*z);
double nX = x / distance;
double nY = y / distance;
double nZ = z / distance;
// Convert from cartesian
// (x, y, z) -> (l, b)
double l = atan2(nY, nX);
double b = asin(nZ);
// 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));
return glm::dvec3(glm::degrees(a), glm::degrees(d), distance);
}
// Return a pair with two formatted strings from the decimal degrees ra and 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
std::pair<std::string, std::string> decimalDegreesToIcrs(double ra, double dec) {
// Radians to degrees
double raDeg = ra;
double decDeg = dec;
// Check input
if (raDeg < 0 || raDeg > 360 || decDeg < -90 || decDeg > 90) {
LWARNING(fmt::format("Given Ra '{}' or Dec '{}' is outside the allowed range, "
"result may be incorrect", ra, dec)
);
}
// Calculate Ra
int raHours = std::trunc(raDeg) / 15.0;
double raMinutesFull = (raDeg - raHours * 15.0) * 60.0 / 15.0;
int raMinutes = std::trunc(raMinutesFull);
double raSeconds = (raMinutesFull - raMinutes) * 60.0;
// Calculate Dec
int decDegrees = std::trunc(decDeg);
double decMinutesFull = (abs(decDeg) - abs(decDegrees)) * 60.0;
int decMinutes = std::trunc(decMinutesFull);
double decSeconds = (decMinutesFull - decMinutes) * 60.0;
// Construct strings
std::pair<std::string, std::string> result;
result.first = std::to_string(raHours) + 'h' +
std::to_string(raMinutes) + 'm' +
std::to_string(raSeconds) + 's';
result.second = std::to_string(decDegrees) + 'd' +
std::to_string(decMinutes) + 'm' +
std::to_string(decSeconds) + 's';
// Check results
const bool isValid = isRaDecValid(raHours,
raMinutes,
raSeconds,
decDegrees,
decMinutes,
decSeconds
);
if (!isValid) {
LWARNING(fmt::format("Resulting Ra '{}' or Dec '{}' is outside the allowed range, "
"result may be incorrect", result.first, result.second)
);
}
return result;
}
} // namespace openspace