James Webb Space Telescope (JWST) ble skutt opp desember 2021, og har vekket mye interesse og oppmerksomhet for astronomi og hvilke fenomener den skal utforske. Men også banen til JWST har vekket manges nysgjerrighet, siden skal ligge i likevektspunktet som heter Lagrange 2 (eller L2). Her er et forsøk på å simulere banen til JWST når det har kommet seg til L2.
Innholdsfortegnelse
Teorien
Likningene og teorien vi bruker her er Newtons uttrykk for tyngdekraft mellom to himmellegemer, og bevegelseslikningene i sin diskrete form. Solen antas å være stasjonær, siden den har mye mer masse enn både Jorden og JWST. Banen til Jorden beregnes kun utifra Solens tyngdekraft, mens banen til JWST beregnes utifra både Solens og Jordens tyngdekraft.
Eulers metode
For å integrere opp bevegelseslikningene, har vi valgt enkleste tilnærming: Eulers metode. Den er rask å kode og lett å lese seg opp på for elevene. Nøyaktighet kan være en utfordring, men etter å ha prøvd med to ulike steglengder for tid, 24 timer og 1 minutt, og oppnådd omtrent samme resultat, er vi trygg på at metoden gir en akseptabel nøyaktighet for dette eksempelet.
Kodebiblioteker
I denne koden bruker bibliotekene matplotlib og numpy.
import numpy as np
import matplotlib.pyplot as plt
Fysiske konstanter og variable
Fysiske konstanter for å regne ut Newtons tyngdekraftslov, objektenes masse og avstand og hastighet som beskriver banen til Jorden og forventet hastighet og posisjon JWST må ha for å kunne ligge i L2, defineres.
# Definer generelle variable
dt = 3600.0 # 1 hour
AE = 149.6 * 10**9 #in m
G = 6.67 * 10**(-11) #unit
# Sett verdier for avstand (baneradius) og hastighet
# for James Webb Space Telescope og Jorden (enhet m og m/s)
Rjord = 1.0 * AE
Vjord = 29.78 * 10**3
Rjames = 151.1 * 10**9
Vjames = Vjord * ( Rjames / Rjord )
# Sett verdier for objektenes masse (enhet kg)
Msol = 1.989 * 10**30
Mjord = 5.972 * 10**24
Mjames = 6500.0
Integrasjon (kode)
Integrasjonen av bevegelseslikningene og kraftberegningene gjøres i den sentrale while-løkken. Vi finner avstand mellom objektene (R) og hvilken retning de ligger i forhold til hverandre (gjennom deres enhetsvektorer, E), før vi beregner kraften (F) mellom dem.
Fra kreftene (F), finner vi akselerasjonen (A) som bi bruker til å oppdatere hastighet (V) og posisjon (P) med Eulers metode.
# La programmet kjøre til det blir avsluttet med Ctrl+C
while True:
# Oppdater tidssteg
tidssteg += 1
tid += dt
#Regn ut avstand (R) og enhetsvektor (E) mellom...
#James Webb og Solen:
Rjames = np.linalg.norm(Pjames)
Ejames = Pjames / Rjames
#Jorden og Solen:
Rjord = np.linalg.norm(Pjord)
Ejord = Pjord / Rjord
#James Webb og Jorden:
Rjames_jord = np.linalg.norm(Pjames-Pjord)
Ejames_jord = (Pjames-Pjord)/Rjames_jord
#Regn ut kraften som virker på Jorden fra Solen
Fjord = -Ejord * G*Msol*Mjord/Rjord**2
#Regn ut kraften som virker på James Webb fra Jorden og Solen
Fjames = -Ejames*G*Msol*Mjames/Rjames**2 \
- Ejames_jord*G*Mjord*Mjames/Rjames_jord**2
#Regn ut aksellerasjonen
Ajord = Fjord/Mjord
Ajames = Fjames/Mjames
#Regn ut hastigheten fra akselerasjonen og forrige tidssteg
Vjord += Ajord*dt
Vjames += Ajames*dt
#Regn ut posisjonen fra hastigheten og forrige tidsteg
Pjord += Vjord*dt
Pjames += Vjames*dt
Plotting
I denne koden har vi valgt å animere beregningene, som dessverre betyr at koden blir ganske mye mer omfattende.
Posisjonene lagres først i egne arrays, som «sendes» til plotte vinduet. Så justeres x- og y-grensene i plotvinduet i vinduet som viser et forstørret utsnitt av Jordens og JWST’s baner.
Resten av plot-kommandoene kan du se i hele koden (bla lenger ned).
# Plott data, men ikke hvert eneste tidssteg. Det tar så langt tid..
if tidssteg % 10 == 0: #Gjør dette hvert 10'ende tidsteg
#Lagre de nye posisjonene i en array (normer størrelsene til AE)
xjord.append(Pjord[0]/AE)
yjord.append(Pjord[1]/AE)
xjames.append(Pjames[0]/AE)
yjames.append(Pjames[1]/AE)
#Kort ned lengden av det som plottes til ett år
xjord = xjord[-plot_data_lengde:]
yjord = yjord[-plot_data_lengde:]
xjames = xjames[-plot_data_lengde:]
yjames = yjames[-plot_data_lengde:]
#Oppdater linjene som plottes
line2.set_xdata(xjord)
line2.set_ydata(yjord)
line4.set_xdata(xjord)
line4.set_ydata(yjord)
line3.set_xdata(xjames)
line3.set_ydata(yjames)
line5.set_xdata(xjames)
line5.set_ydata(yjames)
#Juster vinduet som følger James Webb og Jorden
ax2.set_xlim(Pjord[0]/AE-0.03,Pjord[0]/AE+0.03)
ax2.set_ylim(Pjord[1]/AE-0.03,Pjord[1]/AE+0.03)
#Tegn det nye plottet
fig.canvas.draw()
plt.pause(0.0001)
Last ned kode
Hele koden kan også leses/lastes ned her:
# PYTHONSKOLE.NO
# Simulering av Jorden og James Webb i Lagrange 2.
# 3.januar 2021, Vegard L. Rekaa
# vegard@astronomen.no
import numpy as np
import matplotlib.pyplot as plt
# Definer generelle variable
dt = 3600.0 # 1 hour
AE = 149.6 * 10**9 #in m
G = 6.67 * 10**(-11) #unit
# Sett verdier for avstand (baneradius) og hastighet
# for James Webb Space Telescope og Jorden (enhet m og m/s)
Rjord = 1.0 * AE
Vjord = 29.78 * 10**3
Rjames = 151.1 * 10**9
Vjames = Vjord * ( Rjames / Rjord )
# Sett verdier for objektenes masse (enhet kg)
Msol = 1.989 * 10**30
Mjord = 5.972 * 10**24
Mjames = 6500.0
# Deklarer posisjoner og hastighet som vektorer
# P = posisjon, V = hastighet
# NB Disse overskriver de gamle variablene Vjord og Vjames
Pjord = np.array([Rjord,0.0])
Vjord = np.array([0.0,Vjord])
Pjames = np.array([Rjames,0.0])
Vjames = np.array([0.0,Vjames])
# Deklarer tomme lister hvor vi kan lagre verdier for plotting
tid = []
xjord = []
yjord = []
xjames = []
yjames = []
#Lagre initielle / startverdier for simuleringen
tid = 0.0
xjord.append(Pjord[0]/AE)
yjord.append(Pjord[1]/AE)
xjames.append(Pjames[0]/AE)
yjames.append(Pjames[1]/AE)
# Noen tekniske linjer for å kunne animere plottet
plot_interval = 10
plt.ion()
fig, (ax1,ax2) = plt.subplots(1,2)
ax1.set_title("Normal skala")
ax1.set_xlabel("x (AE)")
ax1.set_ylabel("y (AE)")
ax1.set_aspect('equal')
ax1.set_xlim(-1.2,1.2)
ax1.set_ylim(-1.2,1.2)
ax2.set_title("Forstørret")
ax2.set_aspect('equal')
ax2.set_xlabel("x (AE)")
ax2.set_xlim(Pjord[0]/AE-0.1,Pjord[0]/AE+0.1)
ax2.set_ylim(Pjord[1]/AE-0.1,Pjord[1]/AE+0.1)
line1, = ax1.plot(0,0,'yo',label="Solen")
line2, = ax1.plot(xjord,yjord,'b-',label="Jorden")
line3, = ax1.plot(xjames,yjames,'r-',label="JWST")
line4, = ax2.plot(xjord,yjord,'b-',label="Jorden")
line5, = ax2.plot(xjames,yjames,'r-',label="JWST")
ax1.legend()
ax2.legend()
# Disse bruker vi for å normere tid (utskrift til terminal)
sekunder_per_aar = 60.0*60*24*365.24
sekunder_per_dag = 60.0*60*24
#Disse bruker vi til å bestemme hvor ofte vi skal plotte
tidssteg_per_aar = int(sekunder_per_aar / dt)
plot_data_lengde = int(tidssteg_per_aar/plot_interval)-10
# Denne variabler bruker vi til å holde orden på antall tidssteg
tidssteg = 0
# Alt er klart.
print("SIMULERING AV JAMES WEBB OG JORDEN")
print("pythonskole.no - jameswebb-l2.py")
print("----------------------------------")
print("Steglengde (tid) = ",dt,"s = ",dt/sekunder_per_dag,"dager")
print("Tidssteg per aar = ",tidssteg_per_aar)
print("Jorden")
print(" Posisjon: ",Pjord,"m")
print(" Hastighet: ",Vjord,"m/s")
print("James Webb Space Telescope")
print(" Posisjon: ",Pjames,"m")
print(" Hastighet: ",Vjames,"m/s")
print("")
print("")
# La programmet kjøre til det blir avsluttet med Ctrl+C
while True:
# Oppdater tidssteg
tidssteg += 1
tid += dt
#Regn ut avstand (R) og enhetsvektor (E) mellom...
#James Webb og Solen:
Rjames = np.linalg.norm(Pjames)
Ejames = Pjames / Rjames
#Jorden og Solen:
Rjord = np.linalg.norm(Pjord)
Ejord = Pjord / Rjord
#James Webb og Jorden:
Rjames_jord = np.linalg.norm(Pjames-Pjord)
Ejames_jord = (Pjames-Pjord)/Rjames_jord
#Regn ut kraften som virker på Jorden fra Solen
Fjord = -Ejord * G*Msol*Mjord/Rjord**2
#Regn ut kraften som virker på James Webb fra Jorden og Solen
Fjames = -Ejames*G*Msol*Mjames/Rjames**2 \
- Ejames_jord*G*Mjord*Mjames/Rjames_jord**2
#Regn ut aksellerasjonen
Ajord = Fjord/Mjord
Ajames = Fjames/Mjames
#Regn ut hastigheten fra akselerasjonen og forrige tidssteg
Vjord += Ajord*dt
Vjames += Ajames*dt
#Regn ut posisjonen fra hastigheten og forrige tidsteg
Pjord += Vjord*dt
Pjames += Vjames*dt
# Skriv ut antall dager som er simulert til terminal/konsoll
print("Dager="+str(round(tid/sekunder_per_dag))+
" Tidssteg="+str(tidssteg)+" ",end="\r")
# Plott data, men ikke hvert eneste tidssteg. Det tar så langt tid..
if tidssteg % 10 == 0: #Gjør dette hvert 10'ende tidsteg
#Lagre de nye posisjonene i en array (normer størrelsene til AE)
xjord.append(Pjord[0]/AE)
yjord.append(Pjord[1]/AE)
xjames.append(Pjames[0]/AE)
yjames.append(Pjames[1]/AE)
#Kort ned lengden av det som plottes til ett år
xjord = xjord[-plot_data_lengde:]
yjord = yjord[-plot_data_lengde:]
xjames = xjames[-plot_data_lengde:]
yjames = yjames[-plot_data_lengde:]
#Oppdater linjene som plottes
line2.set_xdata(xjord)
line2.set_ydata(yjord)
line4.set_xdata(xjord)
line4.set_ydata(yjord)
line3.set_xdata(xjames)
line3.set_ydata(yjames)
line5.set_xdata(xjames)
line5.set_ydata(yjames)
#Juster vinduet som følger James Webb og Jorden
ax2.set_xlim(Pjord[0]/AE-0.03,Pjord[0]/AE+0.03)
ax2.set_ylim(Pjord[1]/AE-0.03,Pjord[1]/AE+0.03)
#Tegn det nye plottet
fig.canvas.draw()
plt.pause(0.0001)