10 нояб. 2010 г.

Юлианские даты

Гражданский календарь с его високосными годами и реформами крайне неудобен для расчетов. Сколько дней и часов осталось до начала мото-сезона? В какой день недели произошло памятное событие?... Вместо того, чтобы подсчитывать дни в годах и часы в сутках, удобнее пользоваться действительными числами. Я говорю, конечно, не о вычислениях в уме, а о программировании.

Непрерывный счет дней

Одним из широко известных решений этой проблемы стало так называемое "время Unix", измеряемое количеством секунд, прошедших от полуночи 1 января 1970 года. Прикладные языки программирования предоставляют собственные библиотечные средства. Помнится, в Delphi был тип TDateTime, который хранил даты и время в виде десятичных чисел; целая часть представляла количество дней, прошедших с полуночи 30 декабря 1899 года, а дробная часть -- время суток (стандарт OLE Automation). Базы данных, где проблема расчета времени весьма актуальна, делают это по-своему. Например, TIMESTAMP в MySQL 5.1 отсчитывается от 0ч 1 января 1000 года. Python умеет работать с датами в диапазоне от 1 до 9999 г.г. н.э. В языке Java... впрочем, довольно подробностей. Если нас интересуют события космического масштаба, как-то: восходы и заходы, равноденствия и солнцестояния, фазы Луны, координаты планет, наконец, интервалы времени, не ограничивающиеся нашей эрой -- ни одно из этих средств не подойдет. Что понадобится -- так это юлианские даты, система непрерывного счета времени, которую используют астрономы.

За точку отсчета астрономы приняли средний гринвичский полдень 1 января 4713 г. до н.э. Очередная юлианская дата начинается в 12ч00м UT, что на полсуток расходится с началом гражданских суток. Целая часть числа представляет собой число суток, дробная -- доли суток, прошедшие от полудня. Например, 6ч 17 февраля 1985 года соответствует юлианской дате 2446113.75.

Проблемы с точностью

В эпоху первых персональных компьютеров и программируемых калькуляторов возникали проблемы с точностью. В некоторых системах на хранение вещественных чисел отводилось всего 4 байта (то же самое, что сегодняшнее float), что ограничивало точность 6-7 цифрами. Между тем, юлианские даты часто представлены 11-12 цифрами. Скажем, 11ч. 17 августа 1938 года = JD 2429127.9583

Смещение точки отсчета

Можно принять за точку отсчета вместо даты, отстоящей от нас на 6 тысячелетий, момент поближе. Питер Даффетт-Смит в книге "Astronomy With Your Personal Computer" (Cambridge University Press, 1986) предлагает использовать 12ч GMT 31 декабря 1899 года, т.е. ближайший полдень до начала XXв. Вычисления в рамках ближайших столетий становятся точнее, поскольку на хранение переменных потребуется меньше памяти. Даты до этого момента отрицательные. Чтобы перейти от такой даты (назовем ее вслед за Даффетт-Смитом DJD) к стандартной юлианской (JD), надо прибавить к первой 2415020.

Современное восьмибайтное представление вещественных чисел (double), дает точность до 16 цифр, так что проблема, о которой идет речь, больше неактуальна, однако я намерен следовать методу из книги Даффетт-Смита. Тем более, что момент 12ч UT 31 декабря 1989 года взят не с потолка. Забегая вперед, скажу, что во множестве астрономических расчетов используется время в юлианских столетиях, прошедшее именно от этой даты (T). Формула такова:
T = (JD - 2415020) / 36525
Здесь JD -- юлианская дата. То же самое, что DJD / 36525. Делитель, как нетрудно догадаться, представляет собой количество суток в условном юлианском столетии.

Между прочим, класс DateTime из одноименного пакета Perl содержит метод jd, выдающий стандартную юлианскую дату. В базе SQLite имеется функция для ее получения.

От календаря к юлианской дате

Алгоритм пересчета календарной даты в юлианскую суров и непригляден. Его реализации на любых языках программирования
выглядят примерно одинаково. Ниже представлен ее вариант на языке Perl.

sub date2djd {
my %args = @_;
my $d = $args{date};
my $m = $args{month};
my $y = $args{year};

$y++ if $y < 0;
if ($m < 3) {
$m += 12;
$y--
}
my $b;
if (after_gregorian($y, $m, $d)) {
my $a = floor($y / 100);
$b = 2 - $a + floor($a / 4);
}
else {
$b = 0
}
my $f = 365.25 * $y;
my $c = floor( $y < 0 ? $f - 0.75 : $f ) - 694025;
my $e = floor(30.6001 * ($m + 1));
return $b + $c + $e + $d - 0.5
}
На входе у нее именованные аргументы:
  • year -- номер года, отрицательный для годов до н.э. Нулевое значение недопустимо, т.е. за -1 следует +1 вместо 0.
  • month - номер месяца (1-12)
  • day -- день месяца с долями суток. Например, 07ч 40м первого числа будет соответствовать 1 + (7 + 40 / 60.0) / 24.0
Правильность передаваемых аргументов не проверяется. Все проверки будут реализованы позже, при создании ООП-обертки.

Календарная дата относится к пролептическому григорианскому календарю. Так называется григорианский
календарь, который используется для обозначения дат, предшествовавших его введению. А применен он был впервые 15 октября 1582 года. Вспомогательная функция after_gregorian как раз и помогает определить, приходится ли заданная дата на период до или после календарной реформы

sub after_gregorian {
my ($y, $m, $d) = @_;
return 0 if $y < 1582;
return 1 if $y > 1582;
return 0 if $m < 10;
return 1 if $m > 10;
return 0 if $d <= 15;
return 1;
}
Использовать эту функцию самостоятельно, вне пролептического календаря, не имеет смысла, потому что григорианское летоисчисление не было принято везде одновременно. Наша страна перешла на "новый стиль" в 1918 году.

От юлианской даты к календарю

А вот как выглядит обратная функция, переводящая юлианскую дату в календарную.

sub djd2date {
my $djd = shift;

my $d = $djd + 0.5;
my ($f, $i) = modf($d);

if ($i > -115860) {
my $a = floor($i / 36524.25 + 9.9835726e-1) + 14;
$i += 1 + $a - floor($a / 4);
}
my $b = floor($i / 365.25 + 8.02601e-1);
my $c = $i - floor(365.25 * $b + 7.50001e-1) + 416;
my $g = floor($c / 30.6001);
my $dy = $c - floor(30.6001 * $g) + $f;
my $mn = $g - ($g > 13.5 ? 13 : 1);
my $yr = $b + ($mn < 2.5 ? 1900 : 1899 );
$yr-- if $yr < 1;

return {year => $yr, month => $mn, day => $dy}
}
На входе DJD -- число юлианских дней от начала XX века. На выходе -- хэш того же вида, что и на входе date2jd.

День недели

Система юлианских дней позволяет довольно просто вычислить день недели. Для этого достаточно:
  1. Взять юлианскую дату на начало календарных суток (0ч UT).
  2. Прибавить к ней 1.5.
  3. Найти остаток от деления результата на 7
0 будет соответствовать воскресенью, 1 -- понедельнику, 2 -- вторнику и так до 6 -- субботы.

Добавим две новые функции:
  • jd_midnight для вычисления юлианской даты на календарную полночь, исходя из заданной юлианской даты
  • weekday для вычисления дня недели по юлианской дате на начало суток.

#!/usr/bin/perl -w
use strict;
use warnings;
use POSIX qw/modf floor/;
use Readonly;

Readonly our $DJD_TO_JD => 2415020;
Readonly our @WEEKDAYS => qw/Sun Mon Tue Wed Thi Fri Sat/;

....

sub jd_midnight {
my $j = shift;
my $f = floor($j);
return $f + (abs($j - $f) > 0.5 ? 0.5 : -0.5);
}

sub weekday {
my $j = shift;
my $j0 = jd_midnight($j);
return ($j0 + 1.5) % 7
}


if (@ARGV < 3) {
print "Usage: perl julian.pl YEAR MONTH DAY";
exit(0)
}

my $djd = date2djd(year => $ARGV[0], month => $ARGV[1], day => $ARGV[2]);
my $j = $djd + $DJD_TO_JD;
printf "JD: %f\n", $j;
my $j0 = jd_midnight($j);
printf "JD at 0h: %f\n", $j0;
my $w = weekday($j);
print $WEEKDAYS[$w], "\n";
На месте многоточий в скрипте должны фигурировать функции date2djd и after_gregorian.
Важно: на входе jd_midnight и weekday стандартная юлианская дата, а не DJD. Иначе пришлось бы возиться с отрицательными значениями.

Пример: в какой день недели был запущен первый советский спутник? Известно, что это произошло 4 октября 1957 года 19ч28м34с UT.
Третьим аргументом (день) будет:((34 / 60 + 28) / 60 + 19) / 24 + 4 = 4.8115
$ perl julian.pl 1957 10 4.8115
JD: 2436116.311500
JD at 0h: 2436115.500000
Fri
Пятница, все правильно.

Применение юлианской даты для вычисления дня недели -- стрельба из пушки по воробьям. Но это только начало. Дальше речь пойдет о решении куда более интересных календарных и астрономических задач, а уж там без юлианской даты никуда.

Комментариев нет: