From 788b03b2c636834f1c96f6002edf08d52b87b916 Mon Sep 17 00:00:00 2001 From: Cat Flynn Date: Tue, 24 Sep 2024 00:46:46 +0100 Subject: [PATCH] feat: update orbit based on current velocity wip: planet position is not conserved --- lib/skein/include/skein/orbit.h | 10 +++- lib/skein/include/skein/particlemap.h | 3 ++ lib/skein/src/orbit.cpp | 30 ++++++++--- lib/skein/src/particlemap.cpp | 73 ++++++++++++++++++++++++++- src/main.cpp | 31 ++++++++++-- src/orbitvisualizer.cpp | 2 +- 6 files changed, 135 insertions(+), 14 deletions(-) diff --git a/lib/skein/include/skein/orbit.h b/lib/skein/include/skein/orbit.h index f2e80fd..89000df 100644 --- a/lib/skein/include/skein/orbit.h +++ b/lib/skein/include/skein/orbit.h @@ -25,14 +25,20 @@ public: void setLongitudeOfAscendingNode(double longitudeOfAscendingNode); glm::dvec3 getPosition(double gravitationalParameter, double time) const; - glm::dvec3 getPositionFromMeanAnomaly(double meanAnomaly) 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; + //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 double getMeanAnomaly(const double gravitationalParameter, const double time) const; + private: Vector6 _keplerianElements; const double getEccentricAnomaly(const double meanAnomaly) const; const double getTrueAnomaly(const double meanAnomaly) const; - const Vector6 getCartesianCoordinates(const double meanAnomaly) const; }; diff --git a/lib/skein/include/skein/particlemap.h b/lib/skein/include/skein/particlemap.h index d25fab5..afbb1b5 100644 --- a/lib/skein/include/skein/particlemap.h +++ b/lib/skein/include/skein/particlemap.h @@ -15,11 +15,14 @@ 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); private: std::map _particles; diff --git a/lib/skein/src/orbit.cpp b/lib/skein/src/orbit.cpp index c1ab218..c251162 100644 --- a/lib/skein/src/orbit.cpp +++ b/lib/skein/src/orbit.cpp @@ -88,9 +88,9 @@ const double Orbit::getEccentricAnomaly(const double meanAnomaly) const return eccentricAnomaly; } -glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly) const +glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly, double gravitationalParameter) const { - Vector6 cartesian = getCartesianCoordinates(meanAnomaly); + Vector6 cartesian = getCartesianCoordinates(meanAnomaly, gravitationalParameter); return glm::dvec3( cartesian[astro::xPositionIndex], cartesian[astro::yPositionIndex], @@ -98,10 +98,28 @@ glm::dvec3 Orbit::getPositionFromMeanAnomaly(const double meanAnomaly) const } 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); + return glm::dvec3( + cartesian[astro::xVelocityIndex], + cartesian[astro::yVelocityIndex], + cartesian[astro::zVelocityIndex]); +} + +glm::dvec3 Orbit::getVelocity(double gravitationalParameter, double time) const +{ + return getVelocityFromMeanAnomaly(getMeanAnomaly(gravitationalParameter, time), gravitationalParameter); +} + +const double Orbit::getMeanAnomaly(const double gravitationalParameter, const double time) const { double meanMotion = astro::computeKeplerMeanMotion(getSemiMajorAxis(), gravitationalParameter); - double meanAnomaly = meanMotion * time; - return getPositionFromMeanAnomaly(meanAnomaly); + return meanMotion * time; } const double Orbit::getTrueAnomaly(const double meanAnomaly) const @@ -113,11 +131,11 @@ const double Orbit::getTrueAnomaly(const double meanAnomaly) const eccentricity); } -const Vector6 Orbit::getCartesianCoordinates(const double meanAnomaly) const +const Vector6 Orbit::getCartesianCoordinates(const double meanAnomaly, const double gravitationalParameter) const { Vector6 kepler(_keplerianElements); kepler[astro::trueAnomalyIndex] = getTrueAnomaly(meanAnomaly); - return astro::convertKeplerianToCartesianElements(kepler, 1.0); + return astro::convertKeplerianToCartesianElements(kepler, gravitationalParameter); } //glm::dvec3 Orbit::getTangent(const double meanAnomaly) const diff --git a/lib/skein/src/particlemap.cpp b/lib/skein/src/particlemap.cpp index 477b013..20ed187 100644 --- a/lib/skein/src/particlemap.cpp +++ b/lib/skein/src/particlemap.cpp @@ -1,5 +1,6 @@ #include "skein/particlemap.h" #include "skein/particle.h" +#include "astro/orbitalElementConversions.hpp" const Particle& ParticleMap::getParticle(const std::string& id) const { @@ -17,6 +18,25 @@ const Orbit& ParticleMap::getOrbit(const std::string& id) const return _orbits.at(id); } +glm::dvec3 ParticleMap::getParticleLocalPosition(const std::string& id, double time) const +{ + if (_orbits.find(id) == _orbits.end()) + return {0,0,0}; + + const Particle* particle = &(_particles.at(id)); + glm::dvec3 pos(0,0,0); + + //const std::string& id = particle->getId(); + const std::string& parentId = _relationships.at(id); + const Particle& parent = _particles.at(parentId); + const double gravitationalParameter = parent.getGravitationalParameter(); + const Orbit& orbit = _orbits.at(id); + + pos += orbit.getPosition(gravitationalParameter, time); + + return pos; +} + glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time) const { if (_orbits.find(id) == _orbits.end()) @@ -27,6 +47,7 @@ glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time) do { + // TODO: this is shadowing a parameter with the same name const std::string& id = particle->getId(); const std::string& parentId = _relationships.at(id); const Particle& parent = _particles.at(parentId); @@ -49,6 +70,54 @@ glm::dvec3 ParticleMap::getParticlePosition(const std::string& id, double time) return pos; } +glm::dvec3 ParticleMap::getParticleLocalVelocity(const std::string& id, double time) const +{ + if (_orbits.find(id) == _orbits.end()) + return {0,0,0}; + + const Particle* particle = &(_particles.at(id)); + glm::dvec3 vel(0,0,0); + + //const std::string& id = particle->getId(); + const std::string& parentId = _relationships.at(id); + const Particle& parent = _particles.at(parentId); + const double gravitationalParameter = parent.getGravitationalParameter(); + const Orbit& orbit = _orbits.at(id); + + vel += orbit.getVelocity(gravitationalParameter, time); + + return vel; +} + +void ParticleMap::setParticleLocalVelocity(const std::string& id, const double& time, const glm::dvec3& velocity) +{ + // we want to convert a local position and velocity into a new set of keplerian elements. + + // first, get the local position of the e + const Orbit& orbit = _orbits.at(id); + const std::string& parentId = _relationships.at(id); + const Particle& parent = _particles.at(parentId); + const double gravitationalParameter = parent.getGravitationalParameter(); + const double meanAnomaly = orbit.getMeanAnomaly(gravitationalParameter, time); + Vector6 cartesian(orbit.getCartesianCoordinates(meanAnomaly, gravitationalParameter)); + + // update velocity, leave position as it is + cartesian[astro::xVelocityIndex] = velocity.x; + cartesian[astro::yVelocityIndex] = velocity.y; + cartesian[astro::zVelocityIndex] = velocity.z; + + const Vector6 keplerian = astro::convertCartesianToKeplerianElements(cartesian, gravitationalParameter); + + Orbit newOrbit; + newOrbit.setSemiMajorAxis(keplerian[astro::semiMajorAxisIndex]); + newOrbit.setEccentricity(keplerian[astro::eccentricityIndex]); + newOrbit.setInclination(keplerian[astro::inclinationIndex]); + newOrbit.setArgumentOfPeriapsis(keplerian[astro::argumentOfPeriapsisIndex]); + newOrbit.setLongitudeOfAscendingNode(keplerian[astro::longitudeOfAscendingNodeIndex]); + + setRelationship(parentId, id, newOrbit); +} + void ParticleMap::setParticle(const Particle& particle) { _particles.insert({particle.getId(), particle}); @@ -58,6 +127,6 @@ void ParticleMap::setRelationship(const std::string& parentId, const std::string { // map children to parent - children can only have one parent, so we use the map to // identify the parent directly using the child as the key - _relationships.insert({childId, parentId}); - _orbits.insert({childId, orbit}); + _relationships.insert_or_assign(childId, parentId); + _orbits.insert_or_assign(childId, orbit); } diff --git a/src/main.cpp b/src/main.cpp index 891b0fa..990c2a9 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -35,20 +35,23 @@ struct Input { bool pauseTime; + bool printCartesian; } input; void keyCallback(GLFWwindow* window, int key, int scancode, int action, int mods) { - if (key == GLFW_KEY_SPACE && action == GLFW_PRESS) + if (action == GLFW_PRESS) { - input.pauseTime = true; + input.pauseTime = key == GLFW_KEY_SPACE; + input.printCartesian = key == GLFW_KEY_W; } } void clearInput() { input.pauseTime = false; + input.printCartesian = false; } // now should always increase linearly with real world time, this should not be modified by input @@ -69,6 +72,11 @@ void updateTime() engineTime.now = now; } +void printVector(const glm::dvec3& vector) +{ + std::cout << "(" << vector.x << ", " << vector.y << ", " << vector.z << ")"; +} + int main() { GLFWwindow* window = nullptr; @@ -150,7 +158,24 @@ int main() } else { - printf("maneouvre mode!\n"); + if (input.printCartesian) + { + + const glm::dvec3 p = map.getParticleLocalPosition("station", time); + printVector(p); + std::cout << std::endl; + + // TODO: what unit is this in lol + const glm::dvec3 v = map.getParticleLocalVelocity("station", time); + const glm::dvec3 newVelocity = v * 0.9; + + printVector(v); + std::cout << " -> "; + printVector(newVelocity); + std::cout << std::endl; + + map.setParticleLocalVelocity("station", time, newVelocity); + } } // rendering diff --git a/src/orbitvisualizer.cpp b/src/orbitvisualizer.cpp index 44a39b7..79bf4ea 100644 --- a/src/orbitvisualizer.cpp +++ b/src/orbitvisualizer.cpp @@ -36,7 +36,7 @@ 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); + glm::vec3 pos = orbit.getPositionFromMeanAnomaly(t, 1); pos += basePos; // Vertices come out of the library with X and Y being in the 'flat' plane. Re-order them