<p>Let the two bodies involved have masses <span class="math-container">$m_1, m_2$</span>.  Start with Newton's second law <span class="math-container">$$F = ma$$</span> where <span class="math-container">$a$</span> is acceleration.  The gravitational force on body 2 from body 1 is given by</p>
<p><span class="math-container">$$F_{21} = \frac{G m_1 m_2}{|r_{21}|^3}r_{21}$$</span></p>
<p>where <span class="math-container">$r_{21}$</span> is the relative position vector for the two bodies in question.  The force on body 1 from body two is of course <span class="math-container">$F_{12}=-F_{21}$</span> since <span class="math-container">$r_{12}=-r_{21}$</span>.  If we denote the positions by <span class="math-container">$(x_1,y_1)$</span> and <span class="math-container">$(x_2,y_2)$</span>, then</p>
<p><span class="math-container">$$r_{21} = \begin{pmatrix} x_1-x_2 \\ y_1-y_2 \end{pmatrix}.$$</span></p>
<p>and</p>
<p><span class="math-container">$$|r| = \sqrt{(x_1-x_2)^2 + (y_1-y_2)^2}.$$</span>
Then Newton's law <span class="math-container">$a=F/m$</span> becomes</p>
<p><span class="math-container">\begin{align}
x_1''(t) &amp; = \frac{G m_2 (x_2-x_1)}{|r|^3} \\
y_1''(t) &amp; = \frac{G m_2 (y_2-y_1)}{|r|^3} \\
x_2''(t) &amp; = \frac{G m_1 (x_1-x_2)}{|r|^3} \\
y_2''(t) &amp; = \frac{G m_1 (y_1-y_2)}{|r|^3}.
\end{align}</span></p>
<p>Together with the initial positions and velocities, this system of ordinary differential equations (ODEs) comprises an initial value problem.  The usual approach is to write this as a first-order system of 8 equations and apply a Runge-Kutta or multistep method to solve them.</p>
<p>If you apply something simple like forward Euler or backward Euler, you will see the Earth spiral out to infinity or in toward the sun, respectively, but that is an effect of the numerical errors.  If you use a more accurate method, like the classical 4th-order Runge-Kutta method, you'll find that it stays close to a true orbit for a while but still eventually spirals inward (as long as the time step is small enough) or goes to infinity (if the time step is too large).  The right approach is to use a symplectic method, which will keep the Earth in the correct orbit -- although its phase will still be off due to numerical errors.</p>
<p>For the 2-body problem it's possible to derive a simpler system by basing your coordinate system around the center of mass.  But I think the formulation above is clearer if this is new to you.</p>

CC BY-SA 3.0
2014-10-30T20:31:51.723Z
This will take some time to digest.
statespace
CC BY-SA 3.0
2014-11-10T20:17:23.203Z
Still digesting. Too many unknown words for me but somehow I feel that at some point I will get there. For now my own algorithm is good enough for things to simply work. But when I will plug in simultaneous motion - I will be forced to get to the literature and read about proper algorithms. Given that limitations of modern hardware are much looser, I can afford to fool around with simple equations. Afraid not for long though.
statespace
CC BY-SA 3.0
2014-11-21T10:43:02.337Z
Indeed the symplectic methods are by far the most accurate but I think it is difficult for someone with no background in science to implement them. Instead you can use the very simple Euler method along with the Feynman correction. I don't think that you need something more complex than that for self-education purposes.
Chris Papageorgakis
CC BY-SA 4.0
2014-10-29T17:55:34.067Z
2020-09-28T08:59:24.413Z
2020-09-28T08:59:24.413Z
David Ketcheson
<p>Just a little disclaimer beforehand: I have never studied astronomy or any exact sciences for that matter (not even IT), so I am trying to fill this gap by self-education. Astronomy is one of the areas that has captured my attention and my idea of self-education is head on applied approach. So, straight to the point - this is orbital simulation model that I am casually working on when I have time/mood. My major goal is to create complete solar system in motion and ability to plan spacecraft launches to other planets.</p>

<p>You are all free to pick this project up at any point and have fun experimenting!</p>

<p><strong>update!!! (Nov10)</strong></p>

<ul>
<li>velocity is now proper deltaV and giving additional motion now calculates
sum vector of velocity </li>
<li>you can place as many static objects as you
like, on every time unit object in motion checks for gravity vectors
from all sources (and checks for collision)</li>
<li>greatly improved the performance of calculations</li>
<li>a fix to account for interactive mod in matplotlib. Looks like that this is default option only for ipython. Regular python3 requires that statement explicitly. </li>
</ul>

<p>Basically it is now possible to "launch" a spacecraft from the surface of the Earth and plot a mission to the Moon by making deltaV vector corrections via giveMotion(). Next in line is trying to implement global time variable to enable simultaneous motion e.g. Moon orbits Earth while spacecraft tries out a gravity assist maneuver.</p>

<p>Comments and suggestions for improvements are always welcome!  </p>

<p>Done in Python3 with matplotlib library</p>

<pre><code>import matplotlib.pyplot as plt
import math
plt.ion()

G = 6.673e-11  # gravity constant
gridArea = [0, 200, 0, 200]  # margins of the coordinate grid
gridScale = 1000000  # 1 unit of grid equals 1000000m or 1000km

plt.clf()  # clear plot area
plt.axis(gridArea)  # create new coordinate grid
plt.grid(b="on")  # place grid

class Object:
    _instances = []
    def __init__(self, name, position, radius, mass):
        self.name = name
        self.position = position
        self.radius = radius  # in grid values
        self.mass = mass
        self.placeObject()
        self.velocity = 0
        Object._instances.append(self)

    def placeObject(self):
        drawObject = plt.Circle(self.position, radius=self.radius, fill=False, color="black")
        plt.gca().add_patch(drawObject)
        plt.show()

    def giveMotion(self, deltaV, motionDirection, time):
        if self.velocity != 0:
            x_comp = math.sin(math.radians(self.motionDirection))*self.velocity
            y_comp = math.cos(math.radians(self.motionDirection))*self.velocity
            x_comp += math.sin(math.radians(motionDirection))*deltaV
            y_comp += math.cos(math.radians(motionDirection))*deltaV
            self.velocity = math.sqrt((x_comp**2)+(y_comp**2))

            if x_comp &gt; 0 and y_comp &gt; 0:  # calculate degrees depending on the coordinate quadrant
                self.motionDirection = math.degrees(math.asin(abs(x_comp)/self.velocity))  # update motion direction
            elif x_comp &gt; 0 and y_comp &lt; 0:
                self.motionDirection = math.degrees(math.asin(abs(y_comp)/self.velocity)) + 90
            elif x_comp &lt; 0 and y_comp &lt; 0:
                self.motionDirection = math.degrees(math.asin(abs(x_comp)/self.velocity)) + 180
            else:
                self.motionDirection = math.degrees(math.asin(abs(y_comp)/self.velocity)) + 270

        else:
            self.velocity = self.velocity + deltaV  # in m/s
            self.motionDirection = motionDirection  # degrees
        self.time = time  # in seconds
        self.vectorUpdate()

    def vectorUpdate(self):
        self.placeObject()
        data = []

        for t in range(self.time):
            motionForce = self.mass * self.velocity  # F = m * v
            x_net = 0
            y_net = 0
            for x in [y for y in Object._instances if y is not self]:
                distance = math.sqrt(((self.position[0]-x.position[0])**2) +
                             (self.position[1]-x.position[1])**2)
                gravityForce = G*(self.mass * x.mass)/((distance*gridScale)**2)

                x_pos = self.position[0] - x.position[0]
                y_pos = self.position[1] - x.position[1]

                if x_pos &lt;= 0 and y_pos &gt; 0:  # calculate degrees depending on the coordinate quadrant
                    gravityDirection = math.degrees(math.asin(abs(y_pos)/distance))+90

                elif x_pos &gt; 0 and y_pos &gt;= 0:
                    gravityDirection = math.degrees(math.asin(abs(x_pos)/distance))+180

                elif x_pos &gt;= 0 and y_pos &lt; 0:
                    gravityDirection = math.degrees(math.asin(abs(y_pos)/distance))+270

                else:
                    gravityDirection = math.degrees(math.asin(abs(x_pos)/distance))

                x_gF = gravityForce * math.sin(math.radians(gravityDirection))  # x component of vector
                y_gF = gravityForce * math.cos(math.radians(gravityDirection))  # y component of vector

                x_net += x_gF
                y_net += y_gF

            x_mF = motionForce * math.sin(math.radians(self.motionDirection))
            y_mF = motionForce * math.cos(math.radians(self.motionDirection))
            x_net += x_mF
            y_net += y_mF
            netForce = math.sqrt((x_net**2)+(y_net**2))

            if x_net &gt; 0 and y_net &gt; 0:  # calculate degrees depending on the coordinate quadrant
                self.motionDirection = math.degrees(math.asin(abs(x_net)/netForce))  # update motion direction
            elif x_net &gt; 0 and y_net &lt; 0:
                self.motionDirection = math.degrees(math.asin(abs(y_net)/netForce)) + 90
            elif x_net &lt; 0 and y_net &lt; 0:
                self.motionDirection = math.degrees(math.asin(abs(x_net)/netForce)) + 180
            else:
                self.motionDirection = math.degrees(math.asin(abs(y_net)/netForce)) + 270

            self.velocity = netForce/self.mass  # update velocity
            traveled = self.velocity/gridScale  # grid distance traveled per 1 sec
            self.position = (self.position[0] + math.sin(math.radians(self.motionDirection))*traveled,
                             self.position[1] + math.cos(math.radians(self.motionDirection))*traveled)  # update pos
            data.append([self.position[0], self.position[1]])

            collision = 0
            for x in [y for y in Object._instances if y is not self]:
                if (self.position[0] - x.position[0])**2 + (self.position[1] - x.position[1])**2 &lt;= x.radius**2:
                    collision = 1
                    break
            if collision != 0:
                print("Collision!")
                break

        plt.plot([x[0] for x in data], [x[1] for x in data])

Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)
Moon = Object(name="Moon", position=(100.0, 100.0), radius=1.737, mass = 7.347e22)  # position not to real scale
Craft = Object(name="SpaceCraft", position=(49.0, 40.0), radius=1, mass=1.0e4)

Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)
Craft.giveMotion(deltaV=2000.0, motionDirection=90, time=60000)
plt.show(block=True)
</code></pre>

<p><strong>How it works</strong></p>

<p>It all boils down to two things:</p>

<ol>
<li>Creating object like <code>Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)</code> with parameters of position on grid (1 unit of grid is 1000km by default but this can be changed too), radius in grid units and mass in kg.</li>
<li>Giving object some deltaV such as <code>Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)</code> obviously it requires <code>Craft = Object(...)</code> to be created in the first place as mentioned in previous point. Parameters here are <code>deltaV</code> in m/s (note that for now acceleration is instantaneous), <code>motionDirection</code> is direction of deltaV in degrees (from current position imagine 360 degree circle around object, so direction is a point on that circle) and finally parameter <code>time</code> is how many seconds after the deltaV push trajectory of the object will be monitored. Subsequent <code>giveMotion()</code> start off from last position of previous <code>giveMotion()</code>. </li>
</ol>

<p>Questions:</p>

<ol>
<li>Is this a valid algorithm to calculate orbits?</li>
<li>What are the obvious improvements to be made?</li>
<li>I have been considering "timeScale" variable that will optimize calculations, as it might not be necessary to recalculate vectors and positions for every second. Any thoughts on how it should be implemented or is it generally a good idea? (loss of accuracy vs improved performance)</li>
</ol>

<p>Basically my aim is to start a discussion on the topic and see where it leads. And, if possible, learn (or even better - teach) something new and interesting.</p>

<p><strong>Feel free to experiment!</strong></p>

<p>Try using:</p>

<pre><code>Earth = Object(name="Earth", position=(50.0, 100.0), radius=6.371, mass=5.972e24)
Moon = Object(name="Moon", position=(434.0, 100.0), radius=1.737, mass = 7.347e22)
Craft = Object(name="SpaceCraft", position=(43.0, 100.0), radius=1, mass=1.0e4)

Craft.giveMotion(deltaV=10575.0, motionDirection=180, time=322000)
Craft.giveMotion(deltaV=400.0, motionDirection=180, time=50000)
</code></pre>

<p>With two burns - one prograde at Earth orbit and one retrograde at Moon orbit I achieved stable Moon orbit. Are these close to theoreticaly expected values?</p>

<p>Suggested exercise: Try it in 3 burns - stable Earth orbit from Earth surface, prograde burn to reach Moon, retrograde burn to stabilize orbit around Moon. Then try to minimize deltaV. </p>

<p>Note: I plan to update the code with extensive comments for those not familiar with python3 syntax. </p>

CC BY-SA 3.0
2014-10-28T21:06:35.823Z
Sure, I guess. I will make more extensive comments in the code for those that want to pick it up and summarize the general logic in the question itself.
statespace
CC BY-SA 3.0
2014-10-29T05:50:46.973Z
You can give motion to any object, including Earth. For test purposes I included only object -> Earth relation in main loop. It can easily be converted that every object relates to all other objects that are created. And every object can have their own motion vector. Reason why I didn't do it - very slow calculations even for 1 object. I hope that scaling time units should help alot but I am still not sure how to do it right.
statespace
CC BY-SA 3.0
2014-10-30T20:36:27.793Z
Awesome, thanks. That link alone was worth it to start the discussion. So, first step - I will try to re-write the code so that objects check all other objects for gravity force vector. By the way, most pressing issue of time scaling is still here. And it would be interesting to check how results improve or get worse by changing time scale parameter.
statespace
CC BY-SA 3.0
2014-10-30T20:43:05.753Z
Oh, and F=m*v is really that. Assumption is that by giving motion to the object first of all it goes to absolute halt and by that velocity equals to acceleration. Like instant speed from 0 -> whatever value, and goes on indefinitely with that speed if no other forces act upon the object.
statespace
<orbit><gravity><n-body-simulations>
Exercise: 2D orbital mechanics simulation (python)
statespace
# Exercise: 2D orbital mechanics simulation (python)

(Asked by: username_0 on 2014-10-28)

Just a little disclaimer beforehand: I have never studied astronomy or any exact sciences for that matter (not even IT), so I am trying to fill this gap by self-education. Astronomy is one of the areas that has captured my attention and my idea of self-education is head on applied approach. So, straight to the point - this is orbital simulation model that I am casually working on when I have time/mood. My major goal is to create complete solar system in motion and ability to plan spacecraft launches to other planets.
You are all free to pick this project up at any point and have fun experimenting!
update!!! (Nov10)

velocity is now proper deltaV and giving additional motion now calculates
sum vector of velocity 
you can place as many static objects as you
like, on every time unit object in motion checks for gravity vectors
from all sources (and checks for collision)
greatly improved the performance of calculations
a fix to account for interactive mod in matplotlib. Looks like that this is default option only for ipython. Regular python3 requires that statement explicitly. 

Basically it is now possible to "launch" a spacecraft from the surface of the Earth and plot a mission to the Moon by making deltaV vector corrections via giveMotion(). Next in line is trying to implement global time variable to enable simultaneous motion e.g. Moon orbits Earth while spacecraft tries out a gravity assist maneuver.
Comments and suggestions for improvements are always welcome!  
Done in Python3 with matplotlib library
import matplotlib.pyplot as plt
import math
plt.ion()

G = 6.673e-11  # gravity constant
gridArea = [0, 200, 0, 200]  # margins of the coordinate grid
gridScale = 1000000  # 1 unit of grid equals 1000000m or 1000km

plt.clf()  # clear plot area
plt.axis(gridArea)  # create new coordinate grid
plt.grid(b="on")  # place grid

class Object:
    _instances = []
    def __init__(self, name, position, radius, mass):
        self.name = name
        self.position = position
        self.radius = radius  # in grid values
        self.mass = mass
        self.placeObject()
        self.velocity = 0
        Object._instances.append(self)

    def placeObject(self):
        drawObject = plt.Circle(self.position, radius=self.radius, fill=False, color="black")
        plt.gca().add_patch(drawObject)
        plt.show()

    def giveMotion(self, deltaV, motionDirection, time):
        if self.velocity != 0:
            x_comp = math.sin(math.radians(self.motionDirection))*self.velocity
            y_comp = math.cos(math.radians(self.motionDirection))*self.velocity
            x_comp += math.sin(math.radians(motionDirection))*deltaV
            y_comp += math.cos(math.radians(motionDirection))*deltaV
            self.velocity = math.sqrt((x_comp**2)+(y_comp**2))

            if x_comp > 0 and y_comp > 0:  # calculate degrees depending on the coordinate quadrant
                self.motionDirection = math.degrees(math.asin(abs(x_comp)/self.velocity))  # update motion direction
            elif x_comp > 0 and y_comp < 0:
                self.motionDirection = math.degrees(math.asin(abs(y_comp)/self.velocity)) + 90
            elif x_comp < 0 and y_comp < 0:
                self.motionDirection = math.degrees(math.asin(abs(x_comp)/self.velocity)) + 180
            else:
                self.motionDirection = math.degrees(math.asin(abs(y_comp)/self.velocity)) + 270

        else:
            self.velocity = self.velocity + deltaV  # in m/s
            self.motionDirection = motionDirection  # degrees
        self.time = time  # in seconds
        self.vectorUpdate()

    def vectorUpdate(self):
        self.placeObject()
        data = []

        for t in range(self.time):
            motionForce = self.mass * self.velocity  # F = m * v
            x_net = 0
            y_net = 0
            for x in [y for y in Object._instances if y is not self]:
                distance = math.sqrt(((self.position[0]-x.position[0])**2) +
                             (self.position[1]-x.position[1])**2)
                gravityForce = G*(self.mass * x.mass)/((distance*gridScale)**2)

                x_pos = self.position[0] - x.position[0]
                y_pos = self.position[1] - x.position[1]

                if x_pos <= 0 and y_pos > 0:  # calculate degrees depending on the coordinate quadrant
                    gravityDirection = math.degrees(math.asin(abs(y_pos)/distance))+90

                elif x_pos > 0 and y_pos >= 0:
                    gravityDirection = math.degrees(math.asin(abs(x_pos)/distance))+180

                elif x_pos >= 0 and y_pos < 0:
                    gravityDirection = math.degrees(math.asin(abs(y_pos)/distance))+270

                else:
                    gravityDirection = math.degrees(math.asin(abs(x_pos)/distance))

                x_gF = gravityForce * math.sin(math.radians(gravityDirection))  # x component of vector
                y_gF = gravityForce * math.cos(math.radians(gravityDirection))  # y component of vector

                x_net += x_gF
                y_net += y_gF

            x_mF = motionForce * math.sin(math.radians(self.motionDirection))
            y_mF = motionForce * math.cos(math.radians(self.motionDirection))
            x_net += x_mF
            y_net += y_mF
            netForce = math.sqrt((x_net**2)+(y_net**2))

            if x_net > 0 and y_net > 0:  # calculate degrees depending on the coordinate quadrant
                self.motionDirection = math.degrees(math.asin(abs(x_net)/netForce))  # update motion direction
            elif x_net > 0 and y_net < 0:
                self.motionDirection = math.degrees(math.asin(abs(y_net)/netForce)) + 90
            elif x_net < 0 and y_net < 0:
                self.motionDirection = math.degrees(math.asin(abs(x_net)/netForce)) + 180
            else:
                self.motionDirection = math.degrees(math.asin(abs(y_net)/netForce)) + 270

            self.velocity = netForce/self.mass  # update velocity
            traveled = self.velocity/gridScale  # grid distance traveled per 1 sec
            self.position = (self.position[0] + math.sin(math.radians(self.motionDirection))*traveled,
                             self.position[1] + math.cos(math.radians(self.motionDirection))*traveled)  # update pos
            data.append([self.position[0], self.position[1]])

            collision = 0
            for x in [y for y in Object._instances if y is not self]:
                if (self.position[0] - x.position[0])**2 + (self.position[1] - x.position[1])**2 <= x.radius**2:
                    collision = 1
                    break
            if collision != 0:
                print("Collision!")
                break

        plt.plot([x[0] for x in data], [x[1] for x in data])

Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24)
Moon = Object(name="Moon", position=(100.0, 100.0), radius=1.737, mass = 7.347e22)  # position not to real scale
Craft = Object(name="SpaceCraft", position=(49.0, 40.0), radius=1, mass=1.0e4)

Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000)
Craft.giveMotion(deltaV=2000.0, motionDirection=90, time=60000)
plt.show(block=True)

How it works
It all boils down to two things:

Creating object like Earth = Object(name="Earth", position=(50.0, 50.0), radius=6.371, mass=5.972e24) with parameters of position on grid (1 unit of grid is 1000km by default but this can be changed too), radius in grid units and mass in kg.
Giving object some deltaV such as Craft.giveMotion(deltaV=8500.0, motionDirection=100, time=130000) obviously it requires Craft = Object(...) to be created in the first place as mentioned in previous point. Parameters here are deltaV in m/s (note that for now acceleration is instantaneous), motionDirection is direction of deltaV in degrees (from current position imagine 360 degree circle around object, so direction is a point on that circle) and finally parameter time is how many seconds after the deltaV push trajectory of the object will be monitored. Subsequent giveMotion() start off from last position of previous giveMotion(). 

Questions:

Is this a valid algorithm to calculate orbits?
What are the obvious improvements to be made?
I have been considering "timeScale" variable that will optimize calculations, as it might not be necessary to recalculate vectors and positions for every second. Any thoughts on how it should be implemented or is it generally a good idea? (loss of accuracy vs improved performance)

Basically my aim is to start a discussion on the topic and see where it leads. And, if possible, learn (or even better - teach) something new and interesting.
Feel free to experiment!
Try using:
Earth = Object(name="Earth", position=(50.0, 100.0), radius=6.371, mass=5.972e24)
Moon = Object(name="Moon", position=(434.0, 100.0), radius=1.737, mass = 7.347e22)
Craft = Object(name="SpaceCraft", position=(43.0, 100.0), radius=1, mass=1.0e4)

Craft.giveMotion(deltaV=10575.0, motionDirection=180, time=322000)
Craft.giveMotion(deltaV=400.0, motionDirection=180, time=50000)

With two burns - one prograde at Earth orbit and one retrograde at Moon orbit I achieved stable Moon orbit. Are these close to theoreticaly expected values?
Suggested exercise: Try it in 3 burns - stable Earth orbit from Earth surface, prograde burn to reach Moon, retrograde burn to stabilize orbit around Moon. Then try to minimize deltaV. 
Note: I plan to update the code with extensive comments for those not familiar with python3 syntax.

**Tags:** orbit, gravity, n-body-simulations
> **Comments:**
> **username_0**: Sure, I guess. I will make more extensive comments in the code for those that want to pick it up and summarize the general logic in the question itself.  
> **username_0**: You can give motion to any object, including Earth. For test purposes I included only object -> Earth relation in main loop. It can easily be converted that every object relates to all other objects that are created. And every object can have their own motion vector. Reason why I didn't do it - very slow calculations even for 1 object. I hope that scaling time units should help alot but I am still not sure how to do it right.  
> **username_0**: Awesome, thanks. That link alone was worth it to start the discussion. So, first step - I will try to re-write the code so that objects check all other objects for gravity force vector. By the way, most pressing issue of time scaling is still here. And it would be interesting to check how results improve or get worse by changing time scale parameter.  
> **username_0**: Oh, and F=m*v is really that. Assumption is that by giving motion to the object first of all it goes to absolute halt and by that velocity equals to acceleration. Like instant speed from 0 -> whatever value, and goes on indefinitely with that speed if no other forces act upon the object.  

---

## 1 answers

### Answer by username_1 on 2014-10-29. Score: 11

Let the two bodies involved have masses $m\_1, m\_2$. Start with Newton's second law $$F = ma$$ where $a$ is acceleration. The gravitational force on body 2 from body 1 is given by

$$F\_{21} = \frac{G m\_1 m\_2}{|r\_{21}|^3}r\_{21}$$

where $r\_{21}$ is the relative position vector for the two bodies in question. The force on body 1 from body two is of course $F\_{12}=-F\_{21}$ since $r\_{12}=-r\_{21}$. If we denote the positions by $(x\_1,y\_1)$ and $(x\_2,y\_2)$, then

$$r\_{21} = \begin{pmatrix} x\_1-x\_2 \\ y\_1-y\_2 \end{pmatrix}.$$

and

$$|r| = \sqrt{(x\_1-x\_2)^2 + (y\_1-y\_2)^2}.$$
Then Newton's law $a=F/m$ becomes

\begin{align}
x\_1''(t) & = \frac{G m\_2 (x\_2-x\_1)}{|r|^3} \\
y\_1''(t) & = \frac{G m\_2 (y\_2-y\_1)}{|r|^3} \\
x\_2''(t) & = \frac{G m\_1 (x\_1-x\_2)}{|r|^3} \\
y\_2''(t) & = \frac{G m\_1 (y\_1-y\_2)}{|r|^3}.
\end{align}

Together with the initial positions and velocities, this system of ordinary differential equations (ODEs) comprises an initial value problem. The usual approach is to write this as a first-order system of 8 equations and apply a Runge-Kutta or multistep method to solve them.

If you apply something simple like forward Euler or backward Euler, you will see the Earth spiral out to infinity or in toward the sun, respectively, but that is an effect of the numerical errors. If you use a more accurate method, like the classical 4th-order Runge-Kutta method, you'll find that it stays close to a true orbit for a while but still eventually spirals inward (as long as the time step is small enough) or goes to infinity (if the time step is too large). The right approach is to use a symplectic method, which will keep the Earth in the correct orbit -- although its phase will still be off due to numerical errors.

For the 2-body problem it's possible to derive a simpler system by basing your coordinate system around the center of mass. But I think the formulation above is clearer if this is new to you.

> **Comments:**
> **username_0**: This will take some time to digest.  
> **username_0**: Still digesting. Too many unknown words for me but somehow I feel that at some point I will get there. For now my own algorithm is good enough for things to simply work. But when I will plug in simultaneous motion - I will be forced to get to the literature and read about proper algorithms. Given that limitations of modern hardware are much looser, I can afford to fool around with simple equations. Afraid not for long though.  
> **username_2**: Indeed the symplectic methods are by far the most accurate but I think it is difficult for someone with no background in science to implement them. Instead you can use the very simple Euler method along with the Feynman correction. I don't think that you need something more complex than that for self-education purposes.  

---

