Jump to content
  • Advertisement
cowcow

Weird circular orbit problem

Recommended Posts

Posted (edited)

I am having a strange circular orbit problem. I set up the mass of the Sun, plus Mercury's initial orbit position and velocity:

https://github.com/sjhalayka/mercury_orbit_glut/blob/a0213a84badf238618c19bc7eb08d1754a84cb9b/main.h#L50

I use these parameters to calculate Mercury's path through space, via Euler integration: 

https://github.com/sjhalayka/mercury_orbit_glut/blob/a0213a84badf238618c19bc7eb08d1754a84cb9b/main.cpp#L24

The problem is, if I set Mercury's initial velocity to be sqrt(G*sun_mass/distance), then there is still some eccentricity -- the distance and velocity vector length aren't constant, like you'd expect for a perfectly circular orbit. And the discrepancy is not small, the orbit velocity length varies by 100s of metres per second, so it's not a precision issue. Anyone have any idea why my code produces unwanted results? Thanks for your time and attention.

 

 

Edited by cowcow

Share this post


Link to post
Share on other sites
Advertisement

How did you derive sqrt(G*sun_mass/distance)? (I assume "sum_mass" is a typo.)

Try sqrt(0.5*G*sun_mass/distance) instead.

 

Share this post


Link to post
Share on other sites

Yes. it was a typo.

I got the formula from many sites, including this one: http://orbitsimulator.com/formulas/vcirc.html

Trying the factor of 0.5, I get a very eccentric orbit. :(

Share this post


Link to post
Share on other sites

You can with a very small time step, to see if the orbit is much more circular in that case.

Ideally you would use a symplectic integrator so your planets don't gain or lose energy.

Share this post


Link to post
Share on other sites
Posted (edited)

Can you please provide some code for a symplectic integrator? I need the position vector length and velocity vector length to be constant over time, and RK4 seems to do the trick for the most part, well, much more than Euler integration does anyway. Not sure how much better a symplectic integrator could be.

Edited by cowcow

Share this post


Link to post
Share on other sites

Here's some code that uses a 4-th degree symplectic integrator:

#include <iostream>
#include <cmath>

struct Vector3 {
  double x, y, z;
};

Vector3 operator+(Vector3 const &v1, Vector3 const &v2) {
  return Vector3{v1.x + v2.x, v1.y + v2.y, v1.z + v2.z};
}

Vector3 &operator+=(Vector3 &v1, Vector3 const &v2) {
  v1.x += v2.x;
  v1.y += v2.y;
  v1.z += v2.z;

  return v1;
}

Vector3 operator-(Vector3 const &v1, Vector3 const &v2) {
  return Vector3{v1.x - v2.x, v1.y - v2.y, v1.z - v2.z};
}

Vector3 operator*(double scalar, Vector3 const &v) {
  return Vector3{scalar * v.x, scalar * v.y, scalar * v.z};
}

Vector3 operator*(Vector3 const &v, double scalar) {
  return scalar * v;
}

double length(Vector3 const &v) {
  return std::sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
}

std::ostream &operator<<(std::ostream &os, Vector3 const &v) {
  return os << '(' << v.x << ',' << v.y << ',' << v.z << ')';
}

// Taken from Wikipedia's example of a 4th-order symplectic integrator
void SI_order4(Vector3 &x, Vector3 &v, Vector3 (*Acceleration)(Vector3 const &), double delta_t) {
  double const cr2 = std::pow(2.0, 1.0/3.0);
  double const c[4] = {
    1.0/(2.0*(2.0-cr2)),
    (1.0-cr2)/(2.0*(2.0-cr2)),
    (1.0-cr2)/(2.0*(2.0-cr2)),
    1.0/(2.0*(2.0-cr2))
  };
  double const d[4] = {
    1.0/(2.0-cr2),
    -cr2/(2.0-cr2),
    1.0/(2.0-cr2),
    0.0
  };
  
  for (int s = 0; s < 4; ++s) {
    x += c[s] * delta_t * v;
    v += d[s] * delta_t * Acceleration(x);
  }
}

double const gravitational_constant = 1.0;
double const sun_mass = 1.0;
Vector3 const mercury_initial_position{0.0, 1.0, 0.0};
Vector3 const mercury_initial_velocity{1.0, 0.0, 0.0};
double const delta_t = 0.02;

Vector3 gravitational_acceleration(Vector3 const &x) {
  return - gravitational_constant * sun_mass / std::pow(length(x), 3) * x;
}

int main() {
  Vector3 x = mercury_initial_position;
  Vector3 v = mercury_initial_velocity;
  
  for (int i = 0; i < 158; ++i) {
    std::cout << "x=" << x << " v=" << v << " (distance=" << length(x) << ", speed=" << length(v) << ")\n";
    SI_order4(x, v, gravitational_acceleration, delta_t);
  }

}

 

 

Enjoy!

 

 

Share this post


Link to post
Share on other sites

Holy cow. I'll implement that in my code sometime soon. Thank you so very much for the nice C++ code!

Share this post


Link to post
Share on other sites

For some reason I can now log in to my old account. So long cowcow.

Share this post


Link to post
Share on other sites

Mercury's orbit is the least circular of all the major planets with an eccentricity of about 0.206.

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

GameDev.net is your game development community. Create an account for your GameDev Portfolio and participate in the largest developer community in the games industry.

Sign me up!