Last active
October 19, 2024 07:31
-
-
Save ckhung/5015a61855593dc8d31d21ebbc0afe53 to your computer and use it in GitHub Desktop.
lagrange points as seen using effective potential plot
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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