# 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:

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

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 on other sites

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

##### 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 on other sites

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

##### 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 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 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 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 on other sites

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

## 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

• 17
• 15
• 10
• 10
• 12