Chasing the Terminator on SQ23

I am on SQ23, a direct flight from New York City to Singapore, and I am procrastinating on my writing work by watching the plane chasing the day-night terminator on the in-flight progress screen. (I also watched the sun pseudo-rising for a quite while before being told to shut the window by the cabin crew).

Problem Statement

import matplotlib.pyplot as plt
import numpy as np

SQ23 departs from NYC at 2140h (local time) and lands in Singapore at 0416h (local time)¹ — at least, according to the projected time on the in-flight display. The date at the time of departure is 22nd November.

# for convenience we consider times relative to journey start in units of hours.
# our reference meridian for the problem will also be the starting location.

T = 18 # hours

Unfortunately, the real SQ23 doesn't actually follow the geodesic between JFK and SIN owing to political considerations (and I'm guessing because there are good tailwinds to be had by sticking to specific latitudes). However, supposing it did, and given the information above:

1. Does SQ23 ever catch up with the terminator? (the Sun still hasn't fully risen for me, so I don't actually know this yet; however, odds are looking good right now)²
2. If so, how many hours of astronomical daytime (i.e. altitude of Sun $>0$) am I going to get?

BONUS CHALLENGE: I have to do this problem without internet access because I didn't pay for in-flight WiFi…

Attempt at solution

We need to find the latitudes and longitudes of both the SQ23 and the terminator as a function of time. Unfortunately, since we're not at the equinox, the latitude of the terminator is also a function of longitude, which will have to be supplied by that of SQ23. Therefore, it makes more sense to parameterise the trajectory of SQ23 first.

1. SQ23

Now, I am lazy and don't remember the precise formula to find the appropriate parameterisation of a geodesic passing through two arbitrary latitude and longitude anchor points³. However, we realise that Singapore and New York City have a time difference of 12 hours when both are on Daylight Savings, which is why it is a reasonable approximation that the geodesic passes over the North Pole in the first place. This is obviously not fully correct (and we're deviating from the geodesic anyway), but I will assume that this is the case to greatly simplify our subsequent calculations.

Now, I need to get hold of the latitudes of Singapore and NYC somehow. Unfortunately, Mathematica's location-based functionality all requires Internet connectivity. However, I have quite fortunately set up redshift in two different places, so these numbers are in my redshift.conf file. In particular,

; Configuration of the location-provider:
; type 'redshift -l PROVIDER:help' to see the settings.
; ex: 'redshift -l manual:help'
; Keep in mind that longitudes west of Greenwich (e.g. the Americas)
; are negative numbers.
[manual]
;New Haven
lat=41.3083
lon=-72.9279
;Singapore
;lat=1.3677
;lon=103.7500

This also backs up our intuition about the geodesic being mostly latitudinal. Therefore we are travelling in mostly the latitudinal direction for about 138 degrees (49 deg to the North Pole, then 89 deg down). Assuming we are travelling at a constant speed, we now have latitude and longitude as a function of time.

λ_start = 41 * 2 * np.pi / 360
λ_end = 1 * 2 * np.pi / 360

Θ = np.pi - λ_start - λ_end

def λ_plane(t):
    λ = λ_start + t/T*Θ
    return np.where(λ > np.pi/2, np.pi - λ, λ)

# the longitude of the plane is either the starting longitude or that + π,
# depending on whether it's gone past the North Pole or not.

φ0 = 0

def φ_plane(t):
    λ = λ_start + t/T*Θ
    return np.where(λ > np.pi/2, φ0+np.pi, φ0)

2. The Terminator

Starting to get lazy now. Let's say that the Sun is at a declination of $\delta$ above the equator (which we get from the time of the year). With some algebra we find that the longitude offset $\phi$ between the mean solar noon (which is a meridian) and the terminator at latitude $\lambda$ is

The longitude of the solar noon meridian is itself a function of time, so we can rearrange this to solve for $\lambda$ as a function of time (up to appropriate choices of branch for the arctan).

# Today, Nov 22, is 30 + 31 + 1 days after Sept 21,
# the autumnal equinox.
# The declination of the sun goes as
# sin δ = -sin (23.5 deg) * sin(62 / 365.24 * 2pi)
# (negative because winter)

sin_δ = - np.sin(23.5 * 2 * np.pi / 360) * np.sin(62 / 365.24 * 2 * np.pi)
δ = np.arcsin(sin_δ)

# at t = 0, local time is 2140, meaning it is solar noon
# at a longitude 9.6 hr west (1 hr = 15 deg) of the starting point.
# The solar noon meridian moves west at a rate of 1 hr per hr.

def Φ_sun(t):
    return (φ0-9.6-t)*15/360*2*np.pi

def λ_offset(Δφ):
    return np.arctan(-np.cos(Δφ)/np.tan(δ))

def λ_terminator(t):
    return λ_offset(Φ_sun(t) - φ_plane(t))

Putting it all together

I plot the latitude of our idealised SQ23 and that of the terminator along its current meridian as a function of time.

t = np.linspace(0, T, 500)

plt.plot(t, λ_plane(t), label="Ideal SQ23")
plt.plot(t, λ_terminator(t), label="Terminator geodesic latitude")

plt.ylabel(r"Latitude/rad")
plt.xlabel(r"t/h")

plt.legend()

plt.show()

SQ23 begins its journey north of the terminator at night; we can see that our idealised SQ23 remains in darkness throughout its entire journey. Editing this post now on the ground, I can say that SQ23 actually deviated quite significantly from this idealised trajectory (in particular it passed south of Svalbard!), going around the North Pole in the westerly direction, and therefore briefly experienced daylight. Perhaps this is why the plane landed almost an hour ahead of schedule: it must have benefited quite significantly from tailwinds in the northern polar jet stream.