<p>For those who don't have ready access to a copy of <em>Astronomical Algorithms</em>,
Meeus's first approximation looks like:</p>
<p><span class="math-container">$$ \text{JDE} = 2541547.51 + 365.259636 ~k + 1.6 \times 10^{-8} ~k^2 $$</span></p>
<p>where <em>k</em>, 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:</p>
<p><span class="math-container">$$ \sum_{i=1}^5 a_i \sin A_i $$</span></p>
<p>where</p>
<p><span class="math-container">$$ A_i = c_i + b_i k .$$</span></p>
<p>Inspection of Meeus's <em>b<sub>i</sub></em> values suggests that
<em>A</em><sub>1</sub> is for the Moon,
<em>A</em><sub>2</sub> and <em>A</em><sub>3</sub> are for Venus,
and <em>A</em><sub>4</sub> and <em>A</em><sub>5</sub> 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.</p>
<p>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.</p>
<p>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.</p>

CC BY-SA 4.0
2021-03-19T08:47:01.143Z
I'm interested in the curve fitter you used to find a, b and c; it's publicly available?
Cristiano
CC BY-SA 4.0
2021-03-19T10:21:23.840Z
@Cristiano I used [scipy.optimize](https://docs.scipy.org/doc/scipy/reference/optimize.html).curve_fit.
Mike G
CC BY-SA 4.0
2021-03-19T10:23:33.003Z
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
CC BY-SA 4.0
2021-03-19T21:55:35.807Z
@Cristiano [meeus38.py](https://github.com/mgarraha/astro-se/blob/main/meeus38.py)
Mike G
CC BY-SA 4.0
2021-03-19T01:24:48.883Z
2021-03-19T01:24:48.883Z
Mike G
<p>Do you need 1000 years? 50 years? 1-hour accuracy? 1-second accuracy?</p>
<p>A simple linear interpolation <strong>for the perihelia</strong> 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:</p>
<p>et = 31558511.31638778 * year - 63116806104.00429</p>
<p>et is the so called ephemeris time (used by NAIF team in the SPICE library), it’s the number of TDB seconds past J2000.</p>
<p>MJD<sub>TDB</sub> = et / 86400 + 51544.5</p>
<p>The calculations are done with the SPICE library and DE440.</p>
<p>If you need an improved formula, we have:</p>
<p>t = 1035263.906713132 + year * 810.005204205729</p>
<p>et = -103451133434.8591 + t * 38960.7024090397 + 106967.8662659288 * sin(t * 0.2065792456689355 - 5.15924414294114)</p>
<p>sin in radians. The maximum error is 7.2 hours for the year 2022 and the mean absolute error is 2.4 hours.</p>
<p><strong>EDIT: proof that my formulas work</strong></p>
<p>Let’s do the calculations for the year 2021.<br>
Calculate 31558511.31638778 * 2021 - 63116806104.00429<br>
the result is 662945266.41541338 TDB seconds past J2000.<br>
<br>
Now use a fully authoritative site to convert that time to UTC time.<br>
Go to WebGeocalc <a href="https://wgc.jpl.nasa.gov:8443/webgeocalc/#NewCalculation" rel="noreferrer">https://wgc.jpl.nasa.gov:8443/webgeocalc/#NewCalculation</a>
scroll all the way to the bottom to see “Time Conversion” and click the link.<br>
<br>
Kernel selection: leave it empty<br>
Input Time<br>
Time system: TDB<br>
Time format: Seconds past J2000<br>
Input times: Single time<br>
Time: copy and paste the result (662945266.41541338)<br>
Output Time<br>
Time system: UTC<br>
Time format: Calendar (year-month-day)<br>
<br>
Then press “Calculate”, you get: 2021-01-03 11:26:37.231413 UTC, which is what I wrote.</p>

CC BY-SA 4.0
2021-03-08T09:07:05.320Z
@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
CC BY-SA 4.0
2021-03-09T00:02:38.973Z
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
CC BY-SA 4.0
2021-03-09T00:28:18.483Z
…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
CC BY-SA 4.0
2021-03-09T03:34:17.210Z
I’m sincerely sorry: It seems I can’t unvote my downvote unless your answer is edited. :-(
Pierre Paquette
CC BY-SA 4.0
2021-03-07T22:16:48.017Z
2021-03-09T16:23:55.280Z
2021-03-09T16:23:55.280Z
Cristiano
<p>I post another answer because I think it should work better for you.<br></p>
<p>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.</p>
<p>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):</p>
<pre><code>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
</code></pre>
<p>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.<br>
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.<br>
Then we use again WebGeocalc to verify my claim.<br></p>
<p>Click “Position Finder”<br>
Kernel selection: you should use DE440 to obtain exactly my values, but you could probably leave it empty<br>
Target: 399 (we are absolutely sure that it’s the Earth and not the Earth-Moon barycenter)<br>
Observer: SUN<br>
Reference frame: J2000<br>
Light propagation: none (geometric state)<br>
Input time: choose the format you like<br>
<strong>Start time and Stop time</strong>: 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.<br>
<br>
Coordinate system: Rectangular<br>
Coordinate condition: Distance is local minimum<br>
<br>
Press “Calculate” and you should see: 2041-01-03 21:52:01... UTC<br>
<br>
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.<br></p>
<p><strong>EDIT</strong></p>
<p>Here's the C++ code to do the calculations:</p>
<pre><code>furnsh_c(NAIF_DIR&quot;naif0012.tls.pc&quot;); // Leap seconds
furnsh_c(NAIF_DIR&quot;pck00010.tpc&quot;); // Reference frames
furnsh_c(NAIF_DIR&quot;de440.bsp&quot;); // Ephemeris

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

for(int i= start-1; i &lt; end; i++) {
    double a; char buf[32]; sprintf(buf, &quot;%d-12-15&quot;, i); str2et_c(buf, &amp;a);
    for(int j= 0; j &lt; 2; j++) {
        double b= a + 30 * 86400, pos[3], lt;
        while(1) {
            double c = b - (b - a) / 1.9; spkgps_c(TGT, c, &quot;J2000&quot;, 10, pos, &amp;lt); double fc= vnorm_c(pos);
            double d = a + (b - a) / 1.9; spkgps_c(TGT, d, &quot;J2000&quot;, 10, pos, &amp;lt); double fd= vnorm_c(pos);
            if(j == 0) { if(fc &lt; fd) b = d; else a = c; } // Perihelion
            else       { if(fc &gt; fd) b = d; else a = c; } // Aphelion
            if(b - a &lt; .1) break; // Uncertainty [s]
        }
        timout_c((a+b)/2, &quot;YYYY-MM-DD HR:MN:SC ::UTC ::RND&quot;, sizeof(buf), buf);
        printf(&quot;%s&quot;, buf); if(j == 0) printf(&quot;   &quot;);
        a += 180 * 86400;
    }
    puts(&quot;&quot;);
}
</code></pre>
<p>it's a simple Golden-section search: <a href="https://en.wikipedia.org/wiki/Golden-section_search" rel="nofollow noreferrer">https://en.wikipedia.org/wiki/Golden-section_search</a> to find the minimum and the maximum Earth-Sun distance.<br></p>
<p>After you have created the table, you need an interpolation algorithm to calculate the formula. I use the Levenberg Marquardt Least Squares Fitting algorithm: <a href="https://github.com/mattjr/structured/blob/master/CMVS-PMVS/program/thirdParty/lmfit-3.2/doc/lmfit.pod" rel="nofollow noreferrer">https://github.com/mattjr/structured/blob/master/CMVS-PMVS/program/thirdParty/lmfit-3.2/doc/lmfit.pod</a>, but it's not that easy to use, it will take some time to figure out how to use it.</p>

CC BY-SA 4.0
2021-03-09T23:57:59.503Z
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
CC BY-SA 4.0
2021-03-10T00:47:44.510Z
In other words, may I please have the source of your code instead of just the result? My question was for the formulas (no matter how complex they might be), not the dates.
Pierre Paquette
CC BY-SA 4.0
2021-03-10T23:07:23.250Z
Thank you for your edits. Much appreciated. I see there’s no direct formula for the calculation of the apsides. One would think that Meeus’s formulas could just have been improved on…
Pierre Paquette
CC BY-SA 4.0
2021-03-11T09:21:35.367Z
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
CC BY-SA 4.0
2021-03-11T23:17:06.067Z
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
CC BY-SA 4.0
2021-03-12T08:55:50.537Z
@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
CC BY-SA 4.0
2021-03-13T00:34:46.583Z
Well, I was having a hard time figuring out how to use VSOP2013, but I finally got it!
Pierre Paquette
CC BY-SA 4.0
2021-03-13T09:28:30.257Z
@PierrePaquette I implemented VSOP2013 theory both in C++ and VBA. VSOP2013 doesn't include the Earth, it only include EMB; in your case, it should give worse results.
Cristiano
CC BY-SA 4.0
2021-03-13T19:45:15.933Z
It does, indeed. However, for planetary positions, it’s better. Too bad it’s too big to use as JavaScript in a webpage. Even in PHP, the calculation is slow enough that visitors would click away.
Pierre Paquette
CC BY-SA 4.0
2021-03-09T16:13:27.473Z
2021-03-10T09:35:24.947Z
2021-03-10T09:35:24.947Z
Cristiano
<p>There is no easy formula for what you want. <em>Astronomical Algorithms</em> 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.</p>

CC BY-SA 4.0
2021-03-07T06:33:00.830Z
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
CC BY-SA 4.0
2021-03-07T12:48:41.380Z
@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
CC BY-SA 4.0
2021-03-08T01:14:41.957Z
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
CC BY-SA 4.0
2021-03-07T02:11:32.297Z
2021-03-07T02:11:32.297Z
David Hammen
<p>I’m trying to calculate the precise moment of Earth’s apsides (perihelion and aphelion). The only formulas I find are in <em>Astronomical Algorithms</em> 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.</p>
<p>Can someone point me to more modern formulas for calculating the moments of Earth’s apsides?</p>
<p>NOTE: I don’t want software modules; I want the raw math (unless the raw math is easily found from the software module).</p>

CC BY-SA 4.0
2021-03-06T21:16:50.450Z
What values do you consider authoritative? How large an error do you get? What error would be acceptable?
Mike G
CC BY-SA 4.0
2021-03-06T21:28:45.453Z
I would consider authoritative the values found on Wikipedia or on https://neoprogrammics.com/PHPSL_DIRECTORY/Earth_Perihelions_and_Intervals_Table/ or JPL. I get an error of more than a day. I would like an error of less than an hour.
Pierre Paquette
CC BY-SA 4.0
2021-03-07T01:46:48.903Z
I would not consider values found on Wikipedia authoritative.
David Hammen
CC BY-SA 4.0
2021-03-07T01:53:29.943Z
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
CC BY-SA 4.0
2021-03-07T06:33:26.467Z
DE405 is still more recent than Astronomical Algorithms…
Pierre Paquette
CC BY-SA 4.0
2021-03-07T07:00:20.173Z
Wolframalpha calculates only with a precision up to a day, see https://www.wolframalpha.com/input/?i=Earth+next+periapsis but I could not find the algorithm used, although I know it was possible in Mathematica some years ago.
B--rian
CC BY-SA 4.0
2021-03-07T13:20:41.113Z
@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](https://www.wolframalpha.com/input/?i=earth+perihelion+ten+million+years+from+now). It gives an answer, but the answer is almost certainly is bogus.
David Hammen
<celestial-mechanics><ephemeris>
Apsides Calculation
Pierre Paquette
# Apsides Calculation

(Asked by: username_0 on 2021-03-06)

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

**Tags:** celestial-mechanics, ephemeris
> **Comments:**
> **username_1**: What values do you consider authoritative? How large an error do you get? What error would be acceptable?  
> **username_0**: I would consider authoritative the values found on Wikipedia or on https://neoprogrammics.com/PHPSL_DIRECTORY/Earth_Perihelions_and_Intervals_Table/ or JPL. I get an error of more than a day. I would like an error of less than an hour.  
> **username_2**: I would not consider values found on Wikipedia authoritative.  
> **username_2**: 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.  
> **username_0**: DE405 is still more recent than Astronomical Algorithms…  
> **username_3**: Wolframalpha calculates only with a precision up to a day, see https://www.wolframalpha.com/input/?i=Earth+next+periapsis but I could not find the algorithm used, although I know it was possible in Mathematica some years ago.  
> **username_2**: @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](https://www.wolframalpha.com/input/?i=earth+perihelion+ten+million+years+from+now). It gives an answer, but the answer is almost certainly is bogus.  

---

## 4 answers

### Answer by username_1 on 2021-03-19. Score: 6

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
*A*1 is for the Moon,
*A*2 and *A*3 are for Venus,
and *A*4 and *A*5 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.

> **Comments:**
> **username_4**: I'm interested in the curve fitter you used to find a, b and c; it's publicly available?  
> **username_1**: @Cristiano I used [scipy.optimize](https://docs.scipy.org/doc/scipy/reference/optimize.html).curve_fit.  
> **username_4**: 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.  
> **username_1**: @Cristiano [meeus38.py](https://github.com/mgarraha/astro-se/blob/main/meeus38.py)  

---

### Answer by username_4 on 2021-03-07. Score: 5

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.

> **Comments:**
> **username_4**: @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.  
> **username_0**: 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?  
> **username_0**: …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.)  
> **username_0**: I’m sincerely sorry: It seems I can’t unvote my downvote unless your answer is edited. :-(  

---

### Answer by username_4 on 2021-03-09. Score: 4

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, <); double fc= vnorm_c(pos);
            double d = a + (b - a) / 1.9; spkgps_c(TGT, d, "J2000", 10, pos, <); 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.

> **Comments:**
> **username_0**: 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…)  
> **username_0**: In other words, may I please have the source of your code instead of just the result? My question was for the formulas (no matter how complex they might be), not the dates.  
> **username_0**: Thank you for your edits. Much appreciated. I see there’s no direct formula for the calculation of the apsides. One would think that Meeus’s formulas could just have been improved on…  
> **username_4**: 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).  
> **username_0**: 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?  
> **username_4**: @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! :-)  
> **username_0**: Well, I was having a hard time figuring out how to use VSOP2013, but I finally got it!  
> **username_4**: @PierrePaquette I implemented VSOP2013 theory both in C++ and VBA. VSOP2013 doesn't include the Earth, it only include EMB; in your case, it should give worse results.  
> **username_0**: It does, indeed. However, for planetary positions, it’s better. Too bad it’s too big to use as JavaScript in a webpage. Even in PHP, the calculation is slow enough that visitors would click away.  

---

### Answer by username_2 on 2021-03-07. Score: 3

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.

> **Comments:**
> **username_0**: 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.  
> **username_2**: @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.  
> **username_0**: 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…  

---

