Exercises — Orbit determination — Gauss's method, Gibbs method

Method: Gibbs. Gibbs's input is exactly three position vectors, and it needs no time information (see the parent-note mistake box "I need the times for Gibbs"). The full state of an orbiting body is six numbers: position $\vec{r}$ (which you already have at the middle point) plus velocity $\vec{v}$. So the single vector you are missing is the velocity ${\vec{v}}_2$ at the middle point. Gibbs's whole job is to produce ${\vec{v}}_2$.

• (a) → Gauss. Angles alone give direction ${\hat{\rho }}_i$ but not distance (the slant range ${\rho }_i$). Gauss's cleverness is finding those unknown ranges so that ${\vec{r}}_i={\vec{R}}_i+{\rho }_i{\hat{\rho }}_i$ becomes fully known.
• (b) → Gibbs. The positions are already known, so you skip straight to the velocity-recovery step. Pipeline mnemonic: Angles → Gauss → Positions → Gibbs → Velocity.

The root is $r_2$, the magnitude of the middle position vector (its distance from Earth's centre). Pick the root that is (1) positive and real (a distance cannot be negative or imaginary) and (2) larger than the central body's radius — for Earth, larger than $R_{\oplus }\approx 6378\text{km}$ — and that yields positive slant ranges ${\rho }_i>0$ (the satellite must be above the surface and in front of the telescope, not behind it). Complex, negative, or sub-surface roots (any $r_2<6378\text{km}$) are unphysical and discarded.

The magnitude is $r_i=\sqrt{r_{ix}^2+r_{iy}^2+r_{iz}^2}$ (a scalar length). $r_1=\sqrt{294.32^2+4265.1^2+5986.7^2}=7356.5\text{km}$ $r_2=\sqrt{1365.5^2+3637.6^2+6346.8^2}=7440.3\text{km}$ $r_3=\sqrt{2940.3^2+2473.7^2+6555.8^2}=7598.9\text{km}$ All three exceed $R_{\oplus }=6378$ km ✓ (these are genuine above-surface positions). These scalars are the weights that appear inside $\vec{N}$ and $\vec{S}$, so they must come first.

Why this test: ${\vec{r}}_2\times {\vec{r}}_3$ is a vector perpendicular to the plane containing ${\vec{r}}_2$ and ${\vec{r}}_3$. If ${\vec{r}}_1$ lies in that same plane, it is perpendicular to that normal, so the dot product is (near) zero. See the figure below.

Figure s01 — coplanarity test. Axes are arbitrary in-plane coordinates $x,y$. The mint parallelogram is the orbit plane; the three coloured arrows (lavender ${\vec{r}}_1$, coral ${\vec{r}}_2$, butter ${\vec{r}}_3$) all start at the focus and lie in that plane; the slate arrow is the normal ${\vec{r}}_2\times {\vec{r}}_3$. Coplanarity = all three coloured arrows perpendicular to the slate one. Getting the numbers straight (units matter). The cross product of two vectors in km has units km·km = km². Carrying out ${\vec{r}}_2\times {\vec{r}}_3$ component by component: ${\vec{r}}_2\times {\vec{r}}_3\approx (8.1454\times 10^6,1.9095\times 10^6,7.3908\times 10^6){\text{km}}^2.$ Dotting with ${\vec{r}}_1$ (km) gives a triple product in km³. Each of its three product-terms is of order $r_1\cdot |{\vec{r}}_2\times {\vec{r}}_3|\sim 7.4\times 10^3\cdot 1.1\times 10^7\approx 8\times 10^{10}{\text{km}}^3$, yet they nearly cancel: ${\vec{r}}_1\cdot ({\vec{r}}_2\times {\vec{r}}_3)\approx -8.0\times 10^4{\text{km}}^3.$ So the relative size of the leftover is $\frac{8\times 10^4}{8\times 10^{10}}\approx 1\times 10^{-6}$ — one part in a million. The clean way to read this is the angle test on unit vectors: ${\hat{u}}_{r1}\cdot ({\hat{u}}_{r2}\times {\hat{u}}_{r3})\approx 1.3\times 10^{-6}$, whose arcsine is an out-of-plane angle $\ll 1^{\circ }$. Coplanar ✓ — a single orbit can thread all three.

$\vec{N}=r_1({\vec{r}}_2\times {\vec{r}}_3)+r_2({\vec{r}}_3\times {\vec{r}}_1)+r_3({\vec{r}}_1\times {\vec{r}}_2)$ $\vec{D}=({\vec{r}}_1\times {\vec{r}}_2)+({\vec{r}}_2\times {\vec{r}}_3)+({\vec{r}}_3\times {\vec{r}}_1)$ $\vec{S}={\vec{r}}_1(r_2-r_3)+{\vec{r}}_2(r_3-r_1)+{\vec{r}}_3(r_1-r_2)$ (Here each $r_i$ without an arrow is the scalar magnitude weight.) Carrying out the arithmetic: $\vec{N}\approx (6.0983\times 10^{10},1.4262\times 10^{10},5.5333\times 10^{10}){\text{km}}^3,N\approx 8.3465\times 10^{10}$ $\vec{D}\approx (8.1971\times 10^6,1.9169\times 10^6,7.4370\times 10^6){\text{km}}^2,D\approx 1.1215\times 10^7$ $\vec{S}\approx (6708.3,-10516,-2516.2){\text{km}}^2$ Semi-latus rectum $p=N/D\approx 7442\text{km}$. (That $p$ is close to $r_2$ tells you point 2 is near the semi-latus-rectum crossing — a good sanity feel, not a requirement.)

The prefactor: $\sqrt{\mu /(ND)}=\sqrt{398600/(8.3465\times 10^{10}\cdot 1.1215\times 10^7)}\approx 2.066\times 10^{-7}$. The bracket has two pieces — the transverse part $\vec{D}\times {\vec{r}}_2/r_2$ (velocity along the motion) and the radial part $\vec{S}$ (velocity changing the distance). Adding and scaling: $\boxed{{\vec{v}}_2\approx (-6.218,-4.014,1.599)\text{km/s}}$ matching the parent note. Now $({\vec{r}}_2,{\vec{v}}_2)$ is the complete state.

Setup — write everything in the orbit plane. Put a 2-D polar frame in the orbit plane with the focus at the origin and the polar axis along perigee. A point at true anomaly ${\theta }_i$ then has the vector ${\vec{r}}_i=r_i(\cos {\theta }_i\hat{x}+\sin {\theta }_i\hat{y}),\frac{1}{r_i}=\frac{1+e\cos {\theta }_i}{p}.$ Cross products become sines. Any two of these vectors satisfy ${\vec{r}}_i\times {\vec{r}}_j=r_i{r}_j\sin ({\theta }_j-{\theta }_i)\hat{z},$ because the cross product of two planar vectors is (product of lengths)×(sine of the angle between them) along the normal $\hat{z}$. Write the three anomaly gaps as ${\theta }_{21}={\theta }_2-{\theta }_1$, ${\theta }_{32}={\theta }_3-{\theta }_2$, ${\theta }_{31}={\theta }_3-{\theta }_1$. The two bracket pieces at index 2. Using the above, the transverse term is $\frac{\vec{D}\times {\vec{r}}_2}{r_2}=\frac{1}{r_2}[({\vec{r}}_1\times {\vec{r}}_2)+({\vec{r}}_2\times {\vec{r}}_3)+({\vec{r}}_3\times {\vec{r}}_1)]\times {\vec{r}}_2.$ Each inner cross product is an $\hat{z}$ vector $\propto \sin (\cdot )$; crossing $\hat{z}$ with ${\vec{r}}_2$ rotates ${\vec{r}}_2$ by $90^{\circ }$ inside the plane, so this term is a planar vector whose $\cos {\theta }_2$ and $\sin {\theta }_2$ components carry the tag "${\theta }_2$". Meanwhile $\vec{S}={\vec{r}}_1(r_2-r_3)+{\vec{r}}_2(r_3-r_1)+{\vec{r}}_3(r_1-r_2)$ also expands into $\cos {\theta }_i,\sin {\theta }_i$ pieces. The cancellation: substitute $r_i=p/(1+e\cos {\theta }_i)$ into both. The ${\theta }_1$ and ${\theta }_3$ tags appear in $\vec{S}$ multiplied by factors $\propto \sin ({\theta }_{32})$ and $\sin ({\theta }_{21})$ that are exactly the same sine factors produced by the transverse term at index 2 — but with the opposite sign (one comes from ${\vec{r}}_1\times {\vec{r}}_2$, the other from ${\vec{r}}_2\times {\vec{r}}_3$, which differ in orientation). Adding the two brackets, the ${\theta }_1$-tagged pieces annihilate the ${\theta }_1$-tagged pieces and likewise for ${\theta }_3$; only terms tagged with the middle anomaly ${\theta }_2$, together with the pure constants $p$ and $e$, survive. The result is a clean multiple of $\sqrt{\mu p}$ along the transverse and radial directions — i.e. exactly ${\vec{v}}_2$. Why index 1 fails. Repeat with $\vec{D}\times {\vec{r}}_1/r_1$ instead. Now the transverse term is tagged "${\theta }_1$", but $\vec{S}$ still supplies ${\theta }_1,{\theta }_2,{\theta }_3$ pieces. The would-be partner for the ${\theta }_2$ tag is missing (there is no $\vec{D}\times {\vec{r}}_1$ term that produces a matching $\sin ({\theta }_2-{\theta }_2)=0$), so a residual $\propto \sin ({\theta }_2-{\theta }_1)\ne 0$ survives and pollutes the answer. The symmetric sums are literally centred so the books balance only at the middle index.

Figure s02 — velocity decomposition. Axes are arbitrary in-plane $x,y$. The lavender curve is the orbit, the slate dot the focus, the coral arrow ${\vec{r}}_2$ the radius to the middle point. At the tip, the mint arrow is the transverse direction (perpendicular to ${\vec{r}}_2$, along the motion) from $\vec{D}\times {\vec{r}}_2/r_2$; the butter arrow is the radial direction (along ${\vec{r}}_2$) from $\vec{S}$; the slate arrow is their sum, the true velocity ${\vec{v}}_2$. On any conic, the velocity at a point decomposes into two perpendicular directions:

• Transverse (perpendicular to the radius, along the motion) — supplied by $\frac{\vec{D}\times {\vec{r}}_2}{r_2}$. Since $\vec{D}$ is (near) the orbit-normal direction, crossing it with ${\vec{r}}_2$ rotates ${\vec{r}}_2$ by $90^{\circ }$ within the plane, giving the along-track direction.
• Radial (along ${\vec{r}}_2$, changing the distance) — supplied by $\vec{S}$, which carries the $\dot{r}$ part. Velocity in a plane needs exactly two independent directions; radial + transverse is the natural pair for orbits, so the sum is complete — nothing is missing.

Non-coplanar case. Gibbs assumes a single planar conic threads all three points. If they are genuinely non-coplanar, no such orbit exists, so the formula still returns a number for ${\vec{v}}_2$ but that number is meaningless: the reconstructed $\vec{h}$ and $\vec{e}$ don't correspond to any real orbit through the points. Symptoms: the resulting orbit misses points 1 and 3 badly, and derived elements (eccentricity, inclination) come out inconsistent with a re-propagation. That is why the $<1^{\circ }$ gate is a hard prerequisite. Collinear (singular) case — the other failure mode. Look at $\vec{D}=({\vec{r}}_1\times {\vec{r}}_2)+({\vec{r}}_2\times {\vec{r}}_3)+({\vec{r}}_3\times {\vec{r}}_1)$. Each cross product is proportional to the sine of the angle between two position vectors. If the three points are nearly collinear (all bunched along one line — a very short arc with tiny angular separation between the radius vectors), every cross product shrinks toward zero, so $\vec{D}\to \vec{0}$ and $\vec{N}\to \vec{0}$. The prefactor $\sqrt{\mu /(ND)}$ then blows up ($1/\sqrt{ND}\to \infty$) while the bracket $\to \vec{0}$ — a $0\times \infty$ indeterminate form that is catastrophically noise-sensitive: any measurement jitter in the tiny cross products gets amplified enormously, so the recovered ${\vec{v}}_2$ is dominated by noise rather than signal. Rule of thumb: Gibbs needs the three points to be spread out in true anomaly (well-separated angles). For short, tightly-spaced arcs, prefer Gauss's series treatment or gather more widely-spaced observations.

$\vec{h}={\vec{r}}_2\times {\vec{v}}_2\approx (3.1310\times 10^4,-3.7660\times 10^4,2.7106\times 10^4){\text{km}}^2/\text{s}$ $h=|\vec{h}|\approx 5.5757\times 10^4{\text{km}}^2/\text{s}$ The conservation check done correctly. The authoritative semi-latus rectum comes from the state via $p=h^2/\mu$: $p=\frac{h^2}{\mu }=\frac{(5.5757\times 10^4{)}^2}{398600}\approx 7800\text{km}.$ Then $\sqrt{\mu p}=\sqrt{398600\cdot 7800}=\sqrt{3.109\times 10^9}\approx 5.576\times 10^4=h$ ✓ — the relation $h=\sqrt{\mu p}$ closes to the last digit because both sides now use the same $p$. Why the earlier $N/D\approx 7442$ differed: that $p$ was assembled from hand-carried, rounded helper vectors $\vec{N},\vec{D}$; small rounding in each component propagates into $N/D$. The lesson: once you have $({\vec{r}}_2,{\vec{v}}_2)$, derive conserved quantities from the state, which is internally consistent, rather than from intermediate scratch values.

With ${\hat{r}}_2={\vec{r}}_2/7440.3\approx (-0.18354,0.48889,0.85304)$ and $\vec{h}$ from 4.1: $\vec{e}=\frac{{\vec{v}}_2\times \vec{h}}{\mu }-{\hat{r}}_2\approx (-0.0868,-0.0663,0.0195)$ $e=|\vec{e}|\approx 0.111$ Because $0<e<1$, the orbit is an ellipse (bound). $e\approx 0.11$ is mildly elliptical — nearly circular but clearly not a perfect circle. The direction of $\vec{e}$ points toward perigee (closest approach), fixing the orbit's orientation inside its plane.

flowchart TD
  A["Three angle pairs alpha delta plus observer R_i"] --> B["Line of sight unit vectors rho hat"]
  B --> C["Gauss f and g series about middle time"]
  C --> D["Eighth degree polynomial in r_2"]
  D --> E["Pick positive real root above Earth radius"]
  E --> F["Back out slant ranges rho_i positive"]
  F --> G["Positions r_i equals R_i plus rho_i rho hat"]
  G --> H["Gibbs helper vectors N D S"]
  H --> I["Velocity v_2 from state"]
  I --> J["State r_2 and v_2"]
  J --> K["Compute h and e vectors"]
  K --> L["Six classical elements a e i RAAN arg perigee"]
  L --> M["True anomaly then eccentric then mean anomaly"]
  M --> N["Kepler equation for time since perigee"]
  E --> P["Refine f and g exact then re solve"]
  P --> D

Quadrant / sign care (the part beginners skip):

• Inclination $i=\arccos (h_z/h)$ lands in $[0^{\circ },180^{\circ }]$ cleanly — $\arccos$ already covers it, no ambiguity.
• RAAN $\Omega$ and argument of perigee $\omega$ come from $\arccos$ of a ratio, which only returns $[0^{\circ },180^{\circ }]$. Resolve the other half using a sign test: if the node vector's $y$-component $n_y<0$ then $\Omega =360^{\circ }-\arccos (\dots )$; if $e_z<0$ then $\omega =360^{\circ }-\arccos (\dots )$.
• True anomaly $\theta$: if $\vec{r}\cdot \vec{v}<0$ the body is moving inward (past apogee, heading to perigee), so $\theta =360^{\circ }-\arccos (\dots )$.
• Kepler's equation $M=E-e\sin E$ then relates eccentric anomaly $E$ to mean anomaly $M$ (hence time). Solve for $E$ from $M$ numerically (Newton's method), keeping $E$ in the same revolution as $\theta$. The iterate loop (right branch) is the key subtlety: Gauss's first ${\rho }_i$ come from truncated Lagrange series, so refine them with exact Lagrange Coefficients (f and g) and re-solve the polynomial until ${\rho }_2$ stops moving.

The scalar triple product is the $3\times 3$ determinant of the three vectors stacked as rows — a standard identity. Expand ${\vec{r}}_2\times {\vec{r}}_3$ into its components and dot with ${\vec{r}}_1$; the resulting sum of $3!=6$ signed products is exactly the cofactor expansion of that determinant. Geometry: $|{\vec{r}}_1\cdot ({\vec{r}}_2\times {\vec{r}}_3)|$ is the volume of the parallelepiped spanned by the three vectors. Three vectors lie in a common plane iff that box is flat — zero volume. Hence determinant $=0\iff$ coplanar. This is why the test is exact in principle and near-zero (relative to term size) in noisy practice.

Units: $\vec{N}$ is (km)(km²)=km³, so $N$ is km³. $\vec{D}$ is km², so $D$ is km². Thus $p=N/D$ has units ${\text{km}}^3/{\text{km}}^2=\text{km}$ ✓, and $p\approx 8.35\times 10^{10}/1.12\times 10^7\approx 7.4\times 10^3$ km. Prefactor: $\mu /(ND)$ has units $\frac{{\text{km}}^3{\text{s}}^{-2}}{{\text{km}}^3\cdot {\text{km}}^2}={\text{km}}^{-2}{\text{s}}^{-2}$, so $\sqrt{\mu /(ND)}$ is ${\text{km}}^{-1}{\text{s}}^{-1}$ ✓. The bracket $\frac{\vec{D}\times {\vec{r}}_2}{r_2}+\vec{S}$ has units $\frac{{\text{km}}^2\cdot \text{km}}{\text{km}}={\text{km}}^2$. Multiply: ${\text{km}}^{-1}{\text{s}}^{-1}\cdot {\text{km}}^2=\text{km}{\text{s}}^{-1}$ ✓. Dimensional analysis alone forces the formula's shape — a powerful sanity net.

If ${\vec{r}}_1={\vec{r}}_2$, then several building blocks collapse: ${\vec{r}}_1\times {\vec{r}}_2=\vec{0}$, and the pair contributes nothing new to $\vec{D}$ or $\vec{N}$; the sums degenerate as if you had only two distinct points. Two points on a conic do not determine it uniquely (infinitely many conics of different eccentricities pass through two given points sharing one focus), so the reconstructed ${\vec{v}}_2$ is undetermined/ill-conditioned. Physically: you need three genuinely different samples of the curve to pin the scale ($p$), the plane ($\vec{h}$ direction), and the shape/orientation ($\vec{e}$). This is the same $\vec{D}\to \vec{0}$ pathology as the collinear case in Exercise 3.3 — coincident points are just the extreme of "not spread out": when the angular spread between samples vanishes, the cross-product sums vanish and $\sqrt{\mu /(ND)}$ diverges.