Skip to content

Instantly share code, notes, and snippets.

# FilipDominec/Lagrangian points - force plot Created Jan 25, 2015

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")
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.