#
# Plot the velocity & stuff during the trip from Earth to AC
#
#    
# set terminal x11 0
# set output

# We need Earth time, at the traveler's location, as a function of tau.
# Amazingly, in that thick haze of algebra, we never computed that,
# and it looks a bit messy.  So we'll do it the parametric way instead.
#

set parametric

set key left
set trange [0:11.32]
##set xrange [0:7.04]
set xrange [0:11.32]
set yrange [*:*]
set grid

set title "View Seen via Telescope, From Earth"
set xlabel "Earth Time (on Clock Next to Telescope)"

# Earth time at turnover and arrival
tt = 2.83
tf = 5.66

# Ship time at turnover and arrival
taut = 1.76
tauf = 3.52

# Distance at turnover and arrival
xt = 2
xf = 4

# Complemented time
tc(t) = tf - t

# Velocity before and after turnover, and combined
v0(t) = t/sqrt(1 + t**2)
v1(t) = tc(t)/sqrt(1 + tc(t)**2)
v(t) = t<(tf/2) ? v0(t) : v1(t)

# Ship time before and after turnover, and combined
tau0(t) = log(t + sqrt(t**2 + 1))
tau1(t) = tauf - tau0(tf - t)
tau(t) = (t <= tt)	 ? tau0(t) : \
	 (t <= tf)	 ? tau1(t) : \
	 (t <= tf + tt)	 ? tauf + tau0(t - tf) : \
	 tauf + tau1(t - tf)

# Ship location before and after turnover, and combined
dx0(t) = sqrt(t**2 + 1) - 1
dx1(t) = xf - dx0(tc(t))
distance(t) = (t <= tt) ? dx0(t) : dx1(t)

# Gamma
Gamma0(t) = 1/sqrt(1 - v(t)**2)
Gamma(t) = (t <= tf) ? Gamma0(t) : Gamma0(t - tf)

# Photon arrival time as a function of ship's Earth time coordinate
# Outbound
t_a0(t) = exp(tau(t)) - 1
t_a1(t) = tf + xf - (exp(tauf - tau(t)) - 1)
#
# Returning
t_a2(t) = xf - (sqrt(1 + t**2) - t - 1)
t_a3(t) = t + sqrt(tc(t)**2 + 1) - 1
#
# Combined
t_a(t) = (t <= tt)	? t_a0(t)		: \
       (t <= tf)	? t_a1(t)		: \
       (t <= tf+tt)	? tf + t_a2(t - tf)	: \
       tf + t_a3(t - tf)

# Photon emission time in ship's frame, as a function of Earth time
# Outbound it's just tau(t).
# Returning, it's just tau(t) + tauf
tau_e(t) = (t <= tf)   ? tau(t)	: tauf + tau(t - tf)

# Rate of clock seen in telescope versus clock rate on Earth
# Outbound
dtau_dt0(t) = 1/(Gamma0(t) * (1 + v(t)))
# Returning
dtau_dt1(t) = 1/(Gamma0(t) * (1 - v(t)))
# Combined
dtau_dt(t) = (t <= tf)	? dtau_dt0(t)	: dtau_dt1(t - tf)

set xtics 1
set ytics 1

plot \
[t = ] \
   t_a(t), tau(t) with lines lw 2, \
   t_a(t), dtau_dt(t) with lines lw 2, \
   t_a(t), t, \
   t_a(t), t_a(t), \
   t_a(t), Gamma(t)