Time of the seasons

[cdeaglejr]

This HP Prime program can be used to determine the UTC calendar date and time of the equinoxes and solstices of the Earth. These events are the times when the apparent geocentric longitude of the Sun is an exact multiple of 90 degrees. This script uses Brent’s root-finder and a precision solar ephemeris to calculate these events.

Brent’s method requires an objective function that defines the nonlinear equation to be solved. The objective function for the spring and fall equinoxes is the geocentric declination of the Sun. The spring and fall equinoxes occur whenever the geocentric declination of the Sun is less than or equal to a user-defined convergence criterion.

For the summer and winter solstices, the objective function is theta - lambda, where theta = 90 degrees for the summer solstice, theta = 270 degrees for the winter solstice, and lambda is the geocentric longitude of the Sun. The summer and winter solstices occur whenever the difference delta = theta - lambda is less than or equal to a user-defined convergence criterion.

Brent’s method also requires an initial and final time which bounds the root of the objective function. The initial time for the spring equinox is March 15, for the summer solstice June 15, for the fall equinox September 15 and for the winter solstice December 15. For each event, the final time is equal to these initial dates plus 10 days.

[Gil]

Question

How accurate, say for the present year 2023 or for 2022, are beginnings of the seasons in TT when compared with the official values in TT (for example as given by IMCEE)?

Thanks for your answer.

Gil

[cdeaglejr]

Naturally, the answer to your question depends on the algorithms used to compute the time of the seasons. Since this HP program computes these times geometrically, the prediction depends on the ephemeris used to calculate the position of the sun/earth. It also depends on the fidelity of the root-finder algorithm used to predict the events. Finally, it also depends on the conversion of TDB time to UTC (leap seconds, UTC-UTC1).

Most sources for online information about seasons (IMCEE, US Naval Observatory) provide the data to the nearest minute. This algorithm agrees with these sources to within that minute.

Questions always lead to more questions.

For example, how would you actually measure the time of the seasons?

[Gil]

My reason of my request was because
I try to reproduce the Meeus' algorithm to get the June solstice for the year 1962,with correct answer in Terrestrial Time or (Terrestrial) Dynamical Time within 2-4 seconds.

I work out author's example 27.b (pages 180-181),
using the tables given in his book at pages 418-421.

I get the following results on my HP50G:


:JDE0: 2437837.38588 <Jean Meeus... 35889>
tau = (JDEO-2451545)/365230
:SL0: 176069483.423
:SL1: 628332128985.
:SL2: 47714.5774468
:SL3: -192.115254877
:SL4: -109.4139133
:SL5: -.999998731728
'(SL0+SL1*tau+SL2*tau^2+SL3*tau^3+SL4*tau^4+SL5*tau^5)/100000000'
:L[rad]: -234.048595755 :
L[d.d]: 270.0032623
:R: 1.01630134777

with these L and R values diverging slightly
from the ones given in the book on the top of page 181.

[Gil]

My reason of my request was because
I try to reproduce the Meeus' algorithm to get the June solstice for the year 1962,with correct answer in Terrestrial Time or (Terrestrial) Dynamical Time within 2-4 seconds.

I work out author's example 27.b (pages 180-181),
using the tables given in his book at pages 418-421.

I get the following results on my HP50G:


:JDE0: 2437837.38588 <Jean Meeus... 35889>
tau = (JDEO-2451545)/365230
:SL0: 176069483.423
:SL1: 628332128985.
:SL2: 47714.5774468
:SL3: -192.115254877
:SL4: -109.4139133
:SL5: -.999998731728
'(SL0+SL1*tau+SL2*tau^2+SL3*tau^3+SL4*tau^4+SL5*tau^5)/100000000'
:L[rad]: -234.048595755 :
L[d.d]: 270.0032623
:R: 1.01630134777

with these L and R values diverging slightly
from the ones given in the book on the top of page 181.

[Giuseppe Donnini]

The discrepancy in the last digit between Meeus’ result and yours stems from the fact that you calculated the polynomial directly:

$$
JDE_{0}=\,2\,451\,716.56767 + 365\,241.62603\cdot Y + 0.00325\cdot Y^2 + 0.00888\cdot Y^3 - 0.00030\cdot Y^4
$$

instead of using Horner’s method, as Meeus himself recommends (see the general discussion on pp. 10-11):

$$
\Leftrightarrow \;2\,451\,716.56767 + Y\cdot ( \,365\,241.62603 + Y\cdot ( \,0.00325 + Y\cdot ( \,0.00888 - Y\cdot 0.00030\,) \,) \,)
$$

In general, the latter is both faster and – most importantly – more accurate:

1. If n is the degree of the polynomial, direct evaluation would require n(n+1)/2 multiplications, while Horner’s method only requires n. The larger n, the larger the speed gain.
    
2. Fewer operations automatically translate into less accumulation of roundoff errors.
   

[Gil]

Thanks for your explanation.

I will try it soon, but now I lost some parts of the main the program...

Regards, Gil