Vi tar utgangspunkt i Simulering med pygame 5 som reknar ut banen til eit objekt under påvirkning av gravitasjonskrafta frå eit anna objekt. No skal vi rekna så nøye vi kan på banen til jorda under påvirkning av gravitasjonskrafta frå sola. Men vi har fremdeles begrensningar i forhold til 1) metode (vi bruker Eulermetoden for å finna fart og akselerasjon) og 2) Vi betraktar dette som eit "eit-legeme"-problem. Det vil sei at vi lar sola stå i ro. I verkeligheten er jo også sola påvirka av jorda, slik at dei to vil gå i ellipsebanar rundet sitt felles massesenter. Sjå tolegemeproblemet.
Men vi skal sjå at vi likevel får brukbare tal frå programmet vårt. Då må vi putta inn dei reelle tala for jord-sol-systemet, vi må bruka riktige enheter, og vi må ta med tida i likningane våre. Så vi begynner med startposisjonen og startfarten til jorda. Jorda går i ein ellipsebane, og vi bruker verdiane frå det punktet der jord er nærmast sola (dette kallast for perihel eller perihelium). Avstanden er i km og farten er i km/s:
x = -147.10e6 # avstand i km i perihel
y = 0 # Startverdi i y-retning.
vx = 0
vy = 30.29 # fart i km/s i perihel
Vi lar no sola vera i origo. Det betyr at koordinatane til vektoren frå jorda til sola kan skrivast som:
rx = -x
ry = -y
Den neste forandringen handlar om korleis vi reknar ut krafta på jorda.
Gravitasjonsloven kan skrivast:
$ F=G\frac{{m}{M}}{{r^2}} $, der m er massen til jorda og M er massen til sola. Vi kan rekna ut akselerasjonen til jorda som $ a=\frac{G}{{m}} $. Det gir oss
$ a=G\frac{{M}}{{r^2}} $
Så vi ser at vi ikkje egentlig treng jordmassen i det heile. Vi ser også at vi kan multiplisera saman M og G til ein konstant MG. Så likninga vår er no:
$ a=\frac{{MG}}{{r^2}} $
Når ein reknar avstandar i kilometer, er denne konstanten lik:
MG = 1.32712440018e11
Når vi har rekna ut $r$, $r^3$ og , så kan vi finna komponentane til akselerasjonen:
ax = MG/r3*rx
ay = MG/r3*ry
Så må vi ta hensyn til tida i likningane våre. Så komponentane til farten blir:
vx += ax*dt
vy += ay*dt
Her er dt tidssteget som vi set til dt = 3600. (enhet sekund) Og på same måten kan vi finna posisjonen ved:
x += vx*dt
y += vy*dt
No kan vi også summera tidsstega. Så vi definerer ein variabel T, som vi set lik 0 før simuleringen startar, og så summerer vi:
T += dt
Men no har vi jo koordinater som ikkje lenger passar i vinduet vårt. Så for å få det til må vi rekna ut nye x- og y-verdiar. Disse kallar vi xpos og ypos. Hvis vi tenkjer at ein pixel er det same som ein million km, så må vi dela koordinatane med ein million. Men så må vi også leggja til 500 i x-retning og 200 i y-retning for å få origo i sentrum av vinduet:
xpos = int(x/1e6) + 500
ypos = int(y/1e6) + 200
Så legg vi også inn litt kode som skal skriva ut og nullstilla T for kvart omløp. Dermed skal T liggja i nærheten av eit år i sekunder. I tillegg vil vi gjerne ha tal for kor store store og lille halvakse i ellipsen er, samt eksentrisiteten. Sjå Keplers lover og Ellipsen. Dermed ser heile koden vår slik ut:
import pygame
import math
pygame.init()
BREDDE = 1000 # Bredden på vinduet
HOYDE = 400 # Høyden på vinduet
SOLRADIUS = 10
JORDRADIUS = 4
KOLLDIST = SOLRADIUS + JORDRADIUS
MG = 1.32712440018e11 # Standard gravitational parameter [km3/s2]
dt = 3600 # Delta-t i sekund
JORD = (155, 255, 200)
SOL = (255,255,0)
HIMMEL = (0,0,100)
FPS = 200 #Frames Per Second, dvs kor fort animasjonen skal gå.
vindu = pygame.display.set_mode((BREDDE, HOYDE))
pygame.display.set_caption("Jordbanen: ")
clock = pygame.time.Clock()
def game_loop():
# Koordinatane til jorda. Vi simulerer ikkje sola som ståri ro i origo
x = -147.10e6 # avstand i km i perihel
y = 0 # Startverdi i y-retning.
vx = 0
vy = 30.29 # fart i km/s i perihel
running = True
T = 0
forrigey = 1
xmin = 200e6
xmax = -200e6
ymin = 200e6
ymax = -200e6
vindu.fill(HIMMEL)
while running:
for event in pygame.event.get():
if event.type == pygame.QUIT:
return pygame.quit()
# Finn kvadratet av avstanden (r2) og avstanden (r)
# r-vektor går frå jorda til origo
rx = -x # x1 i origo
ry = -y # y1 i origi
r2 = (rx)**2 + (ry)**2
r = math.sqrt(r2)
r3 = r*r*r
# Finn komponentane til jordens akselerasjon direkte:
ax = MG/r3*rx
ay = MG/r3*ry
# Finn farten i x- og y-retning
vx += ax*dt
vy += ay*dt
#Finn ny posisjon (x- og y-koordinater)
x += vx*dt
y += vy*dt
T += dt
xpos = int(x/1e6) + 500
ypos = int(y/1e6) + 200
# Tegnar systemet
vindu.fill(HIMMEL)
pygame.draw.circle(vindu, SOL, (500, 200), SOLRADIUS, 0)
pygame.draw.circle(vindu, JORD,(xpos,ypos), JORDRADIUS, 0)
pygame.display.update()
#clock.tick(FPS)
if x > xmax:
xmax = x
if x < xmin:
xmin = x
if y > ymax:
ymax = y
if y < ymin:
ymin = y
# Utskrift ein gong per omløp
if y < 0 and forrigey > 0 or y > 0 and forrigey < 0:
if x < 0:
print("Tid",T,"s = ",T/86400,"dagar")
T = 0
forrigey = y
a = (xmax-xmin)/2 # Store halvakse
b = (ymax-ymin)/2 # Lille halvakse
if b < a:
e = math.sqrt(1-(b/a)**2) # Eksentrisitet
else:
e = 0
print("Store halvakse:",a)
print("Lille halvakse:",b)
print("Eksentrisitet:",e) # Skal vera ca 0.0167
game_loop()
while True:
for event in pygame.event.get():
if event.type == pygame.QUIT:
pygame.quit()