[email protected] wrote:
>Perhaps, we should evaluate some of the leap-second estimation algorithms.
If you like. To be clear, I'm not pushing for the conversion to use an
estimation strategy per se, and we're now going beyond what's necessary
to address my original bug report. Documenting the existing behaviour
resolves the bug that I reported qua bug. I do still reckon the existing
behaviour sucks, but with it as a documented API we're in the realm
of differing judgements on a language design question, rather than a
clear bug.
The preferred outcome for which I was pushing was to have either (or
both, in separate methods) of the behaviours that I consider sensible:
signalling an error or making a reasonable estimate. You've made clear
that you're not at all a fan of estimation, and that's fine. It's totally
compatible with my preferences, if one then accepts the conclusion that
the conversion should signal an error. But I'm getting the impression
that you don't find erroring very palatable either.
The starting point for leap estimation is the tidal braking effect by
which the Moon is gradually slowing the rotation of the Earth, affecting
the UT1<->TT relationship. This is a long-term secular change, which
therefore must be taken into account in order to make any reasonable
estimate any significant number of years beyond one's present knowledge.
There are several other effects on the Earth's rotation which affect
UT1<->TT, but they are all oscillations (on periods of a day up to
decades), not secular drift. They can therefore be ignored, at least
for an initial version and for our purposes quite likely forever, even
though on the decadal timescale they swamp the tidal braking effect.
To qualify as a reasonable estimate of UT1<->TT it is both necessary
and sufficient to account for recent length of day and tidal braking.
<http://www.ucolick.org/~sla/leapsecs/deltat.html> has some nice plots
showing differences between time scales. The first plot, with the
3000-year span of UT1<->TT, is the most relevant to our situation.
The roughly-quadratic curve there is what needs to be extrapolated.
The way this historical information has been determined over such a span,
extending way before mechanical time measurement, is pretty clever stuff:
a written record that a solar eclipse was visible from a particular
geographic location tells you which way the Earth was pointing (UT1)
at a time (TT) that can be precisely determined by orbital calculations.
There have been many academic attempts to model and extrapolate this
curve. They differ largely in how closely they attempt to model
the last couple of centuries for which we have much more precise
measurements. Any attempt to model variations on such a short
timescale necessarily ends up modelling some of the oscillations,
not just the secular trend, so ends up a lot more complicated.
<http://www.ucolick.org/~sla/leapsecs/future2100.pdf> is a nice plot
comparing a variety of models against the eclipse observations. As you
can see, there's quite a bit of disagreement between the models, and none
of them is a great fit to the observations. But the models are of value:
there's a rough agreement if one ignores the linear ones.
Excluding the linear models, it looks like none of them is compellingly
superior to another for our purpose. Let's therefore take the simplest
class of these models: a pure linear increase in length of day, giving
a pure parabola of delta-T. Middle-of-the-road values to use are a
LOD increase of 1.7 ms per century (astronomers use the mean Julian
century, 36525 days), with LOD exactly equal to 86400 s at the year 1820.
(1820 is the midpoint of the 1750-to-1890 span of the observations behind
Simon Newcomb's theory of the planetary orbits, which is indirectly what
fixed the length of the second as a modern unit of measurement.)
To fill out our model of projected UTC, let's presume that at the
threshold date (the date of the next possible leap second not yet
scheduled) UT1=UTC, then we'll graft onto that the remaining portion of
the delta-T parabola. That gives us a model of TAI-UT1 for the future.
Then let's suppose that each leap second happens at the end of the
Gregorian month in which the fractional part of TAI-UT1 crosses 0.5.
(Current practice is for leaps to happen only in June and December,
but the rules allow the end of any month. Starting in the 38th century
we require more than one leap per month; it's anybody's guess how UTC
will actually be managed then, so it's not too bad to model it as a
multi-second leap at the end of the month.)
The attached program implements this model, for bidirectional conversions
between TAI and UTC. There is... an amount of support code. I had to
import a bunch of fundamental Gregorian calendar stuff imitating my Perl
5 module Date::ISO8601. The actual leap second logic is only 80 lines
in the middle of the 400 line file, and that covers exact conversions
for the known schedule as well as the estimation for the unknown future.
Each of the two regimes takes about half of the 80 lines. There's then a
bunch of ISO 8601 text formatting and parsing code, which doesn't support
the conversions themselves but is only used for the testing interface.
Invoke like this:
$ perl6 utc_estimate.pl6 '2016-07-28T02:26:33 UTC' '2016-12-01T00:00:20.123 TAI'
2016-07-28T02:27:09 TAI = 2016-07-28T02:26:33 UTC
2016-12-01T00:00:20.123000 TAI = 2016-11-30T23:59:44.123000 UTC
Errors are checked everywhere they should be, but the error messages are
not awesome. Conversions in both directions tick correctly through leap
seconds, both real ones and guessed future ones:
2015-07-01T00:00:34 TAI = 2015-06-30T23:59:59 UTC
2015-07-01T00:00:35 TAI = 2015-06-30T23:59:60 UTC
2015-07-01T00:00:36 TAI = 2015-07-01T00:00:00 UTC
2017-12-01T00:00:36 TAI = 2017-11-30T23:59:59 UTC
2017-12-01T00:00:37 TAI = 2017-11-30T23:59:60 UTC
2017-12-01T00:00:38 TAI = 2017-12-01T00:00:00 UTC
That was satisfying to write.
>There are
>more leap second problems in Rakudo besides the .from-posix,
Sure. Anything else using the leap second table runs into the same
issues in some form.
>that I now
>opened in
>https://rt.perl.org/Ticket/Display.html?id=128752
As written, that ticket is about a more specific idea of exposing the leap
second table explicitly. On its own that doesn't address the analogous
issues for other uses of the table. I don't have very much opinion
about exposing the table per se. If you do expose the current table,
you'd probably want to expose a threshold-of-the-unknown date as well,
because there's more to leap schedule knowledge than just the dates of
actual leaps. If you expose it in writable form, I'd recommend writing
via method rather than via lvalue, to avoid tying yourself to the table's
current format.
-zefram
use v6;
# classes for representing times on particular time scales
class TaiTime {
has FatRat:D $.linear1958 = 0.FatRat;
method from_linear1958(TaiTime:U: FatRat:D $linear1958) {
self.new(:$linear1958)
}
multi method perl(TaiTime:D:) {
"{self.perl}.from_linear1958({$!linear1958.perl})"
}
method mjd(TaiTime:D:) { 36204 + ($!linear1958 / 86400) }
method from_mjd(TaiTime:U: FatRat:D $mjd) {
self.from_linear1958(($mjd - 36204) * 86400)
}
method mjdn(TaiTime:D:) { self.mjd.truncate }
method mjdf(TaiTime:D:) { self.mjd % 1 }
method mjdnf(TaiTime:D:) {
my $mjd = self.mjd;
return ($mjd.truncate, $mjd % 1);
}
method from_mjdnf(TaiTime:U: (Int:D $mjdn,
FatRat:D $mjdf where $mjdf >= 0 && $mjdf < 1)) {
self.from_mjd($mjdn + $mjdf)
}
}
class UtcTime {
has Int:D $.mjdn = 0;
has FatRat:D $.mjdf = 0.FatRat;
method mjdnf(UtcTime:D:) { ($!mjdn, $!mjdf) }
method from_mjdnf(UtcTime:U: (Int:D $mjdn,
FatRat:D $mjdf where $mjdf >= 0)) {
self.new(:$mjdn, :$mjdf)
}
multi method perl(UtcTime:D:) {
"{self.perl}.from_mjdnf(({$!mjdn.perl}, {$!mjdf.perl}))"
}
method mjd(UtcTime:D:) { $!mjdn + $!mjdf }
method from_mjd(UtcTime:U: FatRat:D $mjd) {
self.from_mjdnf(($mjd.truncate, $mjd % 1))
}
}
# Gregorian calendar and 24-hour clock
sub year_leap(Int:D $y) { $y %% 4 && ($y !%% 100 || $y %% 400) }
sub year_days(Int:D $y) { year_leap($y) ?? 366 !! 365 }
my @month_length = (31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31);
sub month_days(Int:D $y, Int:D $m) {
if $m == 2 {
return year_leap($y) ?? 29 !! 28;
} else {
return @month_length[$m - 1];
}
}
my @nonleap_monthstarts =
(0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365);
my @leap_monthstarts =
(0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366);
sub year_monthstarts(Int:D $y) {
year_leap($y) ?? @leap_monthstarts !! @nonleap_monthstarts
}
sub ymd_from_mjdn(Int:D $mjdn) {
my $d = $mjdn - -678941;
my $qc = $d div (365 * 400 + 97);
$d %= (365 * 400 + 97);
my $y = $d div 366;
my $leaps = ($y + 3) div 4;
$leaps -= ($leaps - 1) div 25 unless $leaps == 0;
$d -= 365 * $y + $leaps;
my $yd = year_days($y);
if $d >= $yd {
$d -= $yd;
$y++;
}
$d++;
my @monthstarts = year_monthstarts($y);
my $m = 1;
while $d > @monthstarts[$m] { $m++; }
return ($qc * 400 + $y, $m, $d - @monthstarts[$m - 1]);
}
sub ymd_to_mjdn((Int:D $y, Int:D $m where $m >= 1 && $m <= 12,
Int:D $d where $d >= 1)) {
my @monthstarts = year_monthstarts($y);
my $md = @monthstarts[$m] - @monthstarts[$m - 1];
$d <= $md or die "day number out of range";
my $dd = @monthstarts[$m - 1] + $d - 1;
my $qc = $y div 400;
my $yy = $y % 400;
my $leaps = ($yy + 3) div 4;
$leaps -= ($leaps - 1) div 25 unless $leaps == 0;
return -678941 + (365 * 400 + 97) * $qc + 365 * $yy + $leaps + $dd;
}
sub hms_from_mjdf(FatRat:D $mjdf where $mjdf >= 0) {
my $hf = $mjdf * 24;
my $h = $hf.truncate;
if $h >= 24 {
return (23, 59, ($mjdf * 86400) - 86340);
} else {
my $mf = ($hf % 1) * 60;
my $m = $mf.truncate;
return ($h, $m, ($mf % 1) * 60);
}
}
sub hms_to_mjdf((Int:D $h where $h >= 0 && $h <= 23,
Int:D $m where $m >= 0 && $m <= 59,
FatRat:D $s where $s >= 0)) {
($h == 23 && $m == 59) || $s < 60 or die "seconds out of range";
return ($h * 60 + $m).FatRat / 1440 + $s / 86400;
}
# leap second processing
my @known_leaps = <
1972-06 1972-12 1973-12 1974-12 1975-12 1976-12 1977-12 1978-12 1979-12
1981-06 1982-06 1983-06 1985-06 1987-12 1989-12 1990-12 1992-06 1993-06
1994-06 1995-12 1997-06 1998-12 2005-12 2008-12 2012-06 2015-06 2016-12
>;
my $next_possible_leap = "2017-06";
sub leap_month_to_mjdn(Str:D $str) {
(my $y, my $m) = $str.split("-").map({ .Int });
$m++;
if $m == 13 { $y++; $m = 1; }
return ymd_to_mjdn(($y, $m, 1));
}
sub leap_month_to_tai(Str:D $str, Int:D $dtai) {
return TaiTime.from_mjd(leap_month_to_mjdn($str) +
FatRat.new($dtai, 86400));
}
my @utc_segs = ({
start_utc_mjdn => leap_month_to_mjdn("1971-12"),
start_tai => leap_month_to_tai("1971-12", 10).linear1958,
dtai => 10,
},);
for @known_leaps {
my $dtai = @utc_segs[*-1]<dtai> + 1;
my $bound = leap_month_to_tai($_, $dtai).linear1958;
@utc_segs[*-1]<end_tai> = $bound;
my $end_utc_mjdn = leap_month_to_mjdn($_);
@utc_segs[*-1]<end_utc_mjdn> = $end_utc_mjdn;
@utc_segs.push: {
start_utc_mjdn => $end_utc_mjdn,
start_tai => $bound,
dtai => $dtai,
};
}
@utc_segs[*-1]<end_tai> =
leap_month_to_tai($next_possible_leap, @utc_segs[*-1]<dtai>).linear1958;
@utc_segs[*-1]<end_utc_mjdn> = leap_month_to_mjdn($next_possible_leap);
my $lod_nominal_time = leap_month_to_tai("1819-12", -20);
my $lod_increase_rate = 0.0017;
my $xlod_at_threshold = $lod_increase_rate *
(@utc_segs[*-1]<end_tai> - $lod_nominal_time.linear1958) /
(86400*36525);
sub tai-ut1_for_tai(FatRat:D $lin) {
my $tdiffdays = ($lin - @utc_segs[*-1]<end_tai>) / 86400;
return @utc_segs[*-1]<dtai> + $xlod_at_threshold * $tdiffdays +
$lod_increase_rate * $tdiffdays * $tdiffdays / (2 * 36525);
}
sub dtai_for_month(Int:D $mjdn) {
my $dtai = 0;
loop {
my $lin = TaiTime.from_mjdnf(($mjdn, 0.FatRat)).linear1958 +
$dtai;
my $tai-ut1 = tai-ut1_for_tai($lin);
my $new_dtai = ($tai-ut1 + 0.5).truncate;
return $dtai if $new_dtai == $dtai;
$dtai = $new_dtai;
}
}
sub mjdn_start_month(Int:D $mjdn) {
my $ymd = ymd_from_mjdn($mjdn);
return ymd_to_mjdn(($ymd[0], $ymd[1], 1));
}
sub mjdn_next_month(Int:D $mjdn) {
my $ymd = ymd_from_mjdn($mjdn);
return $mjdn + month_days($ymd[0], $ymd[1]);
}
sub utc_from_tai_best_guess(TaiTime:D $tai) {
my $lin = $tai.linear1958;
if $lin < @utc_segs[0]<start_tai> {
die "can't handle time prior to leap-seconds UTC";
} elsif $lin < @utc_segs[*-1]<end_tai> {
my $l = 0;
my $r = @utc_segs.elems - 1;
until $r == $l {
my $t = ($l + $r) +> 1;
if $lin < @utc_segs[$t]<end_tai> {
$r = $t;
} else {
$l = $t + 1;
}
}
my $dtai = @utc_segs[$l]<dtai>;
my $utc_mjd = $tai.mjd - FatRat.new($dtai, 86400);
my $utc_mjdnf = ($utc_mjd.truncate, $utc_mjd % 1);
if $utc_mjdnf[0] == @utc_segs[$l]<end_utc_mjdn> {
$utc_mjdnf = ($utc_mjdnf[0] - 1, $utc_mjdnf[1] + 1);
}
return UtcTime.from_mjdnf($utc_mjdnf);
} else {
my $tai-ut1 = tai-ut1_for_tai($lin);
my $rough_mjd = $tai.mjd -
FatRat.new(($tai-ut1 + 1).truncate, 86400);
my $rough_mjdn = $rough_mjd.truncate;
my $m0_mjdn = mjdn_start_month($rough_mjdn);
my $m0_dtai = dtai_for_month($m0_mjdn);
my $m0_start = TaiTime.from_mjd($m0_mjdn +
FatRat.new($m0_dtai, 86400)).linear1958;
loop {
my $m1_mjdn = mjdn_next_month($m0_mjdn);
my $m1_dtai = dtai_for_month($m1_mjdn);
my $m1_start = TaiTime.from_mjd($m1_mjdn +
FatRat.new($m1_dtai, 86400)).linear1958;
if $lin < $m1_start {
my $utc_mjd = $tai.mjd -
FatRat.new($m0_dtai, 86400);
my $utc_mjdnf =
($utc_mjd.truncate, $utc_mjd % 1);
while $utc_mjdnf[0] >= $m1_mjdn {
$utc_mjdnf = ($utc_mjdnf[0] - 1,
$utc_mjdnf[1] + 1);
}
return UtcTime.from_mjdnf($utc_mjdnf);
}
$m0_mjdn = $m1_mjdn;
$m0_dtai = $m1_dtai;
$m0_start = $m1_start;
}
}
}
sub utc_to_tai_best_guess(UtcTime:D $utc) {
my $mjdnf = $utc.mjdnf;
if $mjdnf[0] < @utc_segs[0]<start_utc_mjdn> {
die "can't handle time prior to leap-seconds UTC";
} elsif $mjdnf[0] < @utc_segs[*-1]<end_utc_mjdn> {
my $l = 0;
my $r = @utc_segs.elems - 1;
until $r == $l {
my $t = ($l + $r) +> 1;
if $mjdnf[0] < @utc_segs[$t]<end_utc_mjdn> {
$r = $t;
} else {
$l = $t + 1;
}
}
my $dtai = @utc_segs[$l]<dtai>;
if $mjdnf[1] >= 1 &&
$mjdnf[0] != @utc_segs[$l]<end_utc_mjdn> - 1 {
die "given UTC time does not exist";
}
my $tai = TaiTime.from_mjd($mjdnf[0] +
FatRat.new($dtai, 86400) + $mjdnf[1]);
$tai.linear1958 < @utc_segs[$l]<end_tai>
or die "given UTC time does not exist";
return $tai;
} else {
my $m0_mjdn = mjdn_start_month($mjdnf[0]);
my $m0_dtai = dtai_for_month($m0_mjdn);
my $m1_mjdn = mjdn_next_month($m0_mjdn);
if $mjdnf[1] >= 1 && $mjdnf[0] != $m1_mjdn - 1 {
die "given UTC time does not exist";
}
my $m1_dtai = dtai_for_month($m1_mjdn);
my $m1_start = TaiTime.from_mjd($m1_mjdn +
FatRat.new($m1_dtai, 86400)).linear1958;
my $tai = TaiTime.from_mjd($mjdnf[0] +
FatRat.new($m0_dtai, 86400) + $mjdnf[1]);
$tai.linear1958 < $m1_start
or die "given UTC time does not exist";
return $tai;
}
}
# ISO 8601 string format
my regex decdig { <[0123456789]> }
sub fatrat_from_str(Str:D $str) {
/^(<decdig>+)(\.(<decdig>+))?$/.ACCEPTS($str)
or die "malformed numeric string";
return FatRat.new(($0.Str ~ ($1 =:= Nil ?? "" !! $1[0].Str)).Int,
10 ** ($1 =:= Nil ?? 0 !! $1[0].chars));
}
sub iso8601_from_y(Int:D $y) {
sprintf($y < 0 || $y > 9999 ?? "%+05d" !! "%04d", $y)
}
sub iso8601_to_y(Str:D $str) {
/^ ( <[-+]> <decdig>**4..* || <decdig>**4 ) $/.ACCEPTS($str) &&
!/^\-0+$/.ACCEPTS($str)
or die "malformed year string";
return $str.Int;
}
sub iso8601_from_ymd((Int:D $y, Int:D $m, Int:D $d)) {
iso8601_from_y($y) ~ sprintf("-%02d-%02d", $m, $d)
}
sub iso8601_to_ymd(Str:D $str) {
/^(<[-+]>?<decdig>+)\-(<decdig>**2)\-(<decdig>**2)$/.ACCEPTS($str)
or die "malformed date string";
return (iso8601_to_y($0.Str), $1.Int, $2.Int);
}
sub iso8601_from_hms((Int:D $h, Int:D $m, FatRat:D $s)) {
my $res = sprintf("%02d:%02d:%02d", $h, $m, $s.truncate);
my $f = $s % 1;
if $f != 0 {
my $us = $f * 1000000;
$res ~= sprintf(".%06d", $us.truncate);
$us %% 1 or $res ~= "~";
}
return $res;
}
sub iso8601_to_hms(Str:D $str) {
/^(<decdig>**2)\:(<decdig>**2)\:(<decdig>**2 (\.<decdig>+)?)$/\
.ACCEPTS($str)
or die "malformed time string";
return ($0.Int, $1.Int, fatrat_from_str($2.Str));
}
sub iso8601_from_mjdn(Int:D $mjdn) { iso8601_from_ymd(ymd_from_mjdn($mjdn)) }
sub iso8601_to_mjdn(Str:D $str) { ymd_to_mjdn(iso8601_to_ymd($str)) }
sub iso8601_from_mjdf(FatRat:D $mjdf where $mjdf >= 0) {
iso8601_from_hms(hms_from_mjdf($mjdf))
}
sub iso8601_to_mjdf(Str:D $str) { hms_to_mjdf(iso8601_to_hms($str)) }
sub iso8601_from_mjdnf((Int:D $mjdn, FatRat:D $mjdf where $mjdf >= 0)) {
iso8601_from_mjdn($mjdn) ~ "T" ~ iso8601_from_mjdf($mjdf)
}
sub iso8601_to_mjdnf(Str:D $str) {
/^ ((<decdig>||<[-+]>)+) <[tT]> ((<decdig>||<[:.]>)+) $/.ACCEPTS($str)
or die "malformed date/time string";
return (iso8601_to_mjdn($0.Str), iso8601_to_mjdf($1.Str));
}
# main program for test
sub mjdnf_from_str(Str:D $str) {
if /^<decdig>+(\.<decdig>+)?$/.ACCEPTS($str) {
my $mjd = fatrat_from_str($str);
return ($mjd.truncate, $mjd % 1);
} else {
return iso8601_to_mjdnf($str);
}
}
for @*ARGS {
/^(<-[\x20]>+) ' ' (:i (tai|utc))$/.ACCEPTS($_)
or die "bad argument";
my $mjdnf = mjdnf_from_str($0.Str);
my $tai;
my $utc;
if /^(:i tai)$/.ACCEPTS($1.Str) {
$tai = TaiTime.from_mjdnf($mjdnf);
$utc = utc_from_tai_best_guess($tai);
} else {
$utc = UtcTime.from_mjdnf($mjdnf);
$tai = utc_to_tai_best_guess($utc);
}
say "{iso8601_from_mjdnf($tai.mjdnf)} TAI = " ~
"{iso8601_from_mjdnf($utc.mjdnf)} UTC";
}
