Astro-Sunrise
view release on metacpan or search on metacpan
doc/astronomical-notes.pod view on Meta::CPAN
-*- encoding: utf-8; indent-tabs-mode:nil -*-
=encoding utf-8
=head1 Document Status
This text is published under the I<Creative Commons> license
CC-BY-ND.
Copyright (c) 2017-2021 Jean Forget. All rights reserved.
I must precise that I am not a professional astronomer. The text
below may contain errors, be aware of this. I will not be held
responsible for any consequences of your reading of this text.
The "NO WARRANTY" paragraphs from the GPL and the Artistic License
apply not only to Perl code, but also to English (and French)
texts.
The text is often (but irregularly) updated on Github. There are
a French version and an English version. Since I am more at ease
discussing astronomical subjects in French, the English version
will lag behind the French one.
This text is an integral part of the module's distribution package.
So you can read it on web pages generated from CPAN
(for example L<https://metacpan.org>).
But it is not used during the module installation process.
So, I guess it will not appear in C<.deb> or C<.rpm> packages.
Although this text is stored in the L<Astro::Sunrise> repository, it
also documents the L<DateTime::Event::Sunrise> module, which has a
very similar core (astronomical computations) and a different API.
Since both modules move at different speeds, it may happen that the
text you are reading is not synchronised with the L<Astro::Sunrise>
module.
=head1 Why This Text? For Whom?
The main purpose of this text is to explain how the sunrises
and sunsets are computed. These explanations are much too long
to be included into the module's POD section.
=head2 For Whom? For My Teddy Bear
Have you read
L<brian's Guide to Solving Any Perl Problem|http://www252.pair.com/~comdog/brian's_guide.pod>?
While most advices deal with debugging Perl code, a few advices have
a broader scope and apply to any intellectual problem.
One of these advices consists in "talking to the teddy-bear".
And do not pretend to talk while just forming sentences in your
mind. You must really talk in a clear voice in front of your
teddy-bear, with syntactically correct sentences.
Actually, in the present case, the topic is so big that I skipped to the next level and
I prefer writing (on GitHub) for my teddy-bear.
I write this text to tell my teddy-bear which problems I have
encountered while maintaining this module and how I fixed them.
But mainly, I write this to give him a detailed description of the
precise iterative algorithm, because
L<Paul Schlyter's explanations|https://www.stjarnhimlen.se/comp/riset.html#3>
are not detailed enough for my taste and there is no
compilable source available to check this algorithm (unlike the
L<simple version without iteration|https://www.stjarnhimlen.se/comp/sunriset.c>).
=head2 For Whom? For The Next Module Maintainer
Actually, my teddy-bear understands French, so the present
English version is not for him.
The second person for whom I write is the next module maintainer. I have read
L<Neil Bowers' message|http://codeverge.com/perl.module-authors/the-module-authors-pledge/744969>
about I<the module authors pledge>. I agree with him and I declare that
should I stop maintaining my modules for whatever reason, I accept that
any volunteer can take charge of them.
What Neil did not explain, is that the new maintainer must obey a few
criteria and must have three available resources to take over a module maintenance:
be competent in Perl programming, have enough available time to work on
the module and be enthusiastic enough to get around to it.
In the case of astronomical module, the competence in Perl programming is
not enough, you must also be competent in astronomy. So, if you think you might
maintain this module, first read the present text. If you understand why I bother
about such and such question, if you can follow my train of thought without being
lost, then you are competent enough. If you think I am playing
L<Captain Obvious|https://tvtropes.org/pmwiki/pmwiki.php/Main/CaptainObvious>
and if you have instant answers to my questions, then you are the ideal
person that could maintain this module. If you do not understand what all
this is about, and if sines and cosines put you off, do not consider
working on this module's innards.
=head2 For Whom? For Bug Reporters
This text is also for those who think they have found a bug
in the module or who want to offer an idea to improve the module.
doc/astronomical-notes.pod view on Meta::CPAN
I hate these people who, each time snow falls, cry
"Where is this global warming scientists talk about
again and again?" These people seem to ignore that
climate and weather are two different things. When
the temperature from a meteorologic station varies
by 10 degrees C from one day to the following, this
is a mundane meteorological event. When the I<average>
temperature for a decade varies by 2 degrees C from
one century to the next, this is a catastrophic
climate event.
The movements I explain below are more "climatic" and
less "meteorogical" than Earth's spin and orbital
rotation. Their values over a short timespan are so
low that the algorithms computing astronomical positions
over a short timespan do not care about them.
Note: weather (but not climate) will come back in a few
chapters as a real phenomenon, not as a metaphor.
=head3 Equinox Precession
The best known movement with a long timescale is the
equinox precession. Presently, the vernal point lies
within the constellation of Pisces. But actually, it moves
all along the ecliptic, making a whole turn in about
26,000 years.
=head3 Nutation
The angle between the equatorial plane and the ecliptic
plane varies slightly. In Paul Schlyter's C program, the
angle decreases by 356 nanodegrees per day (3.56e-7 °/d,
1.3e-4 °/yr).
=head3 Perihelion Precession
There is also the perihelion precession. This movement
is best known for Mercury, because it is the most apparent,
but all other planets have a perihelion precession, including
the Earth.
=head3 Other Drifts And Fluctuations
The formulas computing the positions of celestial bodies
use some constants. But these values are constant only
on a short timespan (astronomically speaking; or, with the
metaphor above, on a "meteorological" timespan). But they are
variable on a longer timespan (or a "climatic" timespan) For example,
everybody knows that the day lasts 24 hours (the mean
solar day, not the sidereal day). Yet, as I have read it
somewhere, in paleontological times, it used to last
22 hours or so.
The variation of the duration of the day is a tiny
variation, but with our modern measure instruments, we
can measure it. Since the time when scientists abolished the
astronomical standard of time for an atomic
standard, it has been necessary to add 27 leap seconds
over 47 years to synchronise the atomic timescale with
the Earth's spin.
For the moment, all adjustments have consisted in
adding a leap second. But it can happen that we
would have to synchronise in the other direction by removing
a second. So this phenomenon produces fluctuations
rather than a slow drift in a single direction.
=head3 The Equation Of Time
There are other fluctuations, easier to measure and with a more
"meteorological" and less "climatic" timescale. The I<true> solar noon
does not occur on the same precise time as the I<mean> solar noon.
There are two reasons.
=head4 Obliquity of the Earth
First, there is an angle between the ecliptical plane and the equatorial plane,
therefore, a constant-speed rotation on the ecliptical plane does not translate
to a constant-speed rotation when measured by right ascension on the equatorial
plane. The rate of variation of the right ascension is a variable rate.
If we use the same units for the ecliptic longitude and the right ascension
(either degrees or hours), then both values are nearly equal, but still different.
So, when the ecliptic longitude is 46°20'31", the right ascension is 43°52'36",
that is, a 2°27'54" gap. The same happens at longitude 226°20'31". And at
longitude 313°32'52", the right ascension is 316°47", that is a gap of 2°27'54",
but in the other direction. And the same happens at 133°32'52". These are
the maximum values for the gap when using an obliquity of 23°26'.
And if you prefer hours, here are the values:
. longitude right ascension gap longitude right ascension gap
. 3h05mn22s 2h55mn30s -9mn51s 46°20'31" 43°52'36" -2°27'54"
. 8h54mn11s 9h04mn03s 9mn51s 133°32'52" 136°00'47" 2°27'54"
. 15h05mn22s 14h55mn30s -9mn51s 226°20'31" 223°52'36" -2°27'54"
. 20h54mn11s 21h04mn03s 9mn51s 313°32'52" 316°00'47" 2°27'54"
=head4 Kepler's Second Law
Second, the rotational speed of Sun itself on the ecliptical plane is not a constant.
It obeys Kepler's second law, with a rotational speed more or less inversely
proportional to the Earth-Sun distance.
Q: You cannot apply Kepler's second law to a geocentric model!
A: No. Kepler's second law applies to a barycentric model as D above. It applies
I<approximately> to an heliocentric model as C. But once we have computed
Earth's angular speed on its orbit around the Sun in model C, the computation of the
Sun's coordinates and speed in the geocentric model is very simple. Especially,
the Sun's angular speed in a geocentric model is equal to the Earth's speed in
an heliocentric model.
Here are the Sun's positions for 2017, as given by Stellarium. The software
gives the equatorial coordinates and I translate them into ecliptic coordinates.
date right ascension declination ecliptic longitude
4 January 18h59mn1s 284°45'15" -22°44'43" -76°24'58" or -76,4162°
5 January 19h3mn24s 285°51' -22°38'18" -75°23'58" or -75,3996°
3 July 6h48mn 102° 22°58'35" 101°2'7" or 101,0355°
4 July 6h52mn8s 103°02' 22°53'39" 101°59'26" or 101,9907°
This translates as a speed of 1.0166 degree per day in January at perigee (when
in a geocentric model, that is perihelion in an heliocentric model) and a speed
of 0.9552 degree per day in July at apogee (or aphelion). Values are respectively
0.0423°/h and 0.0398°/h.
=head4 Equation of Time
The Earth's spin velocity is constant, that is 360.9856 degrees per day
but the Sun's orbital speed around the Earth is not. The combination of
both speeds is variable and it is I<not> 360 degrees per day. The crossing
of the meridian by the Sun is not exactly every 86400 seconds.
There is a difference between the Solar I<Mean> Time, where noon occurs
every 86400 seconds, no more, no less, and the Solar I<Real> Time, in which
noon is defined by the time when the Sun crosses the meridian. The difference
between the Solar Mean Time and the Solar Real Time is called I<equation
of time>.
Here are the extreme values for the equation of time in 2017, computed
by a script using L<DateTime::Event::Sunrise> and refined with Stellarium.
Date DT::E::S Stellarium
2017-11-02 11:43:33 11:43:37 -16mn23s earliest noon value,
2017-02-10 12:14:12 12:14:14 +14mn14s latest noon value
2017-09-11 11:56:33 11:56:34 -3mn26s
2017-09-12 11:56:11 11:56:13 -3mn47s biggest decrease: 21 or 22 seconds
2017-12-17 11:56:11 11:56:14 -3mn46s
2017-12-18 11:56:41 11:56:44 -3mn16s biggest increase: 30 seconds
And here is the curve for the equation of time.
=for html
<p>
<img src='equ-time.png' alt='Curve of the equation of time during one year' />
</p>
=head4 The Analemma
The irregularity of the Sun's trajectory can be visualised by using the Local Mean Time
as a reference and pinpointing the positions of the Sun at noon in LMT.
The various Sun positions day after day build an 8-shaped curve, called I<analemma>.
=head4 Mean Sun, Virtual Homocinetic Sun
In the following, it is useful to imagine a virtual Sun which would use an constant
angular speed (either in equatorial coordinates or ecliptic coordinates, depending on
which is more convenient).
The concept of I<Mean Sun> is a virtual Sun like this, calibrated so it crosses
the meridian at 12:00 (Local Mean Time) each day, and which minimizes the difference
between the real local noon and the mean local noon.
I will also consider several "virtual homocinetic suns" or VHS (no relation with
magnetic tapes). These virtal suns are synchronised with the real Sun at some
convenient point and then move with a constant angular speed.
=head1 Computing Sunrise and Sunset
Computing sunrise and sunset consists in taking in account both the variation of
day's length and the equation of time to pinpoint when the Sun reaches the
altitude that corresponds to sunrise or sunset.
In the schema below, the variation of day's length results in a bobbing up and
down of the sinusoidal curve (and less obviously, a vertical stretch or compression
of this curve). The equation of time results in a leftward or rightward shift
of the curve.
=for html
<p>
<img src='pseudo-analemma.gif' alt="Evolution of the Sun's trajectory during a year" />
</p>
Q: Wahoo! Impressive!
A: You should not be impressed. I took some liberties with the reality. First,
I figured the Sun's trajectory as a sinusoidal curve, because it is easy to compute,
but I did not check whether it was the real curve. And I would bet that it is
only approximately close to the real curve. Second, the equation of time is very much
increased. Instead of a true solar noon varying from 11:43 to 12:15 (in mean solar time),
here the variation is multiplied by 4 and the solar noon varies from 11:00, or even less,
to 13:00. But without this stretching, you would not have seen anything.
Q: And this figure eight, is this the analemma?
A: No. The analemma gives the position of the Sun as azimuth and height at
I<mean> solar noon. In the curve above, the abscisse is the mean time of
I<true> solar noon and the ordinate is the height of the Sun at this instant.
In other words, the analemma is based on a regular temporal event, the mean solar
noon, and plots the correlation between two variable spatial phenomena, the
azimuth and the height of the Sun. On the other hand, the curve above is based on
a precise spatial event, the azimuth 180°, and plots the correlation between a
variable spatial phenomenon, the height of the Sun and a variable temporal event,
the true solar noon.
I admit that the ordinates of both curves are very similar notions, and it would be
comparing golden apples with Granny Smiths. On the other hand, the abscisses are
a spatial angle in one case and a time of the day in the other case, so it would
be comparing apples with oranges.
Q: And the similarity of the shapes is juste a coincidence?
A: No, this is no coincidence. Let us start with the ordinates.
The curve above, which I will call "pseudo-analemma", gives the
height of the Sun at I<true> solar noon, so the Sun is at its highest for the current day.
Therefore, the height of the Sun in the analemma is obviously lower than on the
pseudo-analemma (except during the 4 days when the curve crosses the Y-axis).
But since we are near a point with an horizontal tangent, the variation
is very small. For example, we consider an observer at Greenwich observatory on
2nd November 2017. At true solar noon (11h 43mn 37s), the Sun is at 23°37'39"
while a quarter of an hour later, at mean solar noon it is at 23°31'40"
(values given by Stellarium).
For the abscisses, it is a bit more complicated. Let us use the same example as
above. At mean solar noon, the Sun's azimuth is 184°19'08", so on the analemma
doc/astronomical-notes.pod view on Meta::CPAN
Then we apply both Earth's spin (360.9856 degrees per day) and
the movement of a VHS ("virtual homocinetic sun"), that is, 0.9856 degrees per
day. The result is a combined rotational speed of 360 degrees per day, that is,
15 degrees per hour. And sunset happens when the VHS reaches the target altitude.
So on 4th Jan, the angle between noon and sunset is 59.9746° (59° 58' 28").
We need 3.9983 hours (3 h 59 mn 53 s) to run this angle and the sunset for
the VHS occurs at 16:04:50.
=head2 Implementation of Precise Algorithm
With the precise algorithm, we keep separate Earth's spin (360.9856 degrees
per day) and the Sun's rotational speed around Earth. In addition, this
rotational speed is the real speed, spanning from 0.9552 to 1.0166
degree per day.
First iteration. We start from the true solar noon at 12:04:56 and we apply Earth's spin
(15.04107 degrees per hour). The first result, a very approximate one, is the
instant when Earth's spin brings the Sun to the target altitude.
On 4th Jan, this first value is 16:04:11.
For iteration 2, we determine the virtual solar noon that corresponds to the
Sun's position at 16:04:11. This virtual solar noon occurs at 12:05:01.
With this reference, we apply the Earth's rotation and we get a second
value for sunset, 16:04:23 (16.0731615074431 in decimal hours).
For iteration 3, we determine the virtual solar noon
that corresponds to the Sun's position at 16:04:23. This new virtual solar
noon occurs at 12:05:01. And one more time we apply the Earth's rotation
and we obtain a third value for sunset, 16:04:23, differing from the previous
value by less than a second: 16.0731642391519 instead of 16.0731615074431.
The difference is 2,73e-6 hours, that is 9 ms, so we leave the computation loop.
This is one way to describe the algorithm: the Sun reaches by anticipation
its position in the evening and stays there, waiting for the spinning Earth
to spin until the Sun disappears below the horizon. Another way to describe
the algorithm is as follows:
During iteration 2, between the real solar noon 12:04:56 and the time given
by iteration 1, 16:04:11, the Sun orbitates with its real speed of 1.0166 degree per day
while the Earth spins at 360.9856 degrees per day.
Then, at 16:04:11, the Sun freezes in its track and after that, we adjust the
position with only the Earth's spin to reach the required altitude.
And the sunset occurs at 16:04:23.
During iteration 3, between the real solar noon 12:04:56 and the time given
by iteration 2, 16:04:23, the Sun orbitates with its real speed 1.0166 degree per day.
Then at 16:04:23, it freezes, letting the Earth continue its spin. And sunset
happens 9 milliseconds later, at 16:04:23. So, there are 3 h 59 mn 27 s when we use
the Sun's real orbital speed and 9 milliseconds when we use an obviously wrong orbital
speed. In the end, it is better than the basic algorithm, which uses an approximate
but still wrong orbital speed, but for the whole span of 3 h, 59 mn and 57 s.
=head2 What Happened in Spring 2020?
Let us look a bit farther into the past. In January 2019, I published
version 0.98 of L<Astro::Sunrise>, with the precise algorithm
implemented as explained in the preceding paragraph. For test data, I
did not know how to cross-check with authoritative sources, so I just
checked they looked plausible and that the iterative computation
stopped after a few iterations. At this time, I did not synchronise
L<DateTime::Event::Sunrise> with L<Astro::Sunrise>, because it was not
the proper time yet.
Then in April 2020, a user created an RT ticket, explaining that he
had compared the results of L<DateTime::Event::Sunrise> with
authoritative websites and that the results were not precise enough.
Of course the results were not precise, I had not yet synchronised
L<DateTime::Event::Sunrise> with L<Astro::Sunrise>. Yet this ticket
gave me two things: first, a website which would provide precise
day-after-day computations of sunrise, sunset and real solar noon, and
second a few round tuits to upgrade L<DateTime::Event::Sunrise> to the
same level as L<Astro::Sunrise>, with real tests values.
There was a glitch. The precise algorithm copied from
L<Astro::Sunrise> and the test data from the NOAA website were not
matching. After tweaking the algorithm and asking for advice on the
DateTime mailing-list
(L<https://www.nntp.perl.org/group/perl.datetime/2020/06/msg8241.html>),
I had to admit that using a 15-degree-per-hour spin speed was giving
better results than the normal 15.04107-degree-per-hour value. So I
coded the C<15> value in L<DateTime::Event::Sunrise> and I published
the module.
=head3 What Went Wrong?
I still think my explanations of the precise algorithm are correct. I
think that the error lies in the implementation of this algorithm.
Using as an example the values above, during iteration 2 we should use
the Sun at its 16:04:11 position and not moving. Without being
certain, I think that actually, the program computes the virtual solar
noon for 16:04:11, which is 12:05:01 and that after that it uses the
Sun at its 12:05:01 position.
How can I check this hypothesis? I need to "create" a new fixed star
at the 16:04:11 position of the Sun. And I do not know how to do that
in Stellarium. A few months later, I read the documentation for
Stellarium 0.20.1-1. In chapter 13, I found that the user can add
unofficial novas to the novas coming from the official star
catalogues. This would answer my need. The problem is that I have
Stellarium 0.18.0 and that I did not succeed in creating "my" nova in
this version. So for the moment, it is impossible to check the
implementation of the precise algorithm.
=head1 More About The Parameters
Below, I give some detailed explanations about the parameter used when
calling the module's functions. These explanations would have been too
long if they had been included in the module's POD and a casual doc reader
would have been drowned in a deluge of informations.
=head2 Choosing The Algorithm, C<precise> Parameter
Q: When should I choose the precise algorithm?
A: The short answer is "Never". The long answer is the following:
=over 4
=item *
If you want some twilight, use the basic algorithm.
=item *
If you live between the polar circles, use the basic algorithm.
=item *
( run in 1.062 second using v1.01-cache-2.11-cpan-ceb78f64989 )