# How to calculate the horizon line of a satellite?

I need an equation to calculate a list of Earth-centered, Earth-fixed (ECEF) XYZ coordinates on the earth that represent the visibility limit of satellite given its ECEF XYZ coordinates.

For any given point on the earth I can tell if it is in view. I can calculate the distance to the horizon from the vehicle and the off nadir angle that points at the horizon, but how do I deterministically calculate the actual "Horizon Line" of the satellite? (I could iterate to find points that have a 0 degree elevation angle on the earth, but that takes too long.)

Here is a simple diagram of what I have so far.

• What do you mean a list of points? It is a continuous ring of points, so a function would be more appropriate. – Jim Dec 9 '14 at 14:55
• Yes, I am looking for a function. I actually need to draw them on a map, so I will need to calculate discrete points for line segment drawing – mlbrink Dec 9 '14 at 15:00
• So then you should be able to use pencil and paper to figure out an appropriate function based on the ECEF XYZ coordinates of the satellite (I've done something similar so I know it's possible). Then you can plug that equation into MatLab or something similar. It doesn't take long at all – Jim Dec 9 '14 at 15:02
• I apologize. I am not being patronizing but it evidently came across that way. I said I have done something similar but not this exactly. When you said you can iterate to find points with 0 degrees elevation, I assumed that meant you could thus find an equation for it and convert those points back to ECEF. Since this is virtually all that need be done, it was my understanding that you had all you need and were stressing the details of it. – Jim Dec 9 '14 at 15:11
• So I have equations for calculating the elevation angle to the satellite from a point on the earth. This doesn't really help since I don't want to iterate over a large number of points to find the actual horizon equation. – mlbrink Dec 9 '14 at 15:15

I'm going to assume you have some familiarity with linear algebra, as the math becomes much less tedious than trying to do three-dimensional trigonometry with the $x$, $y$, $z$-coordinates directly.

You are looking for a function describing the line of the horizon. Since it is a circle and thus a one-dimensional object, I'm going to call it $\vec{h}(t)$, where $t$ is a parameter variable such that $0 \leq t < 2\pi$. What I'm going to do is find the transformations that turn this horizon circle into a unit circle on the x-y plane centered at the origin (which I already know how to do), then apply those transformations in reverse.

First, define $\vec{s}$ as the vector indicating the position of the satellite with respect to the center of the Earth--the line of length $R+h$ in your diagram). It has components $$\vec{s} = \begin{bmatrix}s_x\\s_y\\s_z\end{bmatrix}.$$The horizon circle $\vec{h}(t)$ is centered on the line from the center of the Earth ($\vec{0}$ in ECEF) to the satellite. This line is also perpendicular to the plane that contains the circle.

Imagine a line perpendicular to $\vec{s}$ that connects the horizon point. This line is a radius of the horizon circle, and where this line meets $\vec{s}$ is the center of the circle $\vec{h}(t)$. I'm going to call this point $\vec{h_c}$. To find out this quantity, I'm going to define the angle $\beta$ as the angle between the horizon point and the satellite, centered at Earth's center. It's equal to $90^o-\alpha$ ($\alpha$ should be given by $\sin^{-1}(d/(R+h))$) and calculated by $$\beta = \cos^{-1}\left(\frac{R}{R+h}\right) = \cos^{-1}\left(\frac{R}{\left|\vec{s}\right|}\right).$$ See the picture below.

The distance from Earth's center to $\vec{h_c}$ is $$\left| \vec{h_c} \right| = R \cos\beta.$$ So, the position of the center of $\vec{h}(t)$ is $$\vec{h_c} = \left(R \cos\beta \right)\hat{s},$$ where $\hat{s}$ is the unit vector in the direction of $\vec{s}$. We can also immediately see that the radius of the horizon circle is $R \sin\beta$.

Now, we can turn the unknown $\vec{h}(t)$ into something simple. $$\vec{h}(t) - \vec{h_c}$$ is the horizon circle moved so that it is centered at the center of the Earth. $$\frac{\vec{h}(t) - \vec{h_c}}{R\sin\beta}$$ shrinks the horizon line so it has a radius of one.

Now, we need a transformation that rotates this unit circle onto the $x$-$y$ plane so we can write that circle as $$\begin{bmatrix}\cos{}t\\\sin{}t\\0\end{bmatrix},$$ where this $t$ is the same as $\vec{h}(t)$. In other words, we need an orthogonal transformation $T$ such that $$T \begin{bmatrix}0\\0\\1\end{bmatrix} = \hat{s}$$ I'm going to cheat a little bit here and refer to another of my answers where I derive one possible $T$. $$T = \begin{bmatrix} \frac{-\hat{s}_y}{\sqrt{\hat{s}_x^2 + \hat{s}_y^2}} & \frac{-\hat{s}_x \hat{s}_z}{\sqrt{\hat{s}_x^2 + \hat{s}_y^2}} & \hat{s}_x\\ \frac{\hat{s}_x}{\sqrt{\hat{s}_x^2 + \hat{s}_y^2}} & \frac{-\hat{s}_y \hat{s}_z}{\sqrt{\hat{s}_x^2 + \hat{s}_y^2}} & \hat{s}_y\\ 0 & \sqrt{\hat{s}_x^2 + \hat{s}_y^2} & \hat{s}_z \end{bmatrix}$$ (Notice that this is actually $T^{-1}$ in the other answer since that question wanted the opposite transformation.) The variables $\hat{s}_{x,y,z}$ are the components of the unit vector $\hat{s}$.

Now, we have $$\frac{\vec{h}(t) - \vec{h_c}}{R\sin\beta} = T\begin{bmatrix}\cos{}t\\\sin{}t\\0\end{bmatrix}.$$ Solving for $\vec{h}(t)$, $$\vec{h}(t) = \vec{h_c} + \left( R \sin\beta \right) T\begin{bmatrix}\cos{}t\\\sin{}t\\0\end{bmatrix}.$$ With everything substituted in and multiplied together, I get $$\vec{h}(t) = R \left( \frac{R}{R+h} \begin{bmatrix}\hat{s}_x\\\hat{s}_y\\\hat{s}_z \end{bmatrix} + \sqrt{\frac{1-\left(\frac{R}{R+h}\right)^2}{\hat{s}_x^2 + \hat{s}_y^2}} \begin{bmatrix} -\hat{s}_y \cos{}t - \hat{s}_x \hat{s}_z \sin{}t \\ \hat{s}_x \cos{}t - \hat{s}_y \hat{s}_z \sin{}t \\ \left( \hat{s}_x^2 + \hat{s}_y^2 \right) \sin{}t \end{bmatrix} \right)$$ or $$\vec{h}(t) = R \left( \frac{R}{\left|\vec{s}\right|} \begin{bmatrix}\hat{s}_x\\\hat{s}_y\\\hat{s}_z \end{bmatrix} + \sqrt{\frac{1-\left(\frac{R}{\left|\vec{s}\right|}\right)^2}{\hat{s}_x^2 + \hat{s}_y^2}} \begin{bmatrix} -\hat{s}_y \cos{}t - \hat{s}_x \hat{s}_z \sin{}t \\ \hat{s}_x \cos{}t - \hat{s}_y \hat{s}_z \sin{}t \\ \left( \hat{s}_x^2 + \hat{s}_y^2 \right) \sin{}t \end{bmatrix} \right)$$