Apsides Calculation

Question

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).

Answer 1

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 b_i values suggests that A₁ is for the Moon, A₂ and A₃ are for Venus, and A₄ and A₅ 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.

Answer 2

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.

MJD_TDB = 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.

Answer 3

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.

Answer 4

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.