Instantly share code, notes, and snippets.

Embed
What would you like to do?
there are 5 points in the Sun-Earth system where gravitational and centrifugal forces cancel
#!/usr/bin/env python
#coding:utf8
import numpy as np
#Mearth, Msun = 5.97e24, 1.98e30 ## <-- these are realistic
Mearth, Msun = 1.97e29, 1.98e30 ## <-- unrealistic, but better for visualisation
earth_sun_dist = 149e9
Xsun = - earth_sun_dist * (Mearth/Msun) ## Sun center distance to the center-of-mass
Xearth = earth_sun_dist * (1 - Mearth/Msun) ## Earth center distance to the center-of-mass
kappa = 6.67e-11 ## Gravitational constant
year = 365.25*24*60*60 ## Rotation period [seconds]
## Generate the points in the space (takes a second)
xaxis, yaxis = np.linspace(-300e9, 300e9, 1000), np.linspace(-300e9, 300e9, 1000)
x,y = np.meshgrid(xaxis,yaxis)
## Centrifugal force
rC = (x**2+y**2)**.5
FC = (2*np.pi / year)**2 * rC
FCx = FC * (x/rC)
FCy = FC * (y/rC)
## Gravity of Sun
rS = ((x-Xsun)**2+y**2)**.5
FS = -kappa * Msun / (rS**2)
FSx = FS * (x-Xsun)/rS
FSy = FS * y/rS
## Gravity of Earth
rE = ((x-Xearth)**2+y**2)**.5
FE = -kappa * Mearth / (rE**2)
FEx = FE * (x-Xearth)/rE
FEy = FE * y/rE
## All three forces together
Fax = FCx+FSx+FEx
Fay = FCy+FSy+FEy
Fa = (Fax**2+Fay**2)**.5
import matplotlib.pyplot as plt
d = 15 # decimation factor
## Plot the magnitude of the force:
plt.contourf(xaxis, yaxis, np.log10(Fa), levels=np.linspace(-4, 0), extend='both')
## ... and its direction:
plt.quiver(x[::d, ::d],y[::d, ::d],Fax[::d, ::d]/Fa[::d, ::d], Fay[::d, ::d]/Fa[::d, ::d])
plt.show() # plt.savefig("output.png")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment