Issue
Im am trying to produce the shortest distance between the two orbiting points (Earth and Jupiter) created by this orbital model. I'vebeen working on this for quite some time but have been struggling.
Could someone sugest a way to create this output or possibly create a line connecting the two points and display their distance reative to one another? Ideally in km but AU is fine too.
I'm new to programming so any help would be great.
Thank you in advance!
This is the code ive been using on Google Colab;
`
# -*- coding: utf-8 -*-
"""Untitled2.ipynb
Automatically generated by Colaboratory.
Original file is located at
https://colab.research.google.com/drive/14vpwJ9ixq6YZGSxN2eBWay-hIPe__BxF
"""
#%% plot it
import matplotlib.pyplot as plt
from matplotlib import animation
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 2**128
#matplotlib.use("TkAgg") # for mac M1
from IPython.display import HTML
import matplotlib as mpl
mpl.rcParams.update(mpl.rcParamsDefault)
from tempfile import TemporaryFile
# First to define the Constants
G = 6.67e-11
Ms = 1.988e30 # Sun
Me = 5.972e24 # Earth
Mj = 1.898e27 # Jupiter
AU = 1.5e11
daysec = 24.0*60*60
e_ap_v = 29780 # Earth's velocity
j_ap_v = 13060 # Jupiter's velocity
gravconst_e = G*Me*Ms
gravconst_j = G*Mj*Ms
# To setup the starting conditions and locations on plot of each planet;
# Earth
xe,ye,ze = 0,(1*AU),0
xve,yve,zve = -e_ap_v,0,0
# Jupiter
xj,yj,zj = 0,(5.2*AU),0
xvj,yvj,zvj = -j_ap_v,0,0
# Sun
xs,ys,zs = 0,0,0
xvs,yvs,zvs = 0,0,0
t = 0.0
dt = 1*daysec # every frame move this time
xelist,yelist,zelist = [],[],[]
xslist,yslist,zslist = [],[],[]
xjlist,yjlist,zjlist = [],[],[]
# start simulation
while t<1*365*daysec:
################ earth #############
# compute G force on earth
#rx,ry,rz = xs - xe, ys - ye, zs - ze
rx,ry,rz = xe - xs, ye - ys, ze - zs
modr3_e = (rx**2+ry**2+rz**2)**1.5
fx_e = -gravconst_e*rx/modr3_e
fy_e = -gravconst_e*ry/modr3_e
fz_e = -gravconst_e*rz/modr3_e
# update quantities how is this calculated? F = ma -> a = F/m
xve += fx_e*dt/Me
yve += fy_e*dt/Me
zve += fz_e*dt/Me
# update position
xe += xve*dt
ye += yve*dt
ze += zve*dt
# save the position in list
xelist.append(xe)
yelist.append(ye)
zelist.append(ze)
################ Jupiter ##############
# compute G force on Jupiter
rx_j,ry_j,rz_j = xj - xs, yj - ys, zj - zs
modr3_j = (rx_j**2+ry_j**2+rz_j**2)**1.5
fx_j = -gravconst_j*rx_j/modr3_j
fy_j = -gravconst_j*ry_j/modr3_j
fz_j = -gravconst_j*rz_j/modr3_j
xvj += fx_j*dt/Mj
yvj += fy_j*dt/Mj
zvj += fz_j*dt/Mj
# update position
xj += xvj*dt
yj += yvj*dt
zj += zvj*dt
# add to list
xjlist.append(xj)
yjlist.append(yj)
zjlist.append(zj)
################ the sun ###########
# update quantities how is this calculated? F = ma -> a = F/m
xvs += -(fx_e+fx_j)*dt/Ms
yvs += -(fy_e+fy_j)*dt/Ms
zvs += -(fz_e+fz_j)*dt/Ms
# update position
xs += xvs*dt
ys += yvs*dt
zs += zvs*dt
xslist.append(xs)
yslist.append(ys)
zslist.append(zs)
# update dt
t +=dt
# to plot the data
import matplotlib.pyplot as plt
from matplotlib import animation
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 2**128
from IPython.display import HTML
fig, ax = plt.subplots(figsize=(5.8,5.8))
ax.set_aspect('equal')
# Defining the planets and orbit path attributes on the graph (size & colour)
line_e, = ax.plot([],[],'b',lw=3,)
point_e, = ax.plot([AU], [0], marker="o", markersize=5, markeredgecolor="green", markerfacecolor="blue")
text_e = ax.text(AU,0,'Earth')
line_j, = ax.plot([],[],'r',lw=3)
point_j, = ax.plot([5.2*AU], [0], marker="o", markersize=7, markeredgecolor="red", markerfacecolor="brown")
text_j = ax.text(5.2*AU,0,'Jupiter')
point_s, = ax.plot([0], [0], marker="o", markersize=10, markeredgecolor="orange", markerfacecolor="yellow")
text_s = ax.text(0,0,'Sun')
exdata,eydata = [],[] # earth track
sxdata,sydata = [],[] # sun track
jxdata,jydata = [],[] # Jupiters track
print(len(xelist))
def update(i):
exdata.append(xelist[i])
eydata.append(yelist[i])
jxdata.append(xjlist[i])
jydata.append(yjlist[i])
line_e.set_data(exdata,eydata)
point_e.set_data(xelist[i],yelist[i])
text_e.set_position((xelist[i],yelist[i]))
line_j.set_data(jxdata,jydata)
point_j.set_data(xjlist[i],yjlist[i])
text_j.set_position((xjlist[i],yjlist[i]))
point_s.set_data(xslist[i],yslist[i])
text_s.set_position((xslist[i],yslist[i]))
ax.set_xlim(-5.8*AU,5.8*AU)
ax.set_ylim(-5.8*AU,5.8*AU)
return line_e,point_s,point_e,line_j,point_j,text_e,text_j,text_s,
anim = animation.FuncAnimation(fig,func=update,frames=len(xelist),interval=1,blit=True)
# Showing animation in Jupyter Notebook
from IPython.display import HTML
HTML(anim.to_jshtml())
Solution
The connection line between Jupiter and the Earth can be drawn using the arrow method of matplotlib. This allows to add arrow heads if desired. The distance can be added as a normal label.
annotation = ax.arrow(0, AU, 0, (5.2 - 1) * AU)
text_a = ax.text(0, (5.2 + 1) / 2 * AU, '', ha='center', backgroundcolor='white')
In the update method the distance has to be calculated using the Pythagorean theorem:
annotation.set_data(x=exdata[-1], y=eydata[-1], dx=jxdata[-1] - exdata[-1], dy=jydata[-1] - eydata[-1])
text_a.set_position(((exdata[-1] + jxdata[-1]) / 2, (eydata[-1] + jydata[-1]) / 2))
text_a.set_text(f"{dist/1e9:.1f} Mio km")
Whole code:
# -*- coding: utf-8 -*-
"""Untitled2.ipynb
Automatically generated by Colaboratory.
Original file is located at
https://colab.research.google.com/drive/14vpwJ9ixq6YZGSxN2eBWay-hIPe__BxF
"""
#%% plot it
import matplotlib.pyplot as plt
from matplotlib import animation
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 2**128
#matplotlib.use("TkAgg") # for mac M1
from IPython.display import HTML
from math import sqrt
import matplotlib as mpl
mpl.rcParams.update(mpl.rcParamsDefault)
from tempfile import TemporaryFile
# First to define the Constants
G = 6.67e-11
Ms = 1.988e30 # Sun
Me = 5.972e24 # Earth
Mj = 1.898e27 # Jupiter
AU = 1.5e11
daysec = 24.0*60*60
e_ap_v = 29780 # Earth's velocity
j_ap_v = 13060 # Jupiter's velocity
gravconst_e = G*Me*Ms
gravconst_j = G*Mj*Ms
# To setup the starting conditions and locations on plot of each planet;
# Earth
xe,ye,ze = 0,(1*AU),0
xve,yve,zve = -e_ap_v,0,0
# Jupiter
xj,yj,zj = 0,(5.2*AU),0
xvj,yvj,zvj = -j_ap_v,0,0
# Sun
xs,ys,zs = 0,0,0
xvs,yvs,zvs = 0,0,0
t = 0.0
dt = 1*daysec # every frame move this time
xelist,yelist,zelist = [],[],[]
xslist,yslist,zslist = [],[],[]
xjlist,yjlist,zjlist = [],[],[]
# start simulation
while t<1*365*daysec:
################ earth #############
# compute G force on earth
#rx,ry,rz = xs - xe, ys - ye, zs - ze
rx,ry,rz = xe - xs, ye - ys, ze - zs
modr3_e = (rx**2+ry**2+rz**2)**1.5
fx_e = -gravconst_e*rx/modr3_e
fy_e = -gravconst_e*ry/modr3_e
fz_e = -gravconst_e*rz/modr3_e
# update quantities how is this calculated? F = ma -> a = F/m
xve += fx_e*dt/Me
yve += fy_e*dt/Me
zve += fz_e*dt/Me
# update position
xe += xve*dt
ye += yve*dt
ze += zve*dt
# save the position in list
xelist.append(xe)
yelist.append(ye)
zelist.append(ze)
################ Jupiter ##############
# compute G force on Jupiter
rx_j,ry_j,rz_j = xj - xs, yj - ys, zj - zs
modr3_j = (rx_j**2+ry_j**2+rz_j**2)**1.5
fx_j = -gravconst_j*rx_j/modr3_j
fy_j = -gravconst_j*ry_j/modr3_j
fz_j = -gravconst_j*rz_j/modr3_j
xvj += fx_j*dt/Mj
yvj += fy_j*dt/Mj
zvj += fz_j*dt/Mj
# update position
xj += xvj*dt
yj += yvj*dt
zj += zvj*dt
# add to list
xjlist.append(xj)
yjlist.append(yj)
zjlist.append(zj)
################ the sun ###########
# update quantities how is this calculated? F = ma -> a = F/m
xvs += -(fx_e+fx_j)*dt/Ms
yvs += -(fy_e+fy_j)*dt/Ms
zvs += -(fz_e+fz_j)*dt/Ms
# update position
xs += xvs*dt
ys += yvs*dt
zs += zvs*dt
xslist.append(xs)
yslist.append(ys)
zslist.append(zs)
# update dt
t +=dt
# to plot the data
import matplotlib.pyplot as plt
from matplotlib import animation
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 2**128
from IPython.display import HTML
fig, ax = plt.subplots(figsize=(5.8,5.8))
ax.set_aspect('equal')
# Defining the planets and orbit path attributes on the graph (size & colour)
line_e, = ax.plot([],[],'b',lw=3,)
point_e, = ax.plot([AU], [0], marker="o", markersize=5, markeredgecolor="green", markerfacecolor="blue")
text_e = ax.text(AU,0,'Earth')
line_j, = ax.plot([],[],'r',lw=3)
point_j, = ax.plot([5.2*AU], [0], marker="o", markersize=7, markeredgecolor="red", markerfacecolor="brown")
text_j = ax.text(5.2*AU,0,'Jupiter')
point_s, = ax.plot([0], [0], marker="o", markersize=10, markeredgecolor="orange", markerfacecolor="yellow")
text_s = ax.text(0,0,'Sun')
annotation = ax.arrow(0, AU, 0, (5.2 - 1) * AU)
text_a = ax.text(0, (5.2 + 1) / 2 * AU, '', ha='center', backgroundcolor='white')
exdata,eydata = [],[] # earth track
sxdata,sydata = [],[] # sun track
jxdata,jydata = [],[] # Jupiters track
print(len(xelist))
def update(i):
exdata.append(xelist[i])
eydata.append(yelist[i])
jxdata.append(xjlist[i])
jydata.append(yjlist[i])
line_e.set_data(exdata,eydata)
point_e.set_data(xelist[i],yelist[i])
text_e.set_position((xelist[i],yelist[i]))
line_j.set_data(jxdata,jydata)
point_j.set_data(xjlist[i],yjlist[i])
text_j.set_position((xjlist[i],yjlist[i]))
point_s.set_data(xslist[i],yslist[i])
text_s.set_position((xslist[i],yslist[i]))
ax.set_xlim(-5.8*AU,5.8*AU)
ax.set_ylim(-5.8*AU,5.8*AU)
dist = sqrt((exdata[-1] - jxdata[-1])**2 + (eydata[-1] - jydata[-1])**2)
annotation.set_data(x=exdata[-1], y=eydata[-1], dx=jxdata[-1] - exdata[-1], dy=jydata[-1] - eydata[-1])
text_a.set_position(((exdata[-1] + jxdata[-1]) / 2, (eydata[-1] + jydata[-1]) / 2))
text_a.set_text(f"{dist/1e9:.1f} Mio km")
return line_e,point_s,point_e,line_j,point_j,text_e,text_j,text_s,
anim = animation.FuncAnimation(fig,func=update,frames=len(xelist),interval=1,blit=True)
# Showing animation in Jupyter Notebook
from IPython.display import HTML
HTML(anim.to_jshtml())
Answered By - M472
0 comments:
Post a Comment
Note: Only a member of this blog may post a comment.