Skip to content

Instantly share code, notes, and snippets.

@ckhung
Last active October 19, 2024 07:31
Show Gist options
  • Save ckhung/5015a61855593dc8d31d21ebbc0afe53 to your computer and use it in GitHub Desktop.
Save ckhung/5015a61855593dc8d31d21ebbc0afe53 to your computer and use it in GitHub Desktop.
lagrange points as seen using effective potential plot
# lagrange points as seen using effective potential plot
len(x,y) = sqrt(x*x+y*y)
# mass of planet relative to that of star
alpha = 0.01 # "mu" in wikipedia
# mass of star and mass of planet
Ms = 1000
Mp = Ms * alpha
# Let the barycenter be the origin and the axis of rotation be the z-axis
xb = 0
# Let the star and planet be on the X-axis, separated by distance r12
r12 = 1000
xs = - r12 * Mp / (Ms + Mp)
xp = r12 * Ms / (Ms + Mp)
# potential contributed by gravity
Pg(x,y,m,x0,y0) = - m / len(x-x0,y-y0)
# truncated
gfloor = -3
Pg_t(x,y,m,x0,y0) = Pg(x,y,m,x0,y0) < gfloor ? gfloor : Pg(x,y,m,x0,y0)
# now add centrifugal (pseudo) potential
pot_eff(x,y) = Pg_t(x,y,Ms,xs,0) + Pg_t(x,y,Mp,xp,0) - (Ms+Mp)/2*((x-xb)*(x-xb)+y*y)/r12**3
set pm3d
set contour
set cntrparam levels 30
set isosamples 50,50
set view 0,0
rot0 = pi/6
splot 0
pause -1 "maximize graphics window, then press enter to start\n"
set xrange [-0.5*r12 : 2*r12]
set yrange [-0.5*r12 : 2*r12]
splot pot_eff(x*cos(rot0)+y*sin(rot0), -x*sin(rot0)+cos(rot0)*y)
pause -1 "overview\n"
set xrange [xp*cos(rot0)*0.7 : xp*cos(rot0)*1.3]
set yrange [xp*sin(rot0)*0.7 : xp*sin(rot0)*1.3]
set view 0,0
replot
pause -1 "planet view\n"
# half of view width
hvw = r12 * 0.05
# near L1
delta = (alpha/3) ** (1.0/3)
x_L1 = r12 * (1-delta) + xs
set xrange [x_L1*cos(rot0) - hvw : x_L1*cos(rot0) + hvw]
set yrange [x_L1*sin(rot0) - 0.6*hvw : x_L1*sin(rot0) + 0.6*hvw]
set view 0,0
replot
pause -1 "near L1\n"
delta = (alpha/3) ** (1.0/3)
x_L2 = r12 * (1+delta) + xs
set xrange [x_L2*cos(rot0) - hvw : x_L2*cos(rot0) + hvw]
set yrange [x_L2*sin(rot0) - 0.6*hvw : x_L2*sin(rot0) + 0.6*hvw]
set view 0,0
replot
pause -1 "near L2\n"
x_L3 = - r12 * (1 + 5/12 * alpha) + xs
set xrange [x_L3*cos(rot0) - hvw : x_L3*cos(rot0) + hvw]
set yrange [x_L3*sin(rot0) - 0.6*hvw : x_L3*sin(rot0) + 0.6*hvw]
set view 0,0
replot
pause -1 "near L3\n"
x_L4 = xs * sqrt(3)/2
y_L4 = r12 + xs/2
set xrange [x_L4-hvw : x_L4+hvw]
set yrange [y_L4-0.3*hvw : y_L4+0.3*hvw]
set view 0,0
replot
pause -1 "near L4\n"
set xrange [x_L4-0.1*hvw : x_L4+0.1*hvw]
set yrange [y_L4-0.03*hvw : y_L4+0.03*hvw]
set view 0,0
replot
pause -1 "close up view of L4\n"
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment