1
$\begingroup$

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.

$\endgroup$
  • 6
    $\begingroup$ I'm voting to close this question as off-topic because Astronomy SE is not for debugging code. $\endgroup$ – StephenG Jan 23 '19 at 8:43
  • 1
    $\begingroup$ BTW, if you want to make your code more readable for an international audience you should use English variable names, rather than Portuguese. $\endgroup$ – PM 2Ring Jan 23 '19 at 19:19
  • 1
    $\begingroup$ 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, 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, which is much more precise. $\endgroup$ – PM 2Ring Jan 23 '19 at 19:26
  • $\begingroup$ I suggest you try improved Euler (en.wikipedia.org/wiki/Heun%27s_method). It improves the accuracy a lot. $\endgroup$ – irchans Jan 24 '19 at 3:22
4
$\begingroup$

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.

| improve this answer | |
$\endgroup$
1
$\begingroup$

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

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – Anna Paula Mendes Jan 23 '19 at 20:27
  • $\begingroup$ @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: i.stack.imgur.com/liJ4o.png $\endgroup$ – uhoh Jan 23 '19 at 22:24

Not the answer you're looking for? Browse other questions tagged or ask your own question.