# A gnuplot programme to display equipotential surfaces 
# for a Sun--Jupiter-like system (where both large bodies stay on their 
# places all the time) wrote on the Eastern Friday 2001 
# Jenik Hollan, N. Copernicus Observatory and Planetarium in Brno.

set output 'libr_pnt.ps'
set term postscript
 color

set nokey

set contour surface
set cntrparam levels 60
set cntrparam levels incremental -0.85,0.005

#set hidden3d

 G=1.0
 m0=1.0
 m2=1.0
 x2=5.0
 x1= x2 * m2 / (m0+m2)

 omega2= G * (m0 + m2) / (x2**3)

 Pot(x,y)=  \
    - omega2*((x-x1)**2+y**2)/2 \
    - G*m2/sqrt((x-x2)**2+y**2) \
    - G*m0/sqrt(x**2+y**2)


 xh=x2/2.0
 L4=x2*sqrt(3)/2.0

 d=8.0 # radius to plot

 xhl=xh-d
 xhu=xh+d
 Lyl=L4-d
 Lyu=L4+d/1.5

set view 45, 45
#set hidden3d

set title 'L1,2,4, equal masses'
splot [xhl:xhu] [Lyl:Lyu] [-2:0] Pot(x,y) 
set title 'L1,2 equal masses'
splot [x1-2:10] [-3:3] [-2:0] Pot(x,y) 

 d=6.0 # radius to plot

 xhl=xh-d
 xhu=xh+d
 Lyl=L4-d
 Lyu=L4+d/1.5

set cntrparam levels incremental -0.4,0.002

m2=0.1
x1= x2 * m2 / (m0+m2)
omega2= G * (m0 + m2) / (x2**3)
set title 'L1,2,4, m2=0.1'
splot [xhl:xhu] [Lyl:Lyu] [-1:-0.3] Pot(x,y) 
set title 'L1,2 m2=0.1'
splot [x1:8] [-3:3] [-1:-0.3] Pot(x,y)

set cntrparam levels incremental -0.35,0.001

m2=0.01
x1= x2 * m2 / (m0+m2)
omega2= G * (m0 + m2) / (x2**3)
set title 'L1,2,4, m2=0.01'
splot [xhl:xhu] [Lyl:Lyu] [-1:-0.3] Pot(x,y) 
set title 'L1,2 m2=0.01'
splot [x1+3:6.5] [-2:2] [-0.6:-0.3] Pot(x,y)


m2=0.001
x1= x2 * m2 / (m0+m2)
omega2= G * (m0 + m2) / (x2**3)
set  title 'L1,2,4, m2=0.001'
splot [xhl:xhu] [Lyl:Lyu] [-1:-0.3] Pot(x,y)
set cntrparam levels incremental -0.320,0.0004
set   title 'L1(2) m2=0.001'
splot [x1+4:5.5] [-1.5:1.5] [-0.5:-0.3] Pot(x,y)