I need to calculate the beginning and end points of a line that originates at a point in 3D space and then skims the surface of a sphere at a distance location — the tangent line, of course. I know how to do this in two-dimensions, but am stumped when it comes to three:
Of course, I realize there are an infinite number of such lines, since the point of contact with the sphere could be anywhere in the circular cross-section. So I suppose I need a solution that allows includes a theta value representing the place on the sphere that I want the line to be tangent to.
(The purpose of this exercise is to calculate occultation in space between the sun and the moon, treating the sun as a point instead of a sphere.)
Best Answer
The tangent lines to the sphere form a cone with apex at the given point $\mathbf p$ and axis the line through $\mathbf p$ and the center $\mathbf c$ of the sphere. (To keep things clearer, I’ll mark inhomogeneous Cartesian coordinate vectors with a tilde; if there’s no tilde, this indicates a homogeneous coordinate vector.) The points of tangency lie on the intersection of the sphere with the polar plane of $\mathbf p$. This will be a circle in this polar plane, which is perpendicular to the cone’s axis, with center $\mathbf c'$ at the intersection of the plane and axis. To parameterize this circle, find the radius $r'$ via the Pythagorean theorem, and find a pair of unit vectors $\mathbf u$ and $\mathbf v$ that are orthogonal to the axis and each other. The circle is then $$\mathbf q = \mathbf c'+r'\mathbf u\cos\theta+r'\mathbf v\sin\theta.$$ The line through this point be found in any of the usual ways, such as the parametrization $(1-\lambda)\mathbf p+\lambda\mathbf q$.
To keep things clearer in the following, I’ll mark variables that represent inhomogeneous Cartesian coordinate vectors with a tilde; no tilde indicates a homogeneous coordinate vector. Representing the sphere by the matrix $$S = \left[\begin{array}{c|c} I_3 & -\tilde{\mathbf c} \\ \hline -\tilde{\mathbf c}^T & \|\tilde{\mathbf c}\|^2-r^2\end{array}\right],$$ the polar plane to $\mathbf p$ is $\mathbf\pi=S\mathbf p$. (The components of this vectors correspond to the coefficients of the implicit Cartesian equation of the plane.) Use your favorite method to find its intersection with the line through $\mathbf p$ and $\mathbf c$. Here, I’ll use the Plücker matrix of the line: $$\mathbf c' = (\mathbf p\mathbf c^T-\mathbf c\mathbf p^T)\mathbf\pi = (\mathbf c^T\mathbf\pi)\mathbf p - (\mathbf p^T\mathbf\pi)\mathbf c.$$ (The parenthesized quantities on the r.h.s. are just the dot products of $\mathbf p$ and $\mathbf c$ with $\mathbf\pi$.) The radius of the circle is found from $$r'^2 = r^2-\|\tilde{\mathbf c}-\tilde{\mathbf c}'\|^2.$$ For the unit vectors $\mathbf u$ and $\mathbf v$ of the parametrization, choose one of the cross products of $\tilde{\mathbf c}-\tilde{\mathbf p}$ with the unit basis vectors (at least two of them are guaranteed to be nonzero) and normalize it for $\mathbf u$, and then take $\mathbf u\times(\tilde{\mathbf c}-\tilde{\mathbf p})$, normalized, for $\mathbf v$. You can instead use $\tilde{\mathbf c}-\tilde{\mathbf c}'$ or $\tilde{\mathbf c}'-\tilde{\mathbf p}$ for this calculation; I’m not sure which will provide the best numerical stability.
For your example, I get $$S = \begin{bmatrix} 1&0&0&-100\\0&1&0&50\\0&0&1&-75\\-100&50&-75&17500\end{bmatrix}$$ from which $\mathbf\pi=[400:150:-75:-27500]$. Using the above formula for the axis-plane intersection gives $$\mathbf c' = [400:150:-75:-27500],$$ which corresponds to the inhomogeneous Cartesian coordinates $$\tilde{\mathbf c}' = \left({30500\over301},{14900\over301},{22500\over301}\right)$$ and the circle’s radius works out to be $r'=250\sqrt{\frac3{301}}$.
For the unit vector $\tilde{\mathbf u}$, I'll take the cross product that has the greatest norm, $(-150,400,0)^T$, so $\tilde{\mathbf u} = \frac1{\sqrt{73}}(-3,8,0)$ and (according to Mathematica; the radicals are starting to get a bit ugly) $\tilde{\mathbf v}=\left({24\over\sqrt{21973}},{9\over\sqrt{21973}},2\sqrt{{73\over301}}\right)^T$. The parametrization of the circle finally comes out to be approximately $$x = 101.329 - 8.76349 \cos\theta + 4.04095 \sin\theta \\ y= -49.5017 + 23.3693 \cos\theta + 1.51536 \sin\theta \\ z = 74.7508 + 24.5825 \sin\theta.$$
Depending on what it is you want to do with this, it might be more convenient to work with the tangent cone as a whole. It turns out that its matrix is easily computed from $\mathbf p$ and $C$: $$C_{cone} = (\mathbf p^TC\mathbf p)C - (C\mathbf p)(C\mathbf p)^T.$$ Note that the second term isn’t a dot product; it’s the tensor product of $C\mathbf p$ with itself.