calculate sunset in flight I found an answer at 
https://physics.stackexchange.com/questions/26795/how-to-calculate-the-amount-of-night-time-during-a-flight

but I was stuck at the part:

$\cosθ_1=\sinφ_1\sinφ+\cosφ_1\cosφ\cos(λ−λ_1)$

$\cosθ_2=\sinφ_2\sinφ+\cosφ_2\cosφ\cos(λ−λ_2)$

from which  $ (λ(t),φ(t))$

can be derived (after some tedious calculations).
After choosing an arbitrary time t, you have the two equations you reference and two unknowns: $ λ$ and $φ$. The problem is that the equations cannot be solved directly, so an iterative approach is required. For example, 


assume that $(λ−λ_1)=θ_1$. Note that this also lets you calculate $λ$
solve for $φ$ from the first equation
use that value of $φ$ in the second equation to solve for $λ$
If the values of $λ$ from equation 1 and 2 are different, try another value for $λ$ in the first equation.
repeat until the assumed value for $λ$ in equation 1 and the calculated value from equation 2 are equal.


It has been a long time since I did spherical trig, but I think the solution on the Physics SE is inefficient. It looks like they are calculating the sun's altitude at multiple point (multiple times of t) along the path, so the two equations two unknowns needs to be solved multiple times -- that is inefficient. It should be possible to determine the equation of the great circle through the departure and destination point. Essentially you need to calculate the longitude of where it cross the equator and the inclination. That probably needs to be solved iteratively, but it only needs to be solved once. Then for any time t, the position $λ$ and $φ$ can be calculated directly, and the sun's altitude can be calculated directly.