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

[email protected] (Eric Jensen) Sat, 20 Jun 2020 23:27:20 -0400
Newsgroups perl.datetime
Message-ID <[email protected]>
--Apple-Mail=_BBD57F38-0309-4331-964D-FD4E15763CB6
Content-Transfer-Encoding: quoted-printable
Content-Type: text/plain;
	charset=utf-8

Hi Jean,=20

As an astronomer, this message caught my eye, so I thought I=E2=80=99d =
make a few comments.  Before I do, I want to say that (a) I appreciate =
people working 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 hope you=E2=80=99ll take it as a failure to convey the =
right tone by email rather than anything else.  Some of the background =
below is likely to be things you already know, but I=E2=80=99m trying to =
be complete.=20

That said, I=E2=80=99d encourage you *not* to make the selection of this =
parameter a user choice, simply because it can=E2=80=99t be the case =
that both 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 people who want to use your code.=20

As you note, 15 degrees per hour (where =E2=80=9Chour=E2=80=9D here is a =
mean solar hour, i.e. the time we normally use) is the speed that the =
Earth rotates compared to the average position of the Sun.  (The =E2=80=9C=
mean Sun=E2=80=9D is a convenient fictitious object that moves at a =
constant rate throughout 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 both the Earth=E2=80=99s rotation on its axis and =
its motion around the Sun - because it moves in its orbit each day, it =
has to rotate a little more than 360 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 compared 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 day=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 the =
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/365.25 is the same as the ratio of your two factors.)=20

So if you=E2=80=99re calculating the Sun=E2=80=99s position based on =
some coordinate system fixed to the (distant) sky rather than the Sun =
(which I assume must be what you=E2=80=99re doing), you need to use the =
Earth=E2=80=99s rotation speed compared to that coordinate system, which =
is the larger value (15.04107).

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 taking 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, the whole Sun is already below the =
horizon!  The exact amount of refraction depends on atmospheric pressure =
and temperature, but on average it=E2=80=99s about 0.6 degrees (which =
coincidentally is about the angular diameter of the Sun).  That slightly =
more than half a degree translates to about two minutes of time, so if =
that=E2=80=99s about the difference you=E2=80=99re seeing, that could be =
the cause. =20

The other reason could be different definitions of sunset, i.e. is it =
defined by when the top limb of the Sun is at the horizon?  Or the =
center?  Or is it some other distance below the horizon?  (You list =
0.833 degrees below, which is a combination of the average refraction =
and the Sun=E2=80=99s angular radius, i.e. it=E2=80=99s when the upper =
limb of the Sun would appear to go below the horizon, given average =
refraction.)   But other other codes might make other assumptions that =
could lead to differences, i.e. they could be separately accounting for =
refraction on top of a specified angle of the Sun below the horizon.=20

I suspect that I haven=E2=80=99t told you much (if anything!) that you =
didn=E2=80=99t know already, but maybe the background is useful for =
others reading the thread.  But the bottom line is that again I =
encourage you to dig a little deeper to try to find the source of the =
differences rather than exposing this as a user choice.=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 for 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=99s =
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 refraction - they say they assume 0.833 degrees of refraction, =
though as noted above I think of that value as including both 0.6 =
degrees of refraction *and* the angular radius of the Sun, so maybe =
there=E2=80=99s still a little ambiguity there.=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 on 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 those.  (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

Hope some of this is at least a little helpful, and I=E2=80=99m curious =
to see where you end up -=20

With best wishes,

Eric
=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


--Apple-Mail=_BBD57F38-0309-4331-964D-FD4E15763CB6
Content-Disposition: attachment;
	filename=smime.p7s
Content-Type: application/pkcs7-signature;
	name=smime.p7s
Content-Transfer-Encoding: base64

MIAGCSqGSIb3DQEHAqCAMIACAQExDzANBglghkgBZQMEAgEFADCABgkqhkiG9w0BBwEAAKCCC38w
ggT0MIID3KADAgECAhAA1UesWKfhxwQHxMKu9RHEMA0GCSqGSIb3DQEBCwUAMIGNMQswCQYDVQQG
EwJJVDEQMA4GA1UECAwHQmVyZ2FtbzEZMBcGA1UEBwwQUG9udGUgU2FuIFBpZXRybzEjMCEGA1UE
CgwaQWN0YWxpcyBTLnAuQS4vMDMzNTg1MjA5NjcxLDAqBgNVBAMMI0FjdGFsaXMgQ2xpZW50IEF1
dGhlbnRpY2F0aW9uIENBIEcyMB4XDTIwMDMxMjE1Mjc0NVoXDTIxMDMxMjE1Mjc0NVowIjEgMB4G
A1UEAwwXZWplbnNlbjFAc3dhcnRobW9yZS5lZHUwggEiMA0GCSqGSIb3DQEBAQUAA4IBDwAwggEK
AoIBAQDB1PDa63Eyh5xye599guQBS/ERuotkUO2j4bkyZ6X7zwYkOhTzuSoQl5kCfScMghOB1dGv
ImLn1S59dQ8HW/BpiEzfKjTcMw3N2VAssCiznZR3TefPFHUfTBYbd+Q81GMauYq88C9Bo1wA5Gqq
dq88oUt2/WjmJK4GK2LH2d95Ls1LAHzJbObYG0mvg91gN9bigWCiOG0p1WI6tsJyil4Q6OOQoVM+
sZVsy+4QWCQU6ZB/4v6jo54AtTdvZ8qExhMhCx0co4fLGzvpxdDKACA7GheYLWVtwz5fSVUByLqL
lwkW4vyfB3wixfdYNfixdUjOJGAOA9n0lTqOg97xIUl9AgMBAAGjggG4MIIBtDAMBgNVHRMBAf8E
AjAAMB8GA1UdIwQYMBaAFGvyjZ5owSUEH1E0V/YWXJTqTWkaMH4GCCsGAQUFBwEBBHIwcDA7Bggr
BgEFBQcwAoYvaHR0cDovL2NhY2VydC5hY3RhbGlzLml0L2NlcnRzL2FjdGFsaXMtYXV0Y2xpZzIw
MQYIKwYBBQUHMAGGJWh0dHA6Ly9vY3NwMDkuYWN0YWxpcy5pdC9WQS9BVVRIQ0wtRzIwIgYDVR0R
BBswGYEXZWplbnNlbjFAc3dhcnRobW9yZS5lZHUwRwYDVR0gBEAwPjA8BgYrgR8BGAEwMjAwBggr
BgEFBQcCARYkaHR0cHM6Ly93d3cuYWN0YWxpcy5pdC9hcmVhLWRvd25sb2FkMB0GA1UdJQQWMBQG
CCsGAQUFBwMCBggrBgEFBQcDBDBIBgNVHR8EQTA/MD2gO6A5hjdodHRwOi8vY3JsMDkuYWN0YWxp
cy5pdC9SZXBvc2l0b3J5L0FVVEhDTC1HMi9nZXRMYXN0Q1JMMB0GA1UdDgQWBBQXB5jM8JOC/eyh
Q8WPrAK90GuAvjAOBgNVHQ8BAf8EBAMCBaAwDQYJKoZIhvcNAQELBQADggEBAFkDUom87U+Do3j+
VaTcmfAfjkuBbM8SDuUfGE6SV5Lk3wSNG575xkBi45qdT0f0wo0Y9RPk/IFVKYjtW/J8kPdzEaoq
KgqeqKpeuqCZmcSDvJC3KOq99k8OdMiVfl4PT1QYgABlJgITcSVGWVu3GY2/Ti550LRYoDGSDVuF
1YC4ksFEDHiT763fTJDPh6xsm+4tSpOmUX7jmGwjxnfH9adwj5EKXdYPxqehF21GjqRQYpMefp4b
vki+IfzX/aX944twLsXJdm3ZgOi6igBAETsIEwBHMBPzIiDZNUtcR7LPqodcCj2mcVTDw9Ut3cVs
D3lD9xQBlBiT4hOGvhmgBTgwggaDMIIEa6ADAgECAhBP3hBL7ZVb3outZYfMQV7jMA0GCSqGSIb3
DQEBCwUAMGsxCzAJBgNVBAYTAklUMQ4wDAYDVQQHDAVNaWxhbjEjMCEGA1UECgwaQWN0YWxpcyBT
LnAuQS4vMDMzNTg1MjA5NjcxJzAlBgNVBAMMHkFjdGFsaXMgQXV0aGVudGljYXRpb24gUm9vdCBD
QTAeFw0xOTA5MjAwNzEyMDVaFw0zMDA5MjIxMTIyMDJaMIGNMQswCQYDVQQGEwJJVDEQMA4GA1UE
CAwHQmVyZ2FtbzEZMBcGA1UEBwwQUG9udGUgU2FuIFBpZXRybzEjMCEGA1UECgwaQWN0YWxpcyBT
LnAuQS4vMDMzNTg1MjA5NjcxLDAqBgNVBAMMI0FjdGFsaXMgQ2xpZW50IEF1dGhlbnRpY2F0aW9u
IENBIEcyMIIBIjANBgkqhkiG9w0BAQEFAAOCAQ8AMIIBCgKCAQEAt2hzetk81C/73GfKPc6UfP+J
Gc7aGmPzGUeQJ1go3CdFpsBPonREDXUDdmRCIRkTDroH30RLsTO/0hEFiYjCyvvbSVSm05sXkvfJ
XOXefNqK21fBayr4JCgMRyLVwqRYXlKI7bb42nYSm7YcXGTDmdcydmJuuqcLqFQawWiBMNRRVEi4
uW5uXBZgWGmq8NoKH/+5xGBFbf6tNTWcGhPVceResuwK155+OiH6jTW01Na8aLj7c7IAGJ0Y9e6h
iHtRthfW7SwbU7ys73a3nNXv8Kv9XNr0RvJKHoOsKqxjffew3GKQrMXIHB5tm/je3XEnIxUT8JG3
sEsk7IfF3VirSwIDAQABo4IB/jCCAfowDwYDVR0TAQH/BAUwAwEB/zAfBgNVHSMEGDAWgBRS2Ig6
yJ94Zu2J83s4cJTJAgI20DBBBggrBgEFBQcBAQQ1MDMwMQYIKwYBBQUHMAGGJWh0dHA6Ly9vY3Nw
MDUuYWN0YWxpcy5pdC9WQS9BVVRILVJPT1QwRQYDVR0gBD4wPDA6BgRVHSAAMDIwMAYIKwYBBQUH
AgEWJGh0dHBzOi8vd3d3LmFjdGFsaXMuaXQvYXJlYS1kb3dubG9hZDAnBgNVHSUEIDAeBggrBgEF
BQcDAgYIKwYBBQUHAwQGCCsGAQUFBwMJMIHjBgNVHR8EgdswgdgwgZaggZOggZCGgY1sZGFwOi8v
bGRhcDA1LmFjdGFsaXMuaXQvY24lM2RBY3RhbGlzJTIwQXV0aGVudGljYXRpb24lMjBSb290JTIw
Q0EsbyUzZEFjdGFsaXMlMjBTLnAuQS4lMmYwMzM1ODUyMDk2NyxjJTNkSVQ/Y2VydGlmaWNhdGVS
ZXZvY2F0aW9uTGlzdDtiaW5hcnkwPaA7oDmGN2h0dHA6Ly9jcmwwNS5hY3RhbGlzLml0L1JlcG9z
aXRvcnkvQVVUSC1ST09UL2dldExhc3RDUkwwHQYDVR0OBBYEFGvyjZ5owSUEH1E0V/YWXJTqTWka
MA4GA1UdDwEB/wQEAwIBBjANBgkqhkiG9w0BAQsFAAOCAgEAYES6GaKrcvsOQZpEwboVOb2dri/f
Jrcpb7GSEW9JmA+Kep4GLmp9X50Iv8EK478kwf2aAjnPnsOdiItALcIgecS1qVxN+EY+V5GCNEy4
VAsB5gzlQBmKI9P4PxLt9pnQJneCVEvDnVBMZAllIL5s3uaCiIEb8eYZqG8taOWSM1nqjoCZULcc
hXWYajBqaJg0RUOZ6f5IB0lb26HA/7EUVmh1nSVglDoUeD7elINXHph0z3if1722UydcoH4Jj3Za
Y9dtQ4wJSNhSZOzES72UkS6we/556FOGs7oeJWuQe8Rq2EeeSGmGliZKUbYo4jB/C2omMn0L4QwI
5wMNrWd2FRNUUwxMBmbJYtEaDRTQ72HPA8DnbRkvRDSJkjsToqU6ZpBlBf4s5EwrhXqFVb2rM9mG
CPDZJi7Hw3y8BYD/d3iTL6PW5UjOTSpFcnSIP4HW5PI6MTHXl+ab6ajCnvJw6E1TGLh3zJypv5CQ
8Ftm0z7MKLt5Zr2E4jojZXeZn1sUpSqidZyp9mG/LYMRmHMkthDRnDnO2tHv5+YOO4cUEbTt5Bww
E5RPjqovsnedyd5SijIK+k1MCXFLMTfERz3qUN3i/fwueXcGy4jEf2n/FvYsEY3GBHXZCMVWPffB
fbl/ITjs9Q9NG37bAEm/mg2yNq02NLjDbQIKgt9W0aBU9SsxggOpMIIDpQIBATCBojCBjTELMAkG
A1UEBhMCSVQxEDAOBgNVBAgMB0JlcmdhbW8xGTAXBgNVBAcMEFBvbnRlIFNhbiBQaWV0cm8xIzAh
BgNVBAoMGkFjdGFsaXMgUy5wLkEuLzAzMzU4NTIwOTY3MSwwKgYDVQQDDCNBY3RhbGlzIENsaWVu
dCBBdXRoZW50aWNhdGlvbiBDQSBHMgIQANVHrFin4ccEB8TCrvURxDANBglghkgBZQMEAgEFAKCC
AdcwGAYJKoZIhvcNAQkDMQsGCSqGSIb3DQEHATAcBgkqhkiG9w0BCQUxDxcNMjAwNjIxMDMyNzE5
WjAvBgkqhkiG9w0BCQQxIgQgLH7Kbu+lnlKTbWntjYPOphZ3vqo1VpBo2b6AxREd3bcwgbMGCSsG
AQQBgjcQBDGBpTCBojCBjTELMAkGA1UEBhMCSVQxEDAOBgNVBAgMB0JlcmdhbW8xGTAXBgNVBAcM
EFBvbnRlIFNhbiBQaWV0cm8xIzAhBgNVBAoMGkFjdGFsaXMgUy5wLkEuLzAzMzU4NTIwOTY3MSww
KgYDVQQDDCNBY3RhbGlzIENsaWVudCBBdXRoZW50aWNhdGlvbiBDQSBHMgIQANVHrFin4ccEB8TC
rvURxDCBtQYLKoZIhvcNAQkQAgsxgaWggaIwgY0xCzAJBgNVBAYTAklUMRAwDgYDVQQIDAdCZXJn
YW1vMRkwFwYDVQQHDBBQb250ZSBTYW4gUGlldHJvMSMwIQYDVQQKDBpBY3RhbGlzIFMucC5BLi8w
MzM1ODUyMDk2NzEsMCoGA1UEAwwjQWN0YWxpcyBDbGllbnQgQXV0aGVudGljYXRpb24gQ0EgRzIC
EADVR6xYp+HHBAfEwq71EcQwDQYJKoZIhvcNAQEBBQAEggEAAVtsV9CZLqta3VH1xnk13TlsEaQm
RzuOHY/4efdJU1VN5M2Eb/3Rm64UU3YUyQITRJK7Ukcd63SAc1nN+5GLDrNTCa9MvabeWpDuOncm
U1W4wYI2mCBzEX4wMxqxg/TvDuKpy/wr6UOYmgNDmCXQift5pm0FQ2e9hqHgGKeZr/vEhQfK2f08
tCMtVTaElMOQKAntoeifMezV764blnIK9weegdxZBIWtQrQEe2HRJW8nVmb37yrCdKR0iY8F9Hse
wZ93TuOyGnjI3Vbirsa6y6Ou+aAPoklkzvli4jmvMWFfIE5So7SVRkSwDgwCVPBpGh336Nk1UiYR
EyCIigB43AAAAAAAAA==
--Apple-Mail=_BBD57F38-0309-4331-964D-FD4E15763CB6--