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.

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.

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. :(

I finally found my old code for RK2 and RK4 -- it works way better with RK4 as the integrator.

https://www.gamedev.net/forums/topic/626392-numerical-integration-methods/

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.

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.

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!

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

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

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