#! /bin/sh
mp=$1		# mass of planet
mm=$2		# mass of moon
xm=$3		# initial x-position of moon
ydp=$4		# initial y-velocity of planet
ydm=$5		# initial y-velocity of moon (relative to planet, scaled to vc)
wait=$6

date=`date '+%D %H:%M'`

sm -s << FIN
if ($wait) {xterm_s} else {laser_sp 'moon:mp$mp:mm$mm:xm$xm:ydp$ydp:ydm$ydm.ps' \n lweight 2}

expand 1.15
location 4000 30000 4000 30000 
limits -2 2 -2 2
box

# format of data file is time, then (x,y) for the star, planet, moon
da moon:mp$mp:mm$mm:xm$xm:ydp$ydp:ydm$ydm
read {t 1 x0 2 y0 3 x1 4 y1 5 x2 6 y2 7}

# find the planet period by seeing where it first goes to y<0
define i 2 
while {y1[\$i]>0}  {
    define i (\$i+1)
}
define tper (2*t[\$i])
echo planet period=\$tper

# plot the star orbit
ctype 5
connect x0 y0

# plot the planet orbit
ctype 0
connect x1 y1 

# plot the moon orbit advancing the color type once each planet period
ctype 3
connect x2 y2 if (t <= \$tper)
ctype 4
connect x2 y2 if (t > \$tper && t <= 2*\$tper)
ctype 5
connect x2 y2 if (t > 2*\$tper && t <= 3*\$tper)
ctype 6
connect x2 y2 if (t > 3*\$tper && t <= 4*\$tper)
ctype 7
connect x2 y2 if (t > 4*\$tper)

# label the main box
ctype 0
limits 0 10 0 10
relocate 0 10.2
expand 1.15
putlabel 6 init: m_p=$mp, m_m=$mm, x_m=$xm, \dot{y}_p/v_c=$ydp, \dot{y}_m/v_c=$ydm
define tmax (t[dimen(t)-1])
relocate 0.5 9.5
putlabel 6 T_{planet}=\$tper 
relocate 0.5 9.0
putlabel 6 t_{max}=\$tmax
relocate 5 -1.2
putlabel 5 x
relocate -1.2 5
putlabel 5 y

# plot an inset showing the moon position relative to the planet
location 23000 29500 23000 29500
expand 0.8
define dw (1.5*($xm-1))
limits -\$dw \$dw -\$dw \$dw
box
ctype 3
connect (x2-x1) (y2-y1) if (t <= 1\$tper)
ctype 4
connect (x2-x1) (y2-y1) if (t > \$tper && t <= 2*\$tper)
ctype 5
connect (x2-x1) (y2-y1) if (t > 2*\$tper && t <= 3*\$tper)
ctype 6
connect (x2-x1) (y2-y1) if (t > 3*\$tper && t <= 4*\$tper)
ctype 7
connect (x2-x1) (y2-y1) if (t > 4*\$tper)
ctype 0

# label this inset
limits 0 10 0 10
relocate 0 10.5
relocate 5 -1.8
putlabel 5 \Delta x
relocate -1.8 5
putlabel 5 \Delta y

# plot an inset that is zoomed around the initial position
location 23000 29500 5000 11500
limits 0.85 1.15 -0.15 0.15
box
connect x1 y1

ctype 3
connect x2 y2 if (t <= \$tper)
ctype 4
connect x2 y2 if (t > \$tper && t <= 2*\$tper)
ctype 5
connect x2 y2 if (t > 2*\$tper && t <= 3*\$tper)
ctype 6
connect x2 y2 if (t > 3*\$tper && t <= 4*\$tper)
ctype 7
connect x2 y2 if (t > 4*\$tper)
ctype 0
ctype 0
#putdate "$date"
if ($wait) {!sleep $wait} else {hardcopy}
FIN