Niko's Project Corner

Approximating planets' orbits in closed-form

Description A sufficiently accurate and efficient model for planetary orbits
Languages Matlab
Tags As­tron­omy
Ap­plied math­emat­ics
Duration Fall 2014
Modified 12th October 2014

I wanted to find or cre­ate a for­mula which would ac­cept an epoch times­tamp, lat­itude and lon­gi­tude and it would pro­duce the Sun's ob­served az­imuth and al­ti­tude in ra­di­ans. It needs to take into ac­count de­tails earth's ax­ial tilt and its po­si­tion on its or­bit around the sun. To my sur­prise I wasn't able to find such for­mula, so I had to de­velop it from scratch. Luck­ily earth's or­bit (and or­bits in gen­eral) is a well stud­ied and doc­umented prob­lem, so I could take some short­cuts.

Figure 1: An il­lus­tra­tion of earth's or­bit and its re­la­tion to cal­en­dar day (from Wiki­me­dia com­mons).
Ke­pler's equa­tion does not have an an­alyt­ical so­lu­tion, but it is easy ap­ply in the "op­po­site" di­rec­tion. Here an ap­prox­imate so­lu­tion is cre­ated. First an even sam­pling of true anoma­lies is gen­er­ated, cov­er­ing the whole 360° or­bit ro­ta­tion:
    θ = 0 ... 2 π(1)

Next cor­re­spond­ing ec­cen­tric anomaly E is cal­cu­lated from θ and the or­bit's ec­cen­tric­ity e:

    E = atan2(√1 - e  · sin(θ / 2), √1 + e  · cos(θ / 2))(2)

Mean anomaly M is cal­cu­lated from E by us­ing Ke­pler's equa­tion:

    M = E - e sin E(3)

The M is triv­ial to cal­cu­late for any de­sired cal­en­dar date, so we want to have an ap­prox­ima­tion for­mula θ = fe(M) so that pre­vi­ous equa­tion hold rea­son­ably well. I no­ticed ex­per­imen­tally that it would be a good ap­proach to ap­prox­imate

    f'e(M) = (M - θ) / e(4)
for e > 0, and has a some­what con­stant shape for small val­ues of e. Then the fol­low­ing holds:
    θ = fe(M) = M - e · f'e(M)(5)

The ef­fect of e on f'e(M) shape can be seen in Fig­ure 2. Most im­por­tantly its mag­ni­tude is around ±2, and has val­ues of ex­actly zero at M = 0, M = π and M = 2π. Due to these prop­er­ties, it is ap­prox­imated by the poly­no­mial

    f'e(M) → a(e) · (M - pi) · (M · (2π - M))-b(e)(6)

The shape of f'e(M) for var­ious val­ues of e can be seen in Fig­ure 2. It is no­table that for M → π / 2 it has a value of about -2, which in­di­cates that θ > M around that re­gion. This oc­curs be­cause when M → 0 the planet is clos­est to the Sun, which means it has a higher an­gu­lar ve­loc­ity. Thus θ grows ini­tially faster than the mean anomaly M.

Figure 2: A plot of f'e(M) for var­ious val­ues of e. X-axis has the value of M from 0 to 2π.

The ec­cen­tric­ity e's role on a(e) and b(e) is solved by find­ing op­ti­mal a and b for var­ious dif­fer­ent val­ues of e and fit­ting a model to those find­ings. Pa­ram­eters are es­ti­mated by min­imiz­ing the stan­dard sum of squared er­rors, which is achieved by Mat­lab's lsqnonin func­tion. These re­sults are shown in Fig­ure 3. Luck­ily both pa­ram­eters seem to fol­low the same for­mula y(e) = y0 + exp(k1· log(e) + k2), and k1 and k2 can be reliably estimated by fitting a linear trend to the log-log data points. The constant y0 (which corresponds to a0 and b0, one for both parameters) is estimated so that the best log-linear model can be obtained. If that value is off then the log-log plot is no longer linear, and a more complex model is needed.

Figure 3: Es­ti­mat­ing a(e) and b(e) from sam­pled mod­els.

Model ac­cu­racy is de­ter­mined by cal­cu­lat­ing the max­imum er­ror from the true lo­ca­tion within the whole 360° or­bit. If e = 0 then the or­bit would be ex­actly cir­cu­lar and θ = M. In­ter­est­ing cases are small val­ues of e, for ex­am­ple on earth's or­bit e → 0.0167. Two dif­fer­ent mod­els are com­pared in fig­ures 5 and 6: a sim­pler model which only uses con­stants a0 and b0 regardless of the value of e, and the more accurate model with e-dependent corrections.

The biggest ben­efit of these cor­rec­tions oc­cur when 0.05 < e < 0.3. When e → 0 these corrections are small anyway, and on bigger values of e the model loses its accuracy even with these corrections. The ratio of simple model's maximum error and better model's error is show in Figure 4. When e = 0.1 the maximum ratio of 4 is reached, which means that the simpler model has 300% higher maximum error than the e-corrected version.

Figure 4: The ra­tio of sim­pler model's er­ror to e-cor­rected ver­sion's er­ror for vary­ing val­ues of e. Biggest ben­efit of the more ad­vanced model is gained when e → 0.1.

Er­ror rates (in de­grees) and the fit­ted trend curve for 0 < e < 0.1 are shown in Figure 5. Within this region the error seems to grow approximately linearly in respect to e, and for Earth's e → 0.0167 the maximum error is less than 0.1°. It should be sufficient for most applications.

Figure 5: A trend of the max­imum er­ror (in de­grees) of mod­els for dif­fer­ent small val­ues of e.

Er­ror rates for e-val­ues of up-to 10-0.5 → 0.3 are show in Figure 6. In this context the error grows in proportion to e3.09, reaching a maximum error of 6° at e = 0.3. It might be too high for most applications, but this elliptic orbits are not common in at least our solar system. Our top-3 planets with highest eccentricities are Mercury (e = 0.21), Mars (e = 0.09) and Saturn (e = 0.05). Pluto has a remarkably high eccentricity of 0.25.

Figure 6: A trend of the max­imum er­ror (in de­grees) of mod­els for dif­fer­ent large val­ues of e.

Af­ter this model is in place, it en­ables us to cal­cu­late Earth's lo­ca­tion in re­spect to Sun in po­lar co­or­di­nates for a given times­tamp. It just needs to be con­verted to "sec­onds since Jan­uary 1 of this year", which is fur­ther con­verted to the mean anomaly M (in ra­di­ans). True anomaly θ is cal­cu­lated as

    θ e(M) = M - e · a(e) · (M - pi) · (M · (2π - M))-b(e)
    a(e) = a0 + exp(a1· log(e) + a2)
    b(e) = b0 + exp(b1· log(e) + b2) (7)

Es­ti­mated pa­ram­eters are:

    a0 = 0.0804
    a1 = 0.9869
    a2 = 1.0123
    b0 = -1.3848
    b1 = 0.6143
    b2 = -1.0715 (8)

The next step would be to es­ti­mate the lo­cal co­or­di­nate sys­tem from lat­itude, lon­gi­tude and time, and to use po­lar co­or­di­nates to cal­cu­late the per­son's ori­en­ta­tion with re­spect to the he­lio­cen­tric co­or­di­nate sys­tem where pe­ri­ap­sis (around 3.1. each year) is the "zero" di­rec­tion. Then the Sun's ap­par­ent al­ti­tude and az­imuth can be eas­ily cal­cu­lated.

Related blog posts: