<p>There is no obvious astronomy error, but your use of math.pow is wrong</p>

<p>You write <code>constgrav = math.pow(6.67,-11)</code>.  That means <span class="math-container">$G = 6.67^{-11}$</span>.  You mean <code>constgrav = 6.67e-11</code> or <span class="math-container">$G=6.67\times 10^{-11}$</span>.</p>

<p>With completely different initial conditions, it is not surprising that the Earth flies off into space.</p>

<p>As a suggestion, try not using SI units for this, as it makes for some very small and very large intermediate values, the computers have difficulty in rounding.  Instead measure distance in Astronomical Units, Mass in "solar masses" and time in years. You can convert back to SI units when printing or plotting.</p>

CC BY-SA 4.0
2019-01-23T09:11:26.527Z
2019-01-23T09:11:26.527Z
James K
<p>I think you are using <code>math.pow(x, y)</code> incorrectly. You should write <code>constgrav = 6.67E-11</code> for example, standard scientific notation. Python's <code>math.pow</code> returns <span class="math-container">$x^y$</span> which is not the same thing. </p>

<p>Actually, <code>math.pow(x, y)</code> is similar, perhaps the same as <code>x**y</code>.</p>

<p>Also, you forgot a negative sign in the acceleration: gravity is attractive! Changing it to a negative and fixing the scientific notation, no more hyperbola; <strong>your program works nicely, congratulations!</strong></p>

<p>Add more digits and it should close nicely into a circle.</p>

<pre><code>def aceleracao_equacao(pos_1, pos_2):   
    r_cubo = (pos_1**2 + pos_2**2)**1.5
    k = constgrav*massasol
    aceleracao = -(k*pos_1)/r_cubo        # ¡negative sign!
    return aceleracao
</code></pre>

<p>also</p>

<pre><code>constgrav     = 6.67E-11                #constante gravitacional
massasol      = 1.98E+30                #massa do sol
massaterra    = 5.97E+24                #massa da terra
xo            = 1.5E+11                 #posição inicial no eixo x
</code></pre>

<p>and </p>

<pre><code>t             = 3.1536E+07              #tempo total = tempo de um ano
</code></pre>

<p><a href="https://i.sstatic.net/0pfFz.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/0pfFz.png" alt="enter image description here"></a> </p>

CC BY-SA 4.0
2019-01-23T20:27:15.537Z
Another question... Trying using vyo = sqrt((3*constgrav*massasol)/(2*xo)) it supposed to be an ellipse, but it's a parabola or hyperbole again. What I supposed to do?
Anna Paula Mendes
CC BY-SA 4.0
2019-01-23T22:24:15.733Z
@AnnaPaulaMendes it works for me. The higher velocity means it will orbit farther from the Sun on average, so it will take longer. Try `t  = 3.*3.1536E+7` and add `plt.plot([0], [0], 'or', ms=12)` to `graph()` for a red `o` for the Sun: https://i.sstatic.net/liJ4o.png
uhoh
CC BY-SA 4.0
2019-01-23T10:17:39.317Z
2019-01-23T10:34:56.773Z
2019-01-23T10:34:56.773Z
uhoh
<p>I'm trying to simulate Earth's orbit using 2nd Newton's Law and cartesian coordinates. The program seems ok, but when plotting graph, it's a parabola.</p>

<p>I know the shape of orbit changes according velocity. I'm using vy0 = sqrt(G*M/R) and it supposed to be an circle.</p>

<p>(Sorry the bad english)</p>

<p>Here's the code:</p>

<pre><code>'-------------------Bibliotecas Importadas-------------------'

import math
import matplotlib.pyplot as plt

'-------------------Listas Criadas-------------------'

pos_x = []
pos_y = []
vel_x = []
vel_y = []
teo_virial = []
tempo = []

'-------------------------- Definição de Funções -----------------------------'

def velocidade_atualizador(velocidade_inicial, aceleracao):
    velocidade = velocidade_inicial + (aceleracao * dt)
    return velocidade

def tabela(nome, lista): 
    arquivo = open(nome + '.txt', 'a')
    for i in range(len(lista)):
        arquivo.write(str(lista[i]) + '\n')
    arquivo.close

def graph(x, y, eixo_x, eixo_y): 
    plt.plot(x, y)
    plt.xlabel(eixo_x)
    plt.ylabel(eixo_y)
    plt.show()

def aceleracao_equacao(pos_1, pos_2):   
    r_cubo = (pos_1**2 + pos_2**2)**1.5
    k = constgrav*massasol
    aceleracao = (k*pos_1)/r_cubo
    return aceleracao

def posicao_equacao(posicao_inicial, velocidade_inicial, aceleracao):
    pos = posicao_inicial + (velocidade_inicial*dt) + ((aceleracao*(dt**2))/2)
    return pos

'-------------------Condições Iniciais-------------------'

constgrav = math.pow(6.67, -11)             #constante gravitacional
massasol = math.pow(1.98, 30)               #massa do sol
massaterra = math.pow(5.97, 24)             #massa da terra
xo = math.pow(1.5, 11)                      #posição inicial no eixo x
yo = 0                                    #posição inicial no eixo y
vxo = 0                                   #velocidade inicial no eixo x
vyo = math.sqrt((constgrav*massasol)/xo)                               #velocidade inicial no eixo y
npasso = 10000                                #numero de passos
t = math.pow(3.1536, 7)                               #tempo total = tempo de um ano
conttempo = 0                             #contador de tempo
dt = t/npasso                               #delta t 

'-------------------Laço de Repetição-------------------'

while (t &gt; conttempo):

    ax = aceleracao_equacao(xo, yo)
    ay = aceleracao_equacao(yo, xo)
    a_modulo = math.sqrt(ax**2 + ay**2)

    vx = velocidade_atualizador(vxo, ax)
    vy = velocidade_atualizador(vyo, ay)
    v_modulo = math.sqrt(vx**2 + vy**2)

    x = posicao_equacao(xo, vxo, ax)
    y = posicao_equacao(yo, vyo, ay)
    r = math.sqrt(x**2 + y**2)

    energia_pot = (- constgrav * massasol * massaterra)/math.sqrt(x**2 + y**2)
    energia_cin = (massaterra * (v_modulo**2))/2
    virial = energia_pot / energia_cin    

    print('x =', x)
    print('y =', y)
    print('virial =', virial)

    pos_x.append(x)
    pos_y.append(y)
    vel_x.append(vx)
    vel_y.append(vy)
    teo_virial.append(virial)
    tempo.append(conttempo)

    xo = x
    yo = y
    vxo = vx
    vyo = vy

    conttempo = conttempo + dt

tabela('Lista x', pos_x)
tabela('Lista y', pos_y)
tabela('Lista vx', vel_x)
tabela('Lista vy', vel_y)
tabela('Lista tempo', tempo)

'-------------------Gráfico-------------------'

graph(pos_x, pos_y, 'x (m)', 'y (m)')
graph(vel_x, tempo, 'vx (m/s)', 't (s)')
graph(vel_y, tempo, 'vy (m/s)', 't (s)')
</code></pre>

<hr>

<p>Thank you very much, guys!! I was thinking but don't know what to do anymore. </p>

<p>After changing the code, it works. I can rest in peace now. Hahaha.</p>

CC BY-SA 4.0
2019-01-23T08:43:00.483Z
I'm voting to close this question as off-topic because Astronomy SE is not for debugging code.
StephenG - Help Ukraine
CC BY-SA 4.0
2019-01-23T19:19:44.837Z
BTW, if you want to make your code more readable for an international audience you should use English variable names, rather than Portuguese.
PM 2Ring
CC BY-SA 4.0
2019-01-23T19:26:43.240Z
Your algorithm is ok for orbits with low eccentricity, but it will quickly accumulate errors & need a very small step time for more eccentric orbits. I recommend a symplectic integrator like synchronized [Leapfrog](https://en.wikipedia.org/wiki/Leapfrog_integration), or velocity Verlet. You can do that with only small changes to your code. Also, rather than using mass and G you should use the [standard gravitational parameter](https://en.wikipedia.org/wiki/Standard_gravitational_parameter), which is much more precise.
PM 2Ring
CC BY-SA 4.0
2019-01-24T03:22:02.493Z
I suggest you try improved Euler (https://en.wikipedia.org/wiki/Heun%27s_method).  It improves the accuracy a lot.
irchans
<orbit><python>
Earth Orbit Simulator
Anna Paula Mendes
# Earth Orbit Simulator

(Asked by: username_0 on 2019-01-23)

I'm trying to simulate Earth's orbit using 2nd Newton's Law and cartesian coordinates. The program seems ok, but when plotting graph, it's a parabola.

I know the shape of orbit changes according velocity. I'm using vy0 = sqrt(G\*M/R) and it supposed to be an circle.

(Sorry the bad english)

Here's the code:

```
'-------------------Bibliotecas Importadas-------------------'

import math
import matplotlib.pyplot as plt

'-------------------Listas Criadas-------------------'

pos_x = []
pos_y = []
vel_x = []
vel_y = []
teo_virial = []
tempo = []

'-------------------------- Definição de Funções -----------------------------'

def velocidade_atualizador(velocidade_inicial, aceleracao):
    velocidade = velocidade_inicial + (aceleracao * dt)
    return velocidade

def tabela(nome, lista): 
    arquivo = open(nome + '.txt', 'a')
    for i in range(len(lista)):
        arquivo.write(str(lista[i]) + '\n')
    arquivo.close

def graph(x, y, eixo_x, eixo_y): 
    plt.plot(x, y)
    plt.xlabel(eixo_x)
    plt.ylabel(eixo_y)
    plt.show()

def aceleracao_equacao(pos_1, pos_2):   
    r_cubo = (pos_1**2 + pos_2**2)**1.5
    k = constgrav*massasol
    aceleracao = (k*pos_1)/r_cubo
    return aceleracao

def posicao_equacao(posicao_inicial, velocidade_inicial, aceleracao):
    pos = posicao_inicial + (velocidade_inicial*dt) + ((aceleracao*(dt**2))/2)
    return pos

'-------------------Condições Iniciais-------------------'

constgrav = math.pow(6.67, -11)             #constante gravitacional
massasol = math.pow(1.98, 30)               #massa do sol
massaterra = math.pow(5.97, 24)             #massa da terra
xo = math.pow(1.5, 11)                      #posição inicial no eixo x
yo = 0                                    #posição inicial no eixo y
vxo = 0                                   #velocidade inicial no eixo x
vyo = math.sqrt((constgrav*massasol)/xo)                               #velocidade inicial no eixo y
npasso = 10000                                #numero de passos
t = math.pow(3.1536, 7)                               #tempo total = tempo de um ano
conttempo = 0                             #contador de tempo
dt = t/npasso                               #delta t 

'-------------------Laço de Repetição-------------------'

while (t > conttempo):

    ax = aceleracao_equacao(xo, yo)
    ay = aceleracao_equacao(yo, xo)
    a_modulo = math.sqrt(ax**2 + ay**2)

    vx = velocidade_atualizador(vxo, ax)
    vy = velocidade_atualizador(vyo, ay)
    v_modulo = math.sqrt(vx**2 + vy**2)

    x = posicao_equacao(xo, vxo, ax)
    y = posicao_equacao(yo, vyo, ay)
    r = math.sqrt(x**2 + y**2)

    energia_pot = (- constgrav * massasol * massaterra)/math.sqrt(x**2 + y**2)
    energia_cin = (massaterra * (v_modulo**2))/2
    virial = energia_pot / energia_cin    

    print('x =', x)
    print('y =', y)
    print('virial =', virial)

    pos_x.append(x)
    pos_y.append(y)
    vel_x.append(vx)
    vel_y.append(vy)
    teo_virial.append(virial)
    tempo.append(conttempo)

    xo = x
    yo = y
    vxo = vx
    vyo = vy

    conttempo = conttempo + dt

tabela('Lista x', pos_x)
tabela('Lista y', pos_y)
tabela('Lista vx', vel_x)
tabela('Lista vy', vel_y)
tabela('Lista tempo', tempo)

'-------------------Gráfico-------------------'

graph(pos_x, pos_y, 'x (m)', 'y (m)')
graph(vel_x, tempo, 'vx (m/s)', 't (s)')
graph(vel_y, tempo, 'vy (m/s)', 't (s)')

```

---

Thank you very much, guys!! I was thinking but don't know what to do anymore.

After changing the code, it works. I can rest in peace now. Hahaha.

**Tags:** orbit, python

> **Comments:**
> **username_1**: I'm voting to close this question as off-topic because Astronomy SE is not for debugging code.  
> **username_2**: BTW, if you want to make your code more readable for an international audience you should use English variable names, rather than Portuguese.  
> **username_2**: Your algorithm is ok for orbits with low eccentricity, but it will quickly accumulate errors & need a very small step time for more eccentric orbits. I recommend a symplectic integrator like synchronized [Leapfrog](https://en.wikipedia.org/wiki/Leapfrog_integration), or velocity Verlet. You can do that with only small changes to your code. Also, rather than using mass and G you should use the [standard gravitational parameter](https://en.wikipedia.org/wiki/Standard_gravitational_parameter), which is much more precise.  
> **username_3**: I suggest you try improved Euler (https://en.wikipedia.org/wiki/Heun%27s_method).  It improves the accuracy a lot.  

---

## 2 answers

### Answer by username_4 on 2019-01-23. Score: 5

There is no obvious astronomy error, but your use of math.pow is wrong

You write `constgrav = math.pow(6.67,-11)`. That means $G = 6.67^{-11}$. You mean `constgrav = 6.67e-11` or $G=6.67\times 10^{-11}$.

With completely different initial conditions, it is not surprising that the Earth flies off into space.

As a suggestion, try not using SI units for this, as it makes for some very small and very large intermediate values, the computers have difficulty in rounding. Instead measure distance in Astronomical Units, Mass in "solar masses" and time in years. You can convert back to SI units when printing or plotting.


---

### Answer by username_5 on 2019-01-23. Score: 1

I think you are using `math.pow(x, y)` incorrectly. You should write `constgrav = 6.67E-11` for example, standard scientific notation. Python's `math.pow` returns $x^y$ which is not the same thing.

Actually, `math.pow(x, y)` is similar, perhaps the same as `x**y`.

Also, you forgot a negative sign in the acceleration: gravity is attractive! Changing it to a negative and fixing the scientific notation, no more hyperbola; **your program works nicely, congratulations!**

Add more digits and it should close nicely into a circle.

```
def aceleracao_equacao(pos_1, pos_2):   
    r_cubo = (pos_1**2 + pos_2**2)**1.5
    k = constgrav*massasol
    aceleracao = -(k*pos_1)/r_cubo        # ¡negative sign!
    return aceleracao

```

also

```
constgrav     = 6.67E-11                #constante gravitacional
massasol      = 1.98E+30                #massa do sol
massaterra    = 5.97E+24                #massa da terra
xo            = 1.5E+11                 #posição inicial no eixo x

```

and

```
t             = 3.1536E+07              #tempo total = tempo de um ano

```

[![enter image description here](https://i.sstatic.net/0pfFz.png)](https://i.sstatic.net/0pfFz.png)

> **Comments:**
> **username_0**: Another question... Trying using vyo = sqrt((3*constgrav*massasol)/(2*xo)) it supposed to be an ellipse, but it's a parabola or hyperbole again. What I supposed to do?  
> **username_5**: @username_0 it works for me. The higher velocity means it will orbit farther from the Sun on average, so it will take longer. Try `t  = 3.*3.1536E+7` and add `plt.plot([0], [0], 'or', ms=12)` to `graph()` for a red `o` for the Sun: https://i.sstatic.net/liJ4o.png  

---

