Numerical approximation of planetary orbit

New day new Problem...

I did it the way Eric said and I created a method for the Acceleration!

It looks like this:

``````static Vector2 AccelerationOfTheMoon(Vector2 position)
{
double Mm = 7.349 * Math.Pow(10, 22);

double MagnitudeOfForce()
{
double MagF = position.LengthSquared();
return MagF;
}

Vector2 ForceAsVector()
{
return position.NegativeUnitVector() * MagnitudeOfForce();
}

Vector2 acceleration = ForceAsVector() / Mm;

}
``````

The way i see it, my method which is called AccelerationOfTheMoon receives a Vector 2 with the position. Now i want to work with this vector so i created another Method which should return my Magnitude as a scalar(double). For this I created a Method to calculate the Magnitude of a Vector and squares it. When i call my MagnitudeOfForce-Method inside my ForceAsVector-Method I'm trying to multiply a negative unit vector with a scalar and return it as a Vector2.

Maybe I'm doing horrificly bad here, but I'm trying hard to understand everything Eric did and still need some of your help!

Thanks.

• I don't understand your algorithm at all; the relative positions and velocities of the moon are vectors in 3-space; you have the earth and the moon in a straight line. Can you explain your algorithm in more detail? – Eric Lippert Apr 17 at 13:17
• gravity, acceleration... all of them should be a vector, but what you write are scalars, even the location is a scalar!! How could you get a ellipse in a line? – shingo Apr 17 at 13:22
• Also, looking at your code, it appears that you've modeled gravity as a repelling force, not an attracting force! The earth is pushing the moon away in your system; of course it is not going to slow down. – Eric Lippert Apr 17 at 13:24
• What do you mean exactly by increasing the indent? When I'm editing the code it appears as indented! Furthermore i understand your confusion with the values being as scalars. I wanted to multiply them with trigonometric functions in order to get an ellipse, but I really have no idea how! Finally could you explain more detailed what you mean by the gravity acting as a repelling force? How would you change the code! Thanks in advance and sorry for the inconvenience, I'm a real beginner at programming! – Hannes Hauser Apr 17 at 13:31
• This may be a too-hard problem for a beginner. Orbital mechanics is not the easiest beginner subject! But you are well on your way to modeling spring forces; maybe try that first. – Eric Lippert Apr 17 at 13:46

You've made two obvious mistakes and a poor choice of algorithm.

First, you've modeled the position of the moon as a single number rather than a vector in 3-space, so what you're modeling here is simple harmonic motion, like a spring constrained to a single dimension, and not an orbit. Start by modeling positions and velocities as 3-space vectors, not doubles.

Second, you've made the sign of the gravitational force positive, so that it is a repelling force between two bodies, rather than attractive. The sign of the force has to be in the direction of the earth.

Third, your implementation of Euler's algorithm seems to be correct, but Euler is a poor choice for numerically solving orbital problems because it is not conservative; you can easily get into situations where the moon gains or loses a little bit of momentum, and that adds up, and wrecks your nice elliptical orbit.

Since the Moon's orbit is Hamiltonian, use a symplectic algorithm instead; it's designed to simulate conservative systems.

https://en.wikipedia.org/wiki/Symplectic_integrator

This approach, and your Eulerian approach, are fundamentally about finding the state vectors:

https://en.wikipedia.org/wiki/Orbital_state_vectors

However, for your simple 2-body system there are easier methods. If what you want to do is make a simulation like Kerbal Space Program, where only the body you are orbiting affects your orbit, and planets with multiple moons are "on rails", you don't need to simulate the system on every time unit to work out the sequence of state vectors; you can compute them directly given the Keplerian elements:

https://en.wikipedia.org/wiki/Orbital_elements

We know the six elements for the moon to a high degree of precision; from those you can compute the position in 3-space of the moon at any time, again, assuming that nothing perturbs its orbit. (In reality, the moon's orbit is changed by the sun, by the fact that the earth is not a sphere, by the tides, and so on.)

UPDATE:

First off, since this is for coursework you are required to cite all your sources, and that includes getting help from the internet. Your instructors must know what work is yours and what work you had someone else do for you.

You asked how to do this work in two dimensions; that seems wrong, but whatever, do what the course work says.

The rule that I wish more beginners were taught is: make a type, a method or a variable which solves your problem. In this case, we wish to represent the behavior of a complex value, so it should be a type, and that type should be a value type. Value types are `struct` in C#. So let's do that.

``````struct Vector2
{
double X { get; }
double Y { get; }
public Vector2(double x, double y)
{
this.X = x;
this.Y = y;
}
``````

Notice that vectors are immutable, just like numbers. You never mutate a vector. When you need a new vector, you make a new one.

What operations do we need to perform on vectors? Vector addition, scalar multiplication, and scalar division is just fancy multiplication. Let's implement those:

``````    public static Vector2 operator +(Vector2 a, Vector2 b) =>
new Vector2(a.X + b.X, a.Y + b.Y);
public static Vector2 operator -(Vector2 a, Vector2 b) =>
new Vector2(a.X - b.X, a.Y - b.Y);
public static Vector2 operator *(Vector2 a, double b) =>
new Vector2(a.X * b, a.Y * b);
public static Vector2 operator /(Vector2 a, double b) =>
a * (1.0 / b);
``````

We can do multiplication in the other order too, so let's implement that:

``````    public static Vector2 operator *(double b, Vector2 a) =>
a * b;
``````

Making a vector negative is the same as multiplying it by -1:

``````    public static Vector2 operator -(Vector2 a) => a * -1.0;
``````

And to help us debug:

``````    public override string ToString() => \$"({this.X},{this.Y})";

}
``````

And we're done with vectors.

We are trying to simulate the evolution of orbital state parameters, so again make a type. What are the state parameters? Position and velocity:

``````struct State
{
Vector2 Position { get; }
Vector2 Velocity { get; }
public State(Vector2 position, Vector2 velocity)
{
this.Position = position;
this.Velocity = velocity;
}
``````

Now we get to the core algorithm. What do we have? a state and an acceleration. What do we want? A new state. Make a method:

``````    public State Euler(Vector2 acceleration, double step)
{
Vector2 newVelocity = this.Velocity + acceleration * step;
Vector2 newPosition = this.Position + this.Velocity * step;
return new State(newPosition, newVelocity);
}
}
``````

Super. What's left? We need to work out the acceleration. Make a method:

``````static Vector2 AccelerationOfTheMoon(Vector2 position)
{
// You do this. Acceleration is force divided by mass,
// force is a vector, mass is a scalar. What is the force
// given the position? DO NOT MAKE A SIGN MISTAKE AGAIN.
}
``````

And now you have all the parts you need. Starting from an initial state you can repeatedly call AccelerationOfTheMoon to get the new acceleration, and then call Euler to get the new state, and repeat.

As a beginner, study these techniques carefully. Notice what I did: I made types to represent concepts, and methods to represent operations on those concepts. Once you do that, the program becomes extremely clear to read.

Think about how you would extend this program using these techniques; we made a hard-coded `AccelerationOfTheMoon` method, but what if we wanted to compute accelerations of several bodies? Again, make a type to represent `Body`; a body is characterized by `State` and mass. What if we wanted to solve the n-body problem? Our acceleration computation would have to take the other bodies as a parameter. And so on.

• Ok this one really helped! I tried to do this with 2D vectors but i really have no idea how to work with vectors in C#! If you could show me i would highly appreciate it. Furthermore I HAVE to work with Eulers approximation I could simply increase the accuracy by using RungeKutta Method, which i have no experience with! The only two formulas i am allowed to use are Newtons Graviational Law as well as Newton's 2nd Law of Motion! So i think the links you suggested are a little bit too complex for my simulation here, but are greatly appreciated! – Hannes Hauser Apr 17 at 13:58
• @HannesHauser: Is this for an assignment? You should have said so. It's fine to get help with your assignments, but I've wasted time trying to help solve a problem you don't have. Also, I don't see why you'd solve the problem with 2-D vectors; the moon exists in a 3-dimensional space. – Eric Lippert Apr 17 at 14:08
• Oh i see but as i created this question I already told the StackOverflow editor that I needed this for my work at UNI! But anyways I'm sorry to waste your time, that was definitely not what i wanted! Our lecturer said that we had to do it in 2D since moon and earth are almost perfectly aligned which gives us the possibility to to all calculations and visualizations in 2D. – Hannes Hauser Apr 17 at 14:15
• The moon is not almost perfectly aligned with the ecliptic! If it were then we'd have a lunar eclipse every month! – Eric Lippert Apr 17 at 14:16
• I'm really sorry to bother you with that but our task was to program it in 2D and I don't really know any better, im an absolute beginner at programming and my field of study differs completely from orbital mechanics! – Hannes Hauser Apr 17 at 14:25