# Apsides Calculation

I’m trying to calculate the precise moment of Earth’s apsides (perihelion and aphelion). The only formulas I find are in Astronomical Algorithms by Jean Meeus (1998). However, they give erroneous times for 2021 (and, presumably, other dates “far” from 1998). I looked on ADSABS and on arXiv for possible other formulas, but found nothing.

Can someone point me to more modern formulas for calculating the moments of Earth’s apsides?

NOTE: I don’t want software modules; I want the raw math (unless the raw math is easily found from the software module).

• What values do you consider authoritative? How large an error do you get? What error would be acceptable? – Mike G Mar 6 at 21:16
• I would consider authoritative the values found on Wikipedia or on neoprogrammics.com/PHPSL_DIRECTORY/… or JPL. I get an error of more than a day. I would like an error of less than an hour. – Pierre Paquette Mar 6 at 21:28
• I would not consider values found on Wikipedia authoritative. – David Hammen Mar 7 at 1:46
• I also would not consider neoprogrammics.com authoritative. For one thing, it is based on DE405, which is out of date, and for another, it is a random website by one guy. – David Hammen Mar 7 at 1:53
• @B--rian If WolframAlpha is doing things right, it is using table-driven data provided by some ephemeris for dates in the near past / near future, approximate formulae for dates outside that range, and should report "your question is bogus" for dates outside the validity range of those approximate formulae. It does not do the latter; I can, for example ask it for Earth perihelion ten million years from now. It gives an answer, but the answer is almost certainly is bogus. – David Hammen Mar 7 at 13:20

For those who don't have ready access to a copy of Astronomical Algorithms, Meeus's first approximation looks like:

$$\text{JDE} = 2541547.51 + 365.259636 ~k + 1.6 \times 10^{-8} ~k^2$$

where k, the number of anomalistic years since the January 2000 perihelion, is half of some integer. This neglects the influence of the Moon and other planets, so he adds a correction:

$$\sum_{i=1}^5 a_i \sin A_i$$

where

$$A_i = c_i + b_i k .$$

Inspection of Meeus's bi values suggests that A1 is for the Moon, A2 and A3 are for Venus, and A4 and A5 are for Jupiter. He reports a mean error of 3 hours for years 1980-2019, presumably compared to VSOP87. Comparing his approximation to JPL DE430 data, I compute an RMS error of 3.3 hours not only for that interval but also for years 1800-2199.

If I use a curve fitter to tweak the correction parameters for years 2010-2049, I can shrink the mean error in that interval from 3.8 to 1.6 hours, but then the error outside that interval is worse. An approximation with more terms could achieve better accuracy, but that puts you on the road toward reinventing VSOP87.

Given the requirements of sub-hour accuracy and no multi-megabyte ephemeris tables, VSOP87 or a subset of it seems to be the way to go.

• I'm interested in the curve fitter you used to find a, b and c; it's publicly available? – Cristiano Mar 19 at 8:47
• @Cristiano I used scipy.optimize.curve_fit. – Mike G Mar 19 at 10:21
• I see the option method='lm' (which is Levenberg-Marquardt I used). Since with my implementation it seems very hard to find good Ai's, it would be interesting to see your coefficients and (if it's not too much trouble) also the code for the calculations. – Cristiano Mar 19 at 10:23
• @Cristiano meeus38.py – Mike G Mar 19 at 21:55

Do you need 1000 years? 50 years? 1-hour accuracy? 1-second accuracy?

A simple linear interpolation for the perihelia from the year 2000 to 2050 gives a maximum error of about 1.3 days for the year 2009 and a mean absolute error of about 19.3 hours:

et = 31558511.31638778 * year - 63116806104.00429

et is the so called ephemeris time (used by NAIF team in the SPICE library), it’s the number of TDB seconds past J2000.

MJDTDB = et / 86400 + 51544.5

The calculations are done with the SPICE library and DE440.

If you need an improved formula, we have:

t = 1035263.906713132 + year * 810.005204205729

et = -103451133434.8591 + t * 38960.7024090397 + 106967.8662659288 * sin(t * 0.2065792456689355 - 5.15924414294114)

sin in radians. The maximum error is 7.2 hours for the year 2022 and the mean absolute error is 2.4 hours.

EDIT: proof that my formulas work

Let’s do the calculations for the year 2021.
Calculate 31558511.31638778 * 2021 - 63116806104.00429
the result is 662945266.41541338 TDB seconds past J2000.

Now use a fully authoritative site to convert that time to UTC time.
Go to WebGeocalc https://wgc.jpl.nasa.gov:8443/webgeocalc/#NewCalculation scroll all the way to the bottom to see “Time Conversion” and click the link.

Kernel selection: leave it empty
Input Time
Time system: TDB
Time format: Seconds past J2000
Input times: Single time
Time: copy and paste the result (662945266.41541338)
Output Time
Time system: UTC
Time format: Calendar (year-month-day)

Then press “Calculate”, you get: 2021-01-03 11:26:37.231413 UTC, which is what I wrote.

• @PierrePaquette I wrote: "The calculations are done with the SPICE library and DE440." With the improved formula I get 2021-01-02 09:44:06.212 UTC (the formula gives TDB, here I'm converting from TDB to UTC), while the correct time is 2021-01-02 13:50:34.634 UTC, the error is only 4.1 h. With the linear interpolation I get 2021-01-03 11:26:37.231 UTC, the error is 21.6 h which is well inside the 1.3 days I wrote. It seems a bit unfair to me to downvote just because you don't now ho to do the calculations. – Cristiano Mar 8 at 9:07
• I’m sorry about the downvoting; you are right that it is unfair. I’ll fix that as soon as I can. However, I did download the SPICE library from JPL’s website, and I can’t find those values anywhere… But maybe I just don’t know where to look. Also, is there a paper published with the formulas? – Pierre Paquette Mar 9 at 0:02
• …and I am dumbfounded. I just redid the calculations with your formulas, and I get the same values. But when I tried 23 hours ago, it wasn’t working! Maybe I had a typo (but I copy/pasted?) or I don’t know what. (Side note… I was calculating the positions of Jupiter’s satellites today, and it wasn’t working, then I tried some debugging, found nothing, and all of a sudden, it was working. Ugh! I think I’m going nuts.) – Pierre Paquette Mar 9 at 0:28
• I’m sincerely sorry: It seems I can’t unvote my downvote unless your answer is edited. :-( – Pierre Paquette Mar 9 at 3:34

I post another answer because I think it should work better for you.

Since I read in one your comment that you “only need a few years past and (mostly) future from 2021” and that you need “at least one-hour accuracy”, here’s my proposed method.

I still don't know how many years you really need, but the following table shows the perihelia and aphelia with a 1-second formal accuracy and I claim that the result is authoritative because I did the calculations with the SPICE library (published by the authoritative NAIF team) using the DE440 ephemeris (published by the authoritative JPL):

1990-01-04 17:22:34   1990-07-04 05:04:03
1991-01-03 02:59:07   1991-07-06 15:26:27
1992-01-03 15:02:25   1992-07-03 12:06:58
1993-01-04 03:03:46   1993-07-04 22:20:50
1994-01-02 05:54:11   1994-07-05 19:16:43
1995-01-04 11:05:21   1995-07-04 02:16:43
1996-01-04 07:24:51   1996-07-05 18:59:52
1997-01-01 23:16:02   1997-07-04 19:19:20
1998-01-04 21:15:00   1998-07-03 23:50:13
1999-01-03 13:00:09   1999-07-06 22:50:46
2000-01-03 05:17:41   2000-07-03 23:48:55
2001-01-04 08:52:15   2001-07-04 13:37:08
2002-01-02 14:08:45   2002-07-06 03:46:47
2003-01-04 05:01:44   2003-07-04 05:39:38
2004-01-04 17:41:57   2004-07-05 10:53:43
2005-01-02 00:35:17   2005-07-05 04:57:52
2006-01-04 15:29:38   2006-07-03 23:09:59
2007-01-03 19:42:57   2007-07-06 23:52:35
2008-01-02 23:51:08   2008-07-04 07:40:54
2009-01-04 15:29:40   2009-07-04 01:40:19
2010-01-03 00:09:16   2010-07-06 11:29:58
2011-01-03 18:32:00   2011-07-04 14:53:59
2012-01-05 00:31:51   2012-07-05 03:32:16
2013-01-02 04:37:35   2013-07-05 14:44:24
2014-01-04 11:58:36   2014-07-04 00:13:28
2015-01-04 06:36:11   2015-07-06 19:40:23
2016-01-02 22:48:48   2016-07-04 16:24:14
2017-01-04 14:17:50   2017-07-03 20:11:22
2018-01-03 05:34:44   2018-07-06 16:46:47
2019-01-03 05:20:00   2019-07-04 22:10:49
2020-01-05 07:47:56   2020-07-04 11:34:44
2021-01-02 13:50:35   2021-07-05 22:27:26
2022-01-04 06:54:39   2022-07-04 07:10:44
2023-01-04 16:17:28   2023-07-06 20:06:39
2024-01-03 00:38:37   2024-07-05 05:06:04
2025-01-04 13:28:07   2025-07-03 19:54:43
2026-01-03 17:15:39   2026-07-06 17:30:39
2027-01-03 02:32:46   2027-07-05 05:05:50
2028-01-05 12:28:23   2028-07-03 22:18:06
2029-01-02 18:13:34   2029-07-06 05:11:55
2030-01-03 10:12:35   2030-07-04 12:57:43
2031-01-04 20:47:53   2031-07-06 07:10:07
2032-01-03 05:11:22   2032-07-05 11:53:37
2033-01-04 11:51:21   2033-07-03 20:51:59
2034-01-04 04:46:59   2034-07-06 18:49:15
2035-01-03 00:54:15   2035-07-05 18:21:43
2036-01-05 14:17:09   2036-07-03 21:17:32
2037-01-03 04:00:33   2037-07-06 12:05:28
2038-01-03 05:01:32   2038-07-04 19:46:07
2039-01-05 06:41:38   2039-07-05 13:25:17
2040-01-03 11:32:51   2040-07-05 19:01:46
2041-01-03 21:52:01   2041-07-04 01:38:34


The only non-authoritative thing is the code I wrote to do the calculations, but it’s so simple and I tested it against known results that I can claim that the table is authoritative without any doubt.
But you may think that I am the biggest troll in the world who is here only to waste his time with you, so you are totally legitimate to double check that table.
Then we use again WebGeocalc to verify my claim.

Click “Position Finder”
Kernel selection: you should use DE440 to obtain exactly my values, but you could probably leave it empty
Target: 399 (we are absolutely sure that it’s the Earth and not the Earth-Moon barycenter)
Observer: SUN
Reference frame: J2000
Light propagation: none (geometric state)
Input time: choose the format you like
Start time and Stop time: you may want to verify the distance a few hours before and after my tabulated times; for the perihelion of the year 2041 you could put 2041-01-03, 2041-01-04 12:00 with a 10-minute time step.

Coordinate system: Rectangular
Coordinate condition: Distance is local minimum

Press “Calculate” and you should see: 2041-01-03 21:52:01... UTC

You may want to check a few rows just to realize that I’m not a stupid troll and that my table is authoritative for the reasons I already wrote.

EDIT

Here's the C++ code to do the calculations:

furnsh_c(NAIF_DIR"naif0012.tls.pc"); // Leap seconds
furnsh_c(NAIF_DIR"pck00010.tpc"); // Reference frames
furnsh_c(NAIF_DIR"de440.bsp"); // Ephemeris

const int TGT= 399; // Target body
const int start= 1990, end= 2040; // Years

for(int i= start-1; i < end; i++) {
double a; char buf[32]; sprintf(buf, "%d-12-15", i); str2et_c(buf, &a);
for(int j= 0; j < 2; j++) {
double b= a + 30 * 86400, pos[3], lt;
while(1) {
double c = b - (b - a) / 1.9; spkgps_c(TGT, c, "J2000", 10, pos, &lt); double fc= vnorm_c(pos);
double d = a + (b - a) / 1.9; spkgps_c(TGT, d, "J2000", 10, pos, &lt); double fd= vnorm_c(pos);
if(j == 0) { if(fc < fd) b = d; else a = c; } // Perihelion
else       { if(fc > fd) b = d; else a = c; } // Aphelion
if(b - a < .1) break; // Uncertainty [s]
}
timout_c((a+b)/2, "YYYY-MM-DD HR:MN:SC ::UTC ::RND", sizeof(buf), buf);
printf("%s", buf); if(j == 0) printf("   ");
a += 180 * 86400;
}
puts("");
}


it's a simple Golden-section search: https://en.wikipedia.org/wiki/Golden-section_search to find the minimum and the maximum Earth-Sun distance.

After you have created the table, you need an interpolation algorithm to calculate the formula. I use the Levenberg Marquardt Least Squares Fitting algorithm: https://github.com/mattjr/structured/blob/master/CMVS-PMVS/program/thirdParty/lmfit-3.2/doc/lmfit.pod, but it's not that easy to use, it will take some time to figure out how to use it.

• Thank you for this as well. However, I want to use calculation and not “database polling,” so I would appreciate learning how to use the SPICE library—if only for this!—as opposed to consulting from a list, which I could get anywhere. So… how do I use the SPICE library? (And I never, ever thought you could be a troll…) – Pierre Paquette Mar 9 at 23:57
• You finally said that you need an analytical solution for the apsides. As long as you use the accurate DExxx ephemerides, it' impossible to find an analytical solution (it doesn't exist). An analytical solution might exist if you expressed the position by a relatively simple trigonometric expansion, but in that case you would get largely inaccurate results (if a low order expansion is used). – Cristiano Mar 11 at 9:21
• I ended up using a VSOP87 implementation I have, calculating the Earth–Sun distance on a few days around perihelion and around aphelion, interpolating a function for these values, and finding the minimum (perihelion) or maximum (aphelion) of this function. I get ~1 minute precision compared to the times you provided. Thank you again for your help. Question: I’m curious… What’s your background? – Pierre Paquette Mar 11 at 23:17
• @PierrePaquette Although terribly outdated, I can confirm a surprisingly small error (usually less than 1 minute), even if the calculations are much, much slower. My background? Nothing exceptional, I wrote my first astronomy program in 1981 (I bought Meeus' book) using a Commodore VIC 20 and I never stopped, but now I use a different computer! :-) – Cristiano Mar 12 at 8:55
• Well, I was having a hard time figuring out how to use VSOP2013, but I finally got it! – Pierre Paquette Mar 13 at 0:34

There is no easy formula for what you want. Astronomical Algorithms provides an algorithm for the times of perihelion and aphelion passages of the Earth-Moon barycenter. The times at which the Earth itself is closest to / furthest from the Sun is made much more complex by the fact that the Earth and Moon orbit one another as well as the Sun. Because perihelion and aphelion passage are when the distances between the Earth-Moon barycenter and the Sun change slowest, the phase of the Moon has a significant multi-day impact on the timing of the times of perihelion and aphelion passages of the Earth itself.

• Astronomical Algorithms provides a correction for calculating the moment of Earth’s passage at perihelion or aphelion. But it’s still off from values given by multiple websites. – Pierre Paquette Mar 7 at 6:33
• @PierrePaquette Off by how much? Meeus claims in Astronomical Algorithms that even with his correction, that the error can be as much as six hours from the years 1980-2019. For 2021, I get an error of about five hours using Meeus's formulae, which is not that bad given what Meeus himself wrote about the accuracy of those calculations. – David Hammen Mar 7 at 12:48
• Yes, five hours off for the perihelion, but 47 hours off for the aphelion… That is way too much for a website I’m building that will claim to be giving quality information… If my viewers compare with other sources (and I know they will!), they’ll say I’m not credible and won’t come back… – Pierre Paquette Mar 8 at 1:14