Shadow time for a satellite on an inclined orbit

Let the radius of Earth be $R_e$, and let the radius of the orbit of the satellite be $R \gt R_e$. Choose the Cartesian frame to have its origin at the center of Earth, and with $x$ axis pointing towards the light source. Since the light rays are parallel, they define a cylinder of light whose radius is $R_e$. The equation of this cylinder is

$$ y^2 + z^2 = R_e^2 $$

If $r = [x, y, z]^T $, then the above equation is

$r^T Q r = R_e^2 $

where $Q = \begin{bmatrix} 0 && 0 && 0 \\ 0 && 1 && 0 \\ 0 && 0 && 1 \end{bmatrix}$

Now the orbit of the satellite is a circle centered at the origin (the Earth center) and is spanned by two vectors $u_1, u_2$ which are mutually perpendicular and each having a magnitude of $R$. Thus the parametric equation of the satellite position is

$p(t) = u_1 \cos t + u_2 \sin t = V u $

where $V = [u_1, u_2]$ and $u = [\cos t , \sin t ]^T $

All we need to do is intersect this orbit with the cylinder of shadow and generate the points of intersection. Plugging the equation of the orbit into the equation of the cylinder we obtain,

$ u^T V^T Q V u = R_e^2 $

Let $G = V^T Q V$, then $G$ is a $2 \times 2$ matrix. Expansion of the left hand side results in,

$G_{11} \cos^2 t + G_{22} \sin^2 t + 2 G_{12} \cos t \sin t = R_e^2 $

Using the double angle formulas this reduces to

$ \dfrac{1}{2} \left( G_{11} + G_{22} \right) + \dfrac{1}{2} (G_{11} - G_{22}) \cos 2 t + G_{12} \sin 2 t = R_e^2 $

And this equation can solved using standard methods for solving the equation

$ a \cos \phi + b \sin \phi = c $

After finding the roots of this equation (the one involving $\cos 2 t$ and $\sin 2 t$) which can be up to $4$ roots, one can determine (through evaluating the position of the satellite) the shadow time.