ktyl
/
skein
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

feat: use doubles in lib

This commit is contained in:
Cat Flynn 2024-08-21 14:51:37 +01:00
parent c4edf127e9
commit ae3e5c201b
4 changed files with 87 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
lib/skein/include/skein/orbit.h
View File

@ -4,7 +4,7 @@
#include "glm/glm.hpp" #include "glm/glm.hpp"
typedef std::vector<float> Vector6; typedef std::vector<double> Vector6;
class Orbit class Orbit
{ {
@ -12,28 +12,34 @@ public:
Orbit(); Orbit();
~Orbit() = default; ~Orbit() = default;
float getSemiMajorAxis() const; double getSemiMajorAxis() const;
float getEccentricity() const; double getEccentricity() const;
float getInclination() const; double getInclination() const;
float getArgumentOfPeriapsis() const; double getArgumentOfPeriapsis() const;
float getLongitudeOfAscendingNode() const; double getLongitudeOfAscendingNode() const;
void setSemiMajorAxis(float semiMajorAxis); void setSemiMajorAxis(double semiMajorAxis);
void setEccentricity(float eccentricity); void setEccentricity(double eccentricity);
void setInclination(float inclination); void setInclination(double inclination);
void setArgumentOfPeriapsis(float argumentOfPeriapsis); void setArgumentOfPeriapsis(double argumentOfPeriapsis);
void setLongitudeOfAscendingNode(float longitudeOfAscendingNode); void setLongitudeOfAscendingNode(double longitudeOfAscendingNode);
// TODO: meanAnomaly in all these arguments actually means eccentricMeanAnomaly, // TODO: meanAnomaly in all these arguments actually means ellipticalMeanAnomaly,
// will have to change that when adding non-ellipctical orbits - don't get confused! // will have to change that when adding non-ellipctical orbits - don't get confused!
glm::vec3 getPosition(const float meanAnomaly) const; //
glm::vec3 getTangent(const float meanAnomaly) const; // do we actual want to use mean anomaly for this? we could instead pass in time and
glm::mat4 getLookAlongMatrix(const float meanAnomaly) const; // the gravitational parameter which are easier numbers to come by, i think? at least,
// the mean motion is determined from the mean motion and time time
//glm::vec3 getPosition(double gravitationalParameter, double time) const;
glm::dvec3 getPosition(const double meanAnomaly) const;
glm::dvec3 getTangent(const double meanAnomaly) const;
//glm::mat4 getLookAlongMatrix(const double meanAnomaly) const;
private: private:
Vector6 _keplerianElements; Vector6 _keplerianElements;
const float getEccentricAnomaly(const float meanAnomaly) const; const double getEccentricAnomaly(const double meanAnomaly) const;
const float getTrueAnomaly(const float meanAnomaly) const; const double getTrueAnomaly(const double meanAnomaly) const;
const Vector6 getCartesianCoordinates(const float meanAnomaly) const; const Vector6 getCartesianCoordinates(const double meanAnomaly) const;
}; };

2
lib/skein/include/skein/particlemap.h
View File

@ -14,7 +14,7 @@ class ParticleMap
// providing these as two methods keeps things unambiguous - manipulating particles // providing these as two methods keeps things unambiguous - manipulating particles
// is just done for setup // is just done for setup
const Particle& getParticle(const std::string& id) const; const Particle& getParticle(const std::string& id) const;
glm::vec3 getParticlePosition(const std::string& id, double time) const; glm::dvec3 getParticlePosition(const std::string& id, double time) const;
const Orbit& getOrbit(const std::string& id) const; const Orbit& getOrbit(const std::string& id) const;
void setParticle(const Particle& particle); void setParticle(const Particle& particle);

90
lib/skein/src/orbit.cpp
View File

@ -10,115 +10,115 @@ Orbit::Orbit()
_keplerianElements.resize(6); _keplerianElements.resize(6);
} }
float Orbit::getSemiMajorAxis() const double Orbit::getSemiMajorAxis() const
{ {
return _keplerianElements[astro::semiMajorAxisIndex]; return _keplerianElements[astro::semiMajorAxisIndex];
} }
void Orbit::setSemiMajorAxis(float semiMajorAxis) void Orbit::setSemiMajorAxis(double semiMajorAxis)
{ {
_keplerianElements[astro::semiMajorAxisIndex] = semiMajorAxis; _keplerianElements[astro::semiMajorAxisIndex] = semiMajorAxis;
} }
float Orbit::getEccentricity() const double Orbit::getEccentricity() const
{ {
return _keplerianElements[astro::eccentricityIndex]; return _keplerianElements[astro::eccentricityIndex];
} }
void Orbit::setEccentricity(float eccentricity) void Orbit::setEccentricity(double eccentricity)
{ {
_keplerianElements[astro::eccentricityIndex] = eccentricity; _keplerianElements[astro::eccentricityIndex] = eccentricity;
} }
float Orbit::getInclination() const double Orbit::getInclination() const
{ {
return _keplerianElements[astro::inclinationIndex]; return _keplerianElements[astro::inclinationIndex];
} }
void Orbit::setInclination(float inclination) void Orbit::setInclination(double inclination)
{ {
_keplerianElements[astro::inclinationIndex] = inclination; _keplerianElements[astro::inclinationIndex] = inclination;
} }
float Orbit::getArgumentOfPeriapsis() const double Orbit::getArgumentOfPeriapsis() const
{ {
return _keplerianElements[astro::argumentOfPeriapsisIndex]; return _keplerianElements[astro::argumentOfPeriapsisIndex];
} }
void Orbit::setArgumentOfPeriapsis(float argumentOfPeriapsis) void Orbit::setArgumentOfPeriapsis(double argumentOfPeriapsis)
{ {
_keplerianElements[astro::argumentOfPeriapsisIndex] = argumentOfPeriapsis; _keplerianElements[astro::argumentOfPeriapsisIndex] = argumentOfPeriapsis;
} }
float Orbit::getLongitudeOfAscendingNode() const double Orbit::getLongitudeOfAscendingNode() const
{ {
return _keplerianElements[astro::longitudeOfAscendingNodeIndex]; return _keplerianElements[astro::longitudeOfAscendingNodeIndex];
} }
void Orbit::setLongitudeOfAscendingNode(float longitudeOfAscendingNode) void Orbit::setLongitudeOfAscendingNode(double longitudeOfAscendingNode)
{ {
_keplerianElements[astro::longitudeOfAscendingNodeIndex] = longitudeOfAscendingNode; _keplerianElements[astro::longitudeOfAscendingNodeIndex] = longitudeOfAscendingNode;
} }
glm::mat4 Orbit::getLookAlongMatrix(const float meanAnomaly) const //glm::mat4 Orbit::getLookAlongMatrix(const double meanAnomaly) const
//{
// glm::vec3 position = getPosition(meanAnomaly);
//
// // get the tangent of the orbital ellipse
// glm::vec3 tan = getTangent(meanAnomaly);
// // we want to point along the orbit
// glm::vec3 target = position + tan;
// // TODO: this is not 'up' with respect to the orbited body!
// // 'up' is just the normalized position vector because the orbit is centred at the origin
// glm::vec3 up = glm::normalize(position);
// // easy peasy lookAt matrix
// glm::mat4 look = glm::lookAt(position, target, up);
//
// // invert the lookat matrix because it's meant for cameras, cameras work backwards and
// // we are not a camera
// glm::mat4 lookAlong = glm::inverse(look);
//
// return lookAlong;
//}
const double Orbit::getEccentricAnomaly(const double meanAnomaly) const
{ {
glm::vec3 position = getPosition(meanAnomaly); const double eccentricity = _keplerianElements[astro::eccentricityIndex];
double eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS(
// get the tangent of the orbital ellipse
glm::vec3 tan = getTangent(meanAnomaly);
// we want to point along the orbit
glm::vec3 target = position + tan;
// TODO: this is not 'up' with respect to the orbited body!
// 'up' is just the normalized position vector because the orbit is centred at the origin
glm::vec3 up = glm::normalize(position);
// easy peasy lookAt matrix
glm::mat4 look = glm::lookAt(position, target, up);
// invert the lookat matrix because it's meant for cameras, cameras work backwards and
// we are not a camera
glm::mat4 lookAlong = glm::inverse(look);
return lookAlong;
}
const float Orbit::getEccentricAnomaly(const float meanAnomaly) const
{
const float eccentricity = _keplerianElements[astro::eccentricityIndex];
float eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS(
eccentricity, eccentricity,
meanAnomaly, meanAnomaly,
100); 200);
return eccentricAnomaly; return eccentricAnomaly;
} }
// Interpolate a position around the orbit. // Interpolate a position around the orbit.
// t is in range 0..1 and wraps. // t is in range 0..1 and wraps.
glm::vec3 Orbit::getPosition(const float meanAnomaly) const glm::dvec3 Orbit::getPosition(const double meanAnomaly) const
{ {
Vector6 cartesian = getCartesianCoordinates(meanAnomaly); Vector6 cartesian = getCartesianCoordinates(meanAnomaly);
return glm::vec3( return glm::dvec3(
cartesian[astro::xPositionIndex], cartesian[astro::xPositionIndex],
cartesian[astro::yPositionIndex], cartesian[astro::yPositionIndex],
cartesian[astro::zPositionIndex]); cartesian[astro::zPositionIndex]);
} }
const float Orbit::getTrueAnomaly(const float meanAnomaly) const const double Orbit::getTrueAnomaly(const double meanAnomaly) const
{ {
const float eccentricAnomaly = getEccentricAnomaly(meanAnomaly); const double eccentricAnomaly = getEccentricAnomaly(meanAnomaly);
const float eccentricity = _keplerianElements[astro::eccentricityIndex]; const double eccentricity = _keplerianElements[astro::eccentricityIndex];
return astro::convertEccentricAnomalyToTrueAnomaly( return astro::convertEccentricAnomalyToTrueAnomaly(
eccentricAnomaly, eccentricAnomaly,
eccentricity); eccentricity);
} }
const Vector6 Orbit::getCartesianCoordinates(const float meanAnomaly) const const Vector6 Orbit::getCartesianCoordinates(const double meanAnomaly) const
{ {
Vector6 kepler(_keplerianElements); Vector6 kepler(_keplerianElements);
kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly); kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly);
return astro::convertKeplerianToCartesianElements(kepler, 1.0); return astro::convertKeplerianToCartesianElements(kepler, 1.0);
} }
glm::vec3 Orbit::getTangent(const float meanAnomaly) const glm::dvec3 Orbit::getTangent(const double meanAnomaly) const
{ {
float epsilon = 0.01; double epsilon = 0.01;
glm::vec3 ahead = getPosition(meanAnomaly + epsilon); glm::dvec3 ahead = getPosition(meanAnomaly + epsilon);
glm::vec3 behind = getPosition(meanAnomaly - epsilon); glm::dvec3 behind = getPosition(meanAnomaly - epsilon);
return glm::normalize(ahead - behind); return glm::normalize(ahead - behind);
} }

26
lib/skein/src/particlemap.cpp
View File

@ -12,17 +12,25 @@ const Orbit& ParticleMap::getOrbit(const std::string& id) const
return _orbits.at(id); return _orbits.at(id);
} }
glm::vec3 ParticleMap::getParticlePosition(const std::string& id, double time) const glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time) const
{ {
if (_orbits.find(id) == _orbits.end())
return {0,0,0};
// TODO: actually get stuff based on physics
//const Particle& parent = _relationships.at(id);
// how do we get the gravitational parameter?
// is it just G * M
// where M is the mass of the parent
//const double u = getGravitationalParameter(parent.getId());
//const double mass = parent.getMass();
// how do we get the particle position from the time + gravitational parameter?
// TODO: actually nest stuff so position is determined from all parents // TODO: actually nest stuff so position is determined from all parents
const Orbit& orbit = _orbits.at(id);
if (_orbits.find(id) != _orbits.end()) return orbit.getPosition(time);
{
const Orbit& orbit = _orbits.at(id);
return orbit.getPosition(time);
}
return {0,0,0};
} }
void ParticleMap::setParticle(const Particle& particle) void ParticleMap::setParticle(const Particle& particle)