Re: About the next version of DateTme::Event::Sunrise

Best thread on this list in quite a while, thank you so much, both Jean and E=
ric.

I=E2=80=99m interested to see any info you could provide on the decision-mak=
ing and resulting mplementation!

Thank You,
Troy

> On Jun 20, 2020, at 11:27 PM, Eric Jensen <[email protected]> wrote:=

>=20
> =EF=BB=BFHi Jean,=20
>=20
> As an astronomer, this message caught my eye, so I thought I=E2=80=99d mak=
e a few comments.  Before I do, I want to say that (a) I appreciate people w=
orking on astronomical code, which can be tricky, and (b) my goal here is to=
 be helpful rather than annoyingly pedantic.  If I appear to fail at (b), I h=
ope you=E2=80=99ll take it as a failure to convey the right tone by email ra=
ther than anything else.  Some of the background below is likely to be thing=
s you already know, but I=E2=80=99m trying to be complete.=20
>=20
> That said, I=E2=80=99d encourage you *not* to make the selection of this p=
arameter a user choice, simply because it can=E2=80=99t be the case that bot=
h values are correct - so you=E2=80=99re putting the user in the position of=
 trying to guess which is the right answer, which isn=E2=80=99t helpful to p=
eople who want to use your code.=20
>=20
> As you note, 15 degrees per hour (where =E2=80=9Chour=E2=80=9D here is a m=
ean solar hour, i.e. the time we normally use) is the speed that the Earth r=
otates compared to the average position of the Sun.  (The =E2=80=9Cmean Sun=E2=
=80=9D is a convenient fictitious object that moves at a constant rate throu=
ghout the year, whereas the real Sun doesn=E2=80=99t, due to the Earth=E2=80=
=99s slightly elliptical orbit and tilted axis.)   But that rate combines bo=
th the Earth=E2=80=99s rotation on its axis and its motion around the Sun - b=
ecause it moves in its orbit each day, it has to rotate a little more than 3=
60 degrees (relative to some distant reference frame) to get the Sun back to=
 the same position in the sky.  So the Earth=E2=80=99s rotation rate compare=
d to a distant reference frame (e.g. the stars) is 15.04107 degrees per hour=
, so it rotates 360 degrees in 23 hours and 56 minutes (a =E2=80=9Csidereal d=
ay=E2=80=9D), and rotates about 361 degrees in a mean solar day (24 hours). =
   (To easily derive that =E2=80=9Cmagic value=E2=80=9D, you can note that t=
he Earth will have completed 365.25 rotations with respect to the Sun in one=
 year, but one extra rotation with respect to the stars; the ratio 366.25/36=
5.25 is the same as the ratio of your two factors.)=20
>=20
> So if you=E2=80=99re calculating the Sun=E2=80=99s position based on some c=
oordinate system fixed to the (distant) sky rather than the Sun (which I ass=
ume must be what you=E2=80=99re doing), you need to use the Earth=E2=80=99s r=
otation speed compared to that coordinate system, which is the larger value (=
15.04107).
>=20
> The question still remains of why you=E2=80=99re getting different answers=
 from different other codes you=E2=80=99re comparing to.  While it=E2=80=99s=
 possible that some are simply wrong, my best guess would be that some are t=
aking into account refraction by the Earth=E2=80=99s atmosphere and some are=
 not.  Because the path of a light ray is bent as it travels through air, it=
 turns out that when you see the bottom limb of the Sun touch the horizon, t=
he whole Sun is already below the horizon!  The exact amount of refraction d=
epends on atmospheric pressure and temperature, but on average it=E2=80=99s a=
bout 0.6 degrees (which coincidentally is about the angular diameter of the S=
un).  That slightly more than half a degree translates to about two minutes o=
f time, so if that=E2=80=99s about the difference you=E2=80=99re seeing, tha=
t could be the cause. =20
>=20
> The other reason could be different definitions of sunset, i.e. is it defi=
ned by when the top limb of the Sun is at the horizon?  Or the center?  Or i=
s it some other distance below the horizon?  (You list 0.833 degrees below, w=
hich is a combination of the average refraction and the Sun=E2=80=99s angula=
r radius, i.e. it=E2=80=99s when the upper limb of the Sun would appear to g=
o below the horizon, given average refraction.)   But other other codes migh=
t make other assumptions that could lead to differences, i.e. they could be s=
eparately accounting for refraction on top of a specified angle of the Sun b=
elow the horizon.=20
>=20
> I suspect that I haven=E2=80=99t told you much (if anything!) that you did=
n=E2=80=99t know already, but maybe the background is useful for others read=
ing the thread.  But the bottom line is that again I encourage you to dig a l=
ittle deeper to try to find the source of the differences rather than exposi=
ng this as a user choice.=20
>=20
> As for authoritative sources to compare to, as Bill Ricker notes, the USNO=
 is generally my go-to source.  Unfortunately, their website has been down f=
or months as they reorganize, and unfortunately looks like it will still be a=
 while.  (https://datacenter.iers.org/data/2/message_400.txt - that=E2=80=99=
s about the IERS but the same is true for the astronomical parts of the USNO=
 website.)  So without access to USNO, I=E2=80=99d be inclined to compare to=
 NOAA, who at least are clear on their website that they are including refra=
ction - they say they assume 0.833 degrees of refraction, though as noted ab=
ove I think of that value as including both 0.6 degrees of refraction *and* t=
he angular radius of the Sun, so maybe there=E2=80=99s still a little ambigu=
ity there.=20
>=20
> If you want to compare against a different Perl solution, in my own astro-=
related code I make use of Astro::Coords and Astro::Telescope, which depend o=
n Astro::PAL, a Perl interface to a widely-used positional astronomy library=
.  I=E2=80=99d be curious to hear what you find if you compare your code to t=
hose.  (Or perhaps I can try to make a comparison later, but I need to go to=
 bed now, and this message is too long already!)=20
>=20
> Hope some of this is at least a little helpful, and I=E2=80=99m curious to=
 see where you end up -=20
>=20
> With best wishes,
>=20
> Eric
>=20
>=20
>> On Jun 19, 2020, at 12:03 PM, Jean Forget <[email protected]> wrote:
>>=20
>> Following the creation of RT ticket
>> https://rt.cpan.org/Public/Bug/Display.html?id=3D132336
>> I have decided to implement the proper precise algorithm in
>> DateTime::Event::Sunrise (already implemented in Astro::Sunrise) and
>> to check the values with NOAA's solar calculator and with Stellarium.
>>=20
>> After I have implemented the proper precise algorithm, there are still
>> differences between the values obtained by Stellarium and the NOAA on
>> one side and by DT::E::S on the other side. I have found that these
>> differences are the caused by the "magic" values 15 and 15.04107 used
>> in the module.
>>=20
>> 15.04107 is the angular speed of Earth's spin, while 15 is the
>> combination of Earth's spin with the orbital speed of the
>> Sun-Earth system. According to Paul Schlyter 's description of
>> the precise algorithm, we should use 15.04107, but I have found
>> that using 15 instead gives results closer to Stellarium.
>>=20
>> More abstractly, the problem is: when you are an low-level
>> amateur astronomer, who is an authoritative source for you?
>> Whom do you trust? Paul Schlyter or the Stellarium team or
>> the NOAA?
>>=20
>> I have decided to defer this question to the module user. There is a
>> new parameter, "algo", which lets the user choose between trusting
>> Paul Schlyter more than Stellarium or trusting Stellarium more
>> than Paul Schlyter.
>>=20
>> Now you may code either
>>=20
>> my $sunriseset =3D DateTime::Event::Sunrise->new(
>>                    latitude  =3D> 42.3358,
>>                    longitude =3D> -71.2044,
>>                    altitude  =3D> '-0.833',
>>                    precise   =3D> 1,
>>                    algo      =3D> 'Schlyter');
>>=20
>> or
>>=20
>> my $sunriseset =3D DateTime::Event::Sunrise->new(
>>                    latitude  =3D> 42.3358,
>>                    longitude =3D> -71.2044,
>>                    altitude  =3D> '-0.833',
>>                    precise   =3D> 1,
>>                    algo      =3D> 'Stellarium');
>>=20
>> Can you think of a better name for the new option? The creator of
>> ticket 132336 suggests
>>         spin =3D> 15.04107,
>> or
>>         spin =3D> 15,
>>=20
>> The problem is that the "magic" value 15.04107 is not easy to remember,
>> especially if the module user knows even less astronomy than myself.
>> While "Schlyter" and "Stellarium" are easier to remember.
>>=20
>> I will release the new version of DT::E::S to CPAN next week.
>> If you have a better idea for the name of the new parameter, please
>> tell me before this release.
>>=20
>> In the meantime, the updated code is available on Github.
>>=20
>> Thank you for your advice,
>>=20
>> By the way, this is the second attempt to write this message
>> to the datetime mailing list. As it happened already, my first
>> attempt was lost somewhere and did not reach you or the websites
>> https://www.nntp.perl.org/group/perl.datetime/
>> and https://www.mail-archive.com/[email protected]/.
>>=20
>> Jean Forget
>=20