User:CosineKitty/Calculating planet coordinates from angle measurements

Problem

Suppose you observe a planet (or other celestial object) and you want to determine its equatorial coordinates based on observation alone. One method is to measure the apparent angles between the planet and three or more stars that have known equatorial coordinates. From these measurements, you can calculate the planet's coordinates using the procedure shown here.

Solution strategy

The celestial sphere. The Earth is the point K. The unit vector from the Earth to the planet is u. The unit vectors v and w point to the two stars. The angles measured by the observer are β and γ.

We approach this problem by using angles between the planet and two of the three stars. Imagine the Earth as a point K, and the celestial sphere as a unit sphere centered at K. As seen from K, the set of points on the celestial sphere that are separated from a given star by an angle of β will be a circle. This is the intersection of the unit celestial sphere with an infinitely long cone emanating from K whose angular radius is β.

The set of points on the celestial unit sphere that are an angle γ away from another star will be a different circle. Because we are starting from the assumption that both β and γ are actual angular measurements from two stars to the same planet, we know that the two circles must intersect at the planet's location in the sky. It is possible that the circles are tangent at the planet's location, but in general, there will be a second intersection point also.

Our goal is to find both intersection points, one of which will be the planet's location, the other being a false solution. By applying the same procedure again with the planet, one of the original stars, and a third star, we again obtain a correct solution and another false solution. Unless we are extremely unlucky, it will become apparent at this time which of the solutions is correct: the correct solution will be the same in the first and second attempt, while the two incorrect solutions will not match any of the others.

Another possibility is that the observer can make note of the approximate location of the planet, perhaps drawing a sketch. When a pair of solutions is calculated by the following method, they can both be plotted on a star map. Whichever solution matches the sketch is the correct answer.

Definitions

We concentrate on the planet, two stars, and the angles from the planet to the respective stars. We define the following unit vectors:

u = unit vector from observer to planet

v = unit vector from observer to star #1

w = unit vector from observer to star #2

We also define symbols for the two angles that the observer measured:

β = the angle between u and v

γ = the angle between u and w

The observer must choose stars for which equatorial coordinates are known, thereby making v and w known via calculation. Our goal is to find a pair of solutions for u using the known values of v, w, β, and γ.

Setting up the solution

We look up the equatorial coordinates of the two stars:

α1 = right ascension of star #1	α2 = right ascension of star #2
δ1 = declination of star #1	δ2 = declination of star #2

We transform the equatorial coordinates into unit vectors, and define constants A through F for the components of the unit vectors:

${\begin{aligned}\mathbf {v} &=\left\langle \cos \alpha _{1}\cos \ \delta _{1},\;\sin \alpha _{1}\cos \delta _{1},\;\sin \delta _{1}\right\rangle =\left\langle A,B,C\right\rangle \\\mathbf {w} &=\left\langle \cos \alpha _{2}\cos \ \delta _{2},\;\sin \alpha _{2}\cos \delta _{2},\;\sin \delta _{2}\right\rangle =\left\langle D,E,F\right\rangle \\\end{aligned}}$

We use x, y, z to stand for the unknown components of the unit vector u that we are solving for:

$\mathbf {u} =\left\langle x,y,z\right\rangle$

The dot product of any two vectors is equal to the cosine of the angle between them. We use this fact, along with the constraint that u must be a unit vector, to write a system of three equations with three unknowns:

$x^{2}+y^{2}+z^{2}=1\,$ (1)

$\mathbf {u} \cdot \mathbf {v} =Ax+By+Cz=\cos \beta$ (2)

$\mathbf {u} \cdot \mathbf {w} =Dx+Ey+Fz=\cos \gamma$ (3)

Expressing y and z in terms of x

Multiply equation (2) by F, equation (3) by −C, and add the two results, so as to eliminate z:

$(AF-CD)x+(BF-CE)y=F\cos \beta -C\cos \gamma \,$

Solve for y in terms of x:

$y={\frac {(CD-AF)x+(F\cos \beta -C\cos \gamma )}{BF-CE}}$ (4)

Likewise, we solve for z in terms of x to obtain:

$z={\frac {(AE-BD)x+(B\cos \gamma -E\cos \beta )}{BF-CE}}$ (5)

Equations (4) and (5) require that

$BF\neq CE$

If this constraint is not met, the solution method fails because of a degenerate case in the original angle measurements. For example, if you try to use the same star twice, the two circles are the same circle and overlap in an infinite number of places. Another possibility is that the two stars are exactly on opposite sides of the sky (e.g., on the north and south celestial poles, respectively) and the angles β and γ are both 90°.

For the sake of simplicity, define the following constants:

${\begin{cases}Q=BF-CE\neq 0\\R=CD-AF\\S=AE-BD\\L=F\cos \beta -C\cos \gamma \\M=B\cos \gamma -E\cos \beta \\\end{cases}}$

Now we can rewrite equations (4) and (5) as:

$y={\frac {Rx+L}{Q}}$ (6)

$z={\frac {Sx+M}{Q}}$ (7)

Solving for x using the unit vector constraint

We now have y and z in terms of x, and using equation (1), we can write a quadratic equation in terms of x:

$x^{2}+\left({\frac {Rx+L}{Q}}\right)^{2}+\left({\frac {Sx+M}{Q}}\right)^{2}=1$

Expanding the squared terms and rearranging, we have:

$ax^{2}+bx+c=0\,$

where

${\begin{cases}a=Q^{2}+R^{2}+S^{2}\\b=2(RL+SM)\\c=M^{2}+L^{2}-Q^{2}\end{cases}}$

Now we have two solutions for x, being (in the general case) two distinct x-coordinates where the celestial circles intersect:

$x={\frac {-b\pm {\sqrt {b^{2}-4ac}}}{2a}}\quad :\quad a\neq 0,\quad b^{2}\geq 4ac$

Because a is the sum of three squared real numbers, it must be nonnegative. We already established that Q is not zero, so therefore a must be positive, leaving us safe from worries of dividing by zero. The only other obstacle that could occur here is if the radicand b² − 4ac were negative. If this happens, it means that the celestial circles do not intersect. This can happen only if the observer made a mistake in measuring the planet/star angles in the first place. If the radicand is exactly equal to zero, it means there is only one distinct solution for u because the circles are tangent.

The two solutions for the planet's location

Substitute either of the x values back into equations (6) and (7) to obtain the cartesian coordinates of u:

$\mathbf {u} =\left\langle x,\;{\frac {Rx+L}{Q}},\;{\frac {Sx+M}{Q}}\right\rangle$

At this point it is a good idea to check both solutions for u to make sure they really do satisfy equations (1), (2), and (3).

Convert each of the solutions for u back into equatorial coordinates, using the atan2 function for simplicity of expression:

$\alpha _{p}=\operatorname {atan2} (y,x)$

$\delta _{p}=\arctan \left({\frac {z}{\sqrt {x^{2}+y^{2}}}}\right)$

Eliminating the false solution

As discussed above, the false solution may be eliminated by repeating this method using additional stars, or by comparing the two solutions with an approximate position provided by the observer.

Because of errors in the angle measurement, there will not be exact numerical agreement among all correct solutions. The more stars that are used, the greater is the accuracy that can be achieved by eliminating outliers (excessive errors) and averaging the solutions that remain.

Example

Let's use the method on some actual numbers. We will use the equatorial coordinates for Sirius, Castor, and Regulus, along with the angles these stars formed with Mars on January 1, 2010 at 05:00 UTC.

Star	RA (hours)	DEC (°)	Angle with Mars (°)
Sirius	6.75992	−16.72706	53.757
Castor	7.58732	31.86593	28.883
Regulus	10.14843	11.91800	11.649

Using Sirius and Castor, we find the following pair of solutions for the cartesian coordinates:

u₁ = < 0.153263954977995, 0.826648436983328, 0.541444846441006 >

u₂ = < −0.75065252677078, 0.577055662258427, 0.321756968390946 >

Using Sirius and Regulus, we find:

u₃ = < −0.750651957436305, 0.577055999065544, 0.32175769258793 >

u₄ = < −0.892255001056021, 0.451500291761225, 0.00533850447663287 >

The solutions u₂ and u₃ agree to five places after the decimal point. These are the correct solutions, and u₁ and u₄ are the false solutions. The correct solutions disagree slightly because of the limited precision in the values we started with. In reality, angle measurements β and γ will not be this precise, and there will be less agreement between the correct solutions.

Limiting to the degree of precision implied by the correct solutions, we determine that

u = < −0.750652, 0.577056, 0.321758 >

Converting to equatorial coordinates, we find that the position of Mars was:

RA = 142.44915° = 9^h29^m

DEC = +18.76926° = +18°46'

This result agrees to within one minute of arc with ephemeris values for Mars.