diff --git a/lib/skein/include/skein/orbit.h b/lib/skein/include/skein/orbit.h
index 89000df..050977c 100644
--- a/lib/skein/include/skein/orbit.h
+++ b/lib/skein/include/skein/orbit.h
@@ -17,15 +17,18 @@ public:
double getInclination() const;
double getArgumentOfPeriapsis() const;
double getLongitudeOfAscendingNode() const;
+ double getTrueAnomaly() const;
void setSemiMajorAxis(double semiMajorAxis);
void setEccentricity(double eccentricity);
void setInclination(double inclination);
void setArgumentOfPeriapsis(double argumentOfPeriapsis);
void setLongitudeOfAscendingNode(double longitudeOfAscendingNode);
+ void setTrueAnomaly(double trueAnomaly);
- glm::dvec3 getPosition(double gravitationalParameter, double time) const;
- glm::dvec3 getPositionFromMeanAnomaly(double meanAnomaly, double gravitationalParameter) const;
+ glm::dvec3 getPosition(double gravitationalParameter) const;
+ //glm::dvec3 getPosition(double gravitationalParameter, double time) const;
+ //glm::dvec3 getPositionFromMeanAnomaly(double meanAnomaly, double gravitationalParameter) const;
glm::dvec3 getVelocity(double gravitationalParameter, double time) const;
glm::dvec3 getVelocityFromMeanAnomaly(double meanAnomaly, double gravitationalParameter) const;
@@ -33,12 +36,23 @@ public:
//glm::dvec3 getTangent(const double meanAnomaly) const;
//glm::mat4 getLookAlongMatrix(const double meanAnomaly) const;
- const Vector6 getCartesianCoordinates(const double meanAnomaly, const double gravitationalParameter) const;
+ const Vector6 getCartesianCoordinates(const double gravitationalParameter) const;
const double getMeanAnomaly(const double gravitationalParameter, const double time) const;
+ void update(double time, double gravitationalParameter);
+
private:
+ /*
+ * N.B.: Order of elements is important!
+ * keplerianElements(0) = semiMajorAxis [m]
+ * keplerianElements(1) = eccentricity [-]
+ * keplerianElements(2) = inclination [rad]
+ * keplerianElements(3) = argument of periapsis [rad]
+ * keplerianElements(4) = longitude of [rad]
+ * ascending node [rad]
+ * keplerianElements(5) = true anomaly [rad]
+ */
Vector6 _keplerianElements;
- const double getEccentricAnomaly(const double meanAnomaly) const;
- const double getTrueAnomaly(const double meanAnomaly) const;
+ const double getEccentricAnomaly() const;
};
diff --git a/lib/skein/include/skein/particlemap.h b/lib/skein/include/skein/particlemap.h
index afbb1b5..bc5c55f 100644
--- a/lib/skein/include/skein/particlemap.h
+++ b/lib/skein/include/skein/particlemap.h
@@ -15,18 +15,22 @@ class ParticleMap
// is just done for setup
const Particle& getParticle(const std::string& id) const;
const Particle& getParent(const std::string& id) const;
- glm::dvec3 getParticleLocalPosition(const std::string& id, double time) const;
- glm::dvec3 getParticlePosition(const std::string& id, double time) const;
- glm::dvec3 getParticleLocalVelocity(const std::string& id, double time) const;
const Orbit& getOrbit(const std::string& id) const;
void setParticle(const Particle& particle);
void setRelationship(const std::string& parentId, const std::string& childId, const Orbit& orbit);
- void setParticleLocalVelocity(const std::string& id, const double& time, const glm::dvec3& velocity);
+ //void update(double time);
+ void setParticleLocalVelocity(const std::string& id, double time, const glm::dvec3& velocity);
+
+ glm::dvec3 getParticlePosition(const std::string& id) const;
+ glm::dvec3 getParticleLocalPosition(const std::string& id, double time) const;
+ glm::dvec3 getParticleLocalVelocity(const std::string& id, double time) const;
private:
std::map _particles;
// TODO: improve on this very vague name
std::map _relationships;
std::map _orbits;
+
+ //double _time;
};
diff --git a/lib/skein/src/orbit.cpp b/lib/skein/src/orbit.cpp
index c251162..de2d632 100644
--- a/lib/skein/src/orbit.cpp
+++ b/lib/skein/src/orbit.cpp
@@ -56,6 +56,15 @@ void Orbit::setLongitudeOfAscendingNode(double longitudeOfAscendingNode)
_keplerianElements[astro::longitudeOfAscendingNodeIndex] = longitudeOfAscendingNode;
}
+double Orbit::getTrueAnomaly() const
+{
+ return _keplerianElements[astro::trueAnomalyIndex];
+}
+void Orbit::setTrueAnomaly(double trueAnomaly)
+{
+ _keplerianElements[astro::trueAnomalyIndex] = trueAnomaly;
+}
+
//glm::mat4 Orbit::getLookAlongMatrix(const double meanAnomaly) const
//{
// glm::vec3 position = getPosition(meanAnomaly);
@@ -77,34 +86,48 @@ void Orbit::setLongitudeOfAscendingNode(double longitudeOfAscendingNode)
// return lookAlong;
//}
-const double Orbit::getEccentricAnomaly(const double meanAnomaly) const
+const double Orbit::getEccentricAnomaly() const
{
- const double eccentricity = _keplerianElements[astro::eccentricityIndex];
- double eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS(
- eccentricity,
- meanAnomaly,
- 200);
+ // TODO: get the eccentric anomaly from the true anomaly?
+//Real convertTrueAnomalyToEllipticalEccentricAnomaly(const Real trueAnomaly,
+ //const Real eccentricity)
+
+ const double trueAnomaly = getTrueAnomaly();
+ const double eccentricity = getEccentricity();
+ //double eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS(
+ // eccentricity,
+ // meanAnomaly,
+ // 200);
+ double eccentricAnomaly = astro::convertTrueAnomalyToEllipticalEccentricAnomaly(trueAnomaly, eccentricity);
return eccentricAnomaly;
}
-glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly, double gravitationalParameter) const
+//glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly, double gravitationalParameter) const
+//{
+// Vector6 cartesian = getCartesianCoordinates(meanAnomaly, gravitationalParameter);
+// return glm::dvec3(
+// cartesian[astro::xPositionIndex],
+// cartesian[astro::yPositionIndex],
+// cartesian[astro::zPositionIndex]);
+//}
+//
+//glm::dvec3 Orbit::getPosition(double gravitationalParameter, double time) const
+//{
+// return getPositionFromMeanAnomaly(getMeanAnomaly(gravitationalParameter, time), gravitationalParameter);
+//}
+glm::dvec3 Orbit::getPosition(double gravitationalParameter) const
{
- Vector6 cartesian = getCartesianCoordinates(meanAnomaly, gravitationalParameter);
+ const Vector6 cartesian = getCartesianCoordinates(gravitationalParameter);
return glm::dvec3(
cartesian[astro::xPositionIndex],
cartesian[astro::yPositionIndex],
cartesian[astro::zPositionIndex]);
}
-glm::dvec3 Orbit::getPosition(double gravitationalParameter, double time) const
-{
- return getPositionFromMeanAnomaly(getMeanAnomaly(gravitationalParameter, time), gravitationalParameter);
-}
-
glm::dvec3 Orbit::getVelocityFromMeanAnomaly(const double meanAnomaly, double gravitationalParameter) const
{
- Vector6 cartesian = getCartesianCoordinates(meanAnomaly, gravitationalParameter);
+ Vector6 cartesian = getCartesianCoordinates(gravitationalParameter);
return glm::dvec3(
cartesian[astro::xVelocityIndex],
cartesian[astro::yVelocityIndex],
@@ -122,20 +145,41 @@ const double Orbit::getMeanAnomaly(const double gravitationalParameter, const do
return meanMotion * time;
}
-const double Orbit::getTrueAnomaly(const double meanAnomaly) const
+//const double Orbit::getTrueAnomaly(const double meanAnomaly) const
+//{
+// const double eccentricAnomaly = getEccentricAnomaly(meanAnomaly);
+// const double eccentricity = _keplerianElements[astro::eccentricityIndex];
+// return astro::convertEccentricAnomalyToTrueAnomaly(
+// eccentricAnomaly,
+// eccentricity);
+//}
+
+const Vector6 Orbit::getCartesianCoordinates(const double gravitationalParameter) const
{
- const double eccentricAnomaly = getEccentricAnomaly(meanAnomaly);
- const double eccentricity = _keplerianElements[astro::eccentricityIndex];
- return astro::convertEccentricAnomalyToTrueAnomaly(
- eccentricAnomaly,
- eccentricity);
+ //Vector6 kepler(_keplerianElements);
+ //kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly);
+ return astro::convertKeplerianToCartesianElements(_keplerianElements, gravitationalParameter);
}
-const Vector6 Orbit::getCartesianCoordinates(const double meanAnomaly, const double gravitationalParameter) const
+void Orbit::update(double time, double gravitationalParameter)
{
- Vector6 kepler(_keplerianElements);
- kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly);
- return astro::convertKeplerianToCartesianElements(kepler, gravitationalParameter);
+ // compute mean motion from time
+ // we assume true anomaly is 0 at t=0 so we this motion is the mean anomaly
+ const double semiMajorAxis = getSemiMajorAxis();
+ const double meanMotion = astro::computeKeplerMeanMotion(semiMajorAxis, gravitationalParameter);
+ const double meanAnomaly = meanMotion * time;
+
+ // compute elliptical eccentric anomaly from mean anomaly
+ const double eccentricity = getEccentricity();
+ const double eccentricAnomaly = astro::convertEllipticalMeanAnomalyToEccentricAnomalyBS(
+ eccentricity,
+ meanAnomaly,
+ 200);
+
+ // compute true anomaly from eccentric anomaly
+ const double trueAnomaly = astro::convertEccentricAnomalyToTrueAnomaly(eccentricAnomaly, eccentricity);
+
+ setTrueAnomaly(trueAnomaly);
}
//glm::dvec3 Orbit::getTangent(const double meanAnomaly) const
diff --git a/lib/skein/src/particlemap.cpp b/lib/skein/src/particlemap.cpp
index 20ed187..f9086af 100644
--- a/lib/skein/src/particlemap.cpp
+++ b/lib/skein/src/particlemap.cpp
@@ -32,12 +32,13 @@ glm::dvec3 ParticleMap::getParticleLocalPosition(const std::string& id, double t
const double gravitationalParameter = parent.getGravitationalParameter();
const Orbit& orbit = _orbits.at(id);
- pos += orbit.getPosition(gravitationalParameter, time);
+ //pos += orbit.getPosition(gravitationalParameter, time);
+ pos += orbit.getPosition(gravitationalParameter);
return pos;
}
-glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time) const
+glm::dvec3 ParticleMap::getParticlePosition(const std::string& id) const
{
if (_orbits.find(id) == _orbits.end())
return {0,0,0};
@@ -54,7 +55,8 @@ glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time)
const double gravitationalParameter = parent.getGravitationalParameter();
const Orbit& orbit = _orbits.at(id);
- pos += orbit.getPosition(gravitationalParameter, time);
+ //pos += orbit.getPosition(gravitationalParameter, time);
+ pos += orbit.getPosition(gravitationalParameter);
auto it = _relationships.find(parentId);
if (it != _relationships.end())
@@ -89,7 +91,7 @@ glm::dvec3 ParticleMap::getParticleLocalVelocity(const std::string& id, double t
return vel;
}
-void ParticleMap::setParticleLocalVelocity(const std::string& id, const double& time, const glm::dvec3& velocity)
+void ParticleMap::setParticleLocalVelocity(const std::string& id, double time, const glm::dvec3& velocity)
{
// we want to convert a local position and velocity into a new set of keplerian elements.
@@ -99,7 +101,8 @@ void ParticleMap::setParticleLocalVelocity(const std::string& id, const double&
const Particle& parent = _particles.at(parentId);
const double gravitationalParameter = parent.getGravitationalParameter();
const double meanAnomaly = orbit.getMeanAnomaly(gravitationalParameter, time);
- Vector6 cartesian(orbit.getCartesianCoordinates(meanAnomaly, gravitationalParameter));
+ //Vector6 cartesian(orbit.getCartesianCoordinates(meanAnomaly, gravitationalParameter));
+ Vector6 cartesian(orbit.getCartesianCoordinates(gravitationalParameter));
// update velocity, leave position as it is
cartesian[astro::xVelocityIndex] = velocity.x;
@@ -130,3 +133,13 @@ void ParticleMap::setRelationship(const std::string& parentId, const std::string
_relationships.insert_or_assign(childId, parentId);
_orbits.insert_or_assign(childId, orbit);
}
+
+//void ParticleMap::update(double time)
+//{
+// //const double interval = time - _time;
+// //_time = time;
+//
+// // what should happen when we update the time?
+//
+// // assume that all particles' true anomaly is 0 at t=0
+//}
diff --git a/src/main.cpp b/src/main.cpp
index 990c2a9..0048622 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -94,6 +94,7 @@ int main()
moonOrbit.setInclination(glm::radians(5.15)); // radians
moonOrbit.setArgumentOfPeriapsis(318.15); // in the case of the moon these last two values are
moonOrbit.setLongitudeOfAscendingNode(60.0); // pretty much constantly changing so use whatever
+ moonOrbit.setTrueAnomaly(0.0);
// set parameters of satellite orbit around moon
Orbit stationOrbit;
@@ -102,6 +103,7 @@ int main()
stationOrbit.setInclination(glm::radians(89.0));
stationOrbit.setArgumentOfPeriapsis(43.2);
stationOrbit.setLongitudeOfAscendingNode(239.7);
+ stationOrbit.setTrueAnomaly(0.0);
ParticleMap map;
map.setParticle({"earth", 5.9e24});
@@ -140,11 +142,17 @@ int main()
// apply input
if (input.pauseTime)
{
- mode++;
- if (mode > MODE_MANEOUVRE)
+ if (mode == MODE_MANEOUVRE)
{
- mode = 0;
+ mode = MODE_ORBIT;
+ std::cout << "play";
}
+ else if (mode == MODE_ORBIT)
+ {
+ mode = MODE_MANEOUVRE;
+ std::cout << "pause";
+ }
+ std::cout << std::endl;
}
// the moon takes like 27 days to orbit earth. to see this in motion, then, we need to
@@ -167,7 +175,7 @@ int main()
// TODO: what unit is this in lol
const glm::dvec3 v = map.getParticleLocalVelocity("station", time);
- const glm::dvec3 newVelocity = v * 0.9;
+ const glm::dvec3 newVelocity = v * 0.99;
printVector(v);
std::cout << " -> ";
@@ -182,6 +190,16 @@ int main()
glClearColor(0.2, 0.3, 0.3, 1.0);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+ // update orbits with time
+ const Particle& earth = map.getParticle("earth");
+ const Particle& moon = map.getParticle("moon");
+ moonOrbit = map.getOrbit("moon");
+ stationOrbit = map.getOrbit("station");
+ moonOrbit.update(time, earth.getGravitationalParameter());
+ stationOrbit.update(time, moon.getGravitationalParameter());
+ map.setRelationship("earth", "moon", moonOrbit);
+ map.setRelationship("moon", "station", stationOrbit);
+
earthVis.render(time);
moonVis.render(time);
stationVis.render(time);
diff --git a/src/orbitvisualizer.cpp b/src/orbitvisualizer.cpp
index 79bf4ea..a45c4a1 100644
--- a/src/orbitvisualizer.cpp
+++ b/src/orbitvisualizer.cpp
@@ -12,7 +12,7 @@ void OrbitVisualizer::render(const float time)
{
// Orbit needs to be drawn relative to its parent
const Particle& parent = _map.getParent(_particleId);
- const glm::vec3 parentPos = _map.getParticlePosition(parent.getId(), time);
+ const glm::vec3 parentPos = _map.getParticlePosition(parent.getId());
regenerateVertices(parentPos);
@@ -27,7 +27,7 @@ void OrbitVisualizer::render(const float time)
void OrbitVisualizer::regenerateVertices(const glm::vec3& basePos)
{
- const Orbit& orbit = _map.getOrbit(_particleId);
+ Orbit orbit(_map.getOrbit(_particleId));
_vertices.clear();
for (int i = 0; i < _vertexCount; i++)
@@ -36,7 +36,9 @@ void OrbitVisualizer::regenerateVertices(const glm::vec3& basePos)
// better to actually create a first-class ellipse object and use that to generate
// a nice continuous mesh, instead of using orbital positions.
float t = (float)i / (float)_vertexCount * 2.0 * _pi;
- glm::vec3 pos = orbit.getPositionFromMeanAnomaly(t, 1);
+ orbit.setTrueAnomaly(t);
+ //glm::vec3 pos = orbit.getPositionFromMeanAnomaly(t, 1);
+ glm::vec3 pos = orbit.getPosition(1);
pos += basePos;
// Vertices come out of the library with X and Y being in the 'flat' plane. Re-order them
diff --git a/src/particlevisualizer.cpp b/src/particlevisualizer.cpp
index 82fa1bf..3f1fd22 100644
--- a/src/particlevisualizer.cpp
+++ b/src/particlevisualizer.cpp
@@ -10,8 +10,8 @@ ParticleVisualizer::ParticleVisualizer(const ParticleMap& map, const std::string
void ParticleVisualizer::render(float time)
{
// TODO: get mean anomly from particle which has the mass!!
- const float meanAnomaly = time;
- glm::vec3 pos = _map.getParticlePosition(_particleId, meanAnomaly);
+ //const float meanAnomaly = time;
+ glm::vec3 pos = _map.getParticlePosition(_particleId);
float y = pos.z;
pos.z = pos.y;
pos.y = y;