![[personal profile]](https://www.dreamwidth.org/img/silk/identity/user.png){kind=small}
waspagvВсякий, кому приходилось заниматься астрономическими расчетами, наверняка сталкивался с необходимостью переводить календарную дату в юлианскую и обратно. Для этого можной найти множество формул и алгоритмов, но все они имеют один существенный недостаток: ограниченную область дат. Как правило, они неприменимы для далекого прошлого - отрицательных юлианских дней. А необходимость забираться так далеко имеется. Например, численная эфемерида DE431 простирается в прошлое аж до 13200 г. до н.э.
Почему же имеющиеся алгоритмя врут в таких случаях? Исследование показывает, что главная причина связана с понятиями целочисленного частного и остатка. Дональд Кнут в своем фундаментальном труде "Искусство программирования" определял целочисленное частное чисел x и y как наибольшее целое, не превосходящее x/y. Иными словами, частное чисел 10 и 3 ожидаемо равно 3, но вот для чисел -10 и 3 результат должен быть -4, а вовсе не -3, потому что -4 - как раз наибольшее целое, не превосходящее -10/3=-3.333... Сообразно с этим, остаток определяется как неотрицательное число r, такое, что x=q·y+r, где q - определенное выше целочисленное частное.
Проблема в том, что популярные языки программирования этим определениям не следуют. Как правило, при вычислении частного дробная часть просто отбрасывается, и получается -10÷3=-3. Остаток тоже сплошь и рядом получается отрицательным. Из известных мне языков только Модула-2 реализует эти операции в том же смысле, как их понимает Кнут. Но стоит только корректно выразить эти операции, как отыскание календарных формул, годных для любых дат, становится простым.
Ниже приводится алгоритм вычисления юлианской даты по календарной и обратно, пригодный для любых дат. Формулы и их вывод можно найти здесь. Алгоритм написан на ныне мертвом, но очень популярном в прошлом языке Алгол-60. Для любителей живых даю две реализации: на языках Паскаль (ISO 7185) и Фортран. Для языка Паскаль нужно обратить внимание, что стандарт ISO 7185 определяет операцию mod точно в смысле Кнута, чего не скажешь про операцию div. Если ваш компилятор не следует стандарту и реализует mod иначе - выкручивайтесь самостоятельно. Любители всяких Сей и Яв, тоже пишите сами!
Вывод для тестового примера такой:
integer JulFlag;
comment Вычисление частного q и остатка r двух целых чисел с условием r≥0;
integer procedure mod(x,y);
value x,y; integer x,y;
begin
integer q,r;
if y<0 then
begin
outstring(1,‘mod: отрицательный делитель: ’); outinteger(1,y);
outstring(1,‘‘NL’’);
r := 0
end
else
begin
q := x÷y;
r := x - y×q;
if x<0 ∧ r≠0 then
r := r + y
end;
mod := r
end mod;
integer procedure div(x,y);
value x,y; integer x,y;
begin
integer q,r;
if y<0 then
begin
outstring(1,‘div: отрицательный делитель: ’); outinteger(1,y);
outstring(1,‘‘NL’’);
r := 0
end
else
begin
q := x÷y;
r := x - y×q;
if x<0 ∧ r≠0 then
q := q - 1
end;
div := q
end div;
procedure divmod(x,y,q,r);
value x,y; integer x,y,q,r;
begin
if y<0 then
begin
outstring(1,‘divmod: отрицательный делитель: ’); outinteger(1,y);
outstring(1,‘‘NL’’);
q := 0;
r := 0
end
else
begin
q := x÷y;
r := x - y×q;
if x<0 ∧ r≠0 then
begin
q := q - 1;
r := r + y
end
end
end divmod;
comment Вычисление юлианской даты по календарной. Глобальная переменная
JulFlag определяет используемый календарь: JulFlag<0 - юлианский,
JulFlag>0 - григорианский, при JulFlag=0 с 15.10.1582 года
используется григорианский календарь, иначе - юлианский.
Функция дает корректный результат для любых дат;
real procedure julday(j,m,d,h,min,s);
value j,m,d,h,min,s; integer j,m,d,h,min; real s;
begin
integer c0,j0,j1,j2,x1,x2,x3,x4,jd;
if JulFlag<0 ∨ JulFlag=0 ∧ (j<1582 ∨ (j=1582 ∧ m<10) ∨ (j=1582 ∧ m=10 ∧ d<15)) then
begin
comment юлианский календарь;
j0 := 1721117;
c0 := div(m-3,12);
j1 := div(1461×(j+c0),4);
j2 := div(153×m-1836×c0-457,5);
jd := j1 + j2 + d + j0
end
else
begin
comment григорианский календарь;
c0 := div(m-3,12);
x4 := j + c0;
divmod(x4,100,x3,x2);
x1 := m - 12×c0 - 3;
jd := div(146097×x3,4) + div(36525×x2,100) + div(153×x1+2,5) + d + 1721119
end;
julday := jd + (h+(min+s/60)/60)/24 - 0.5
end julday;
comment Вычисление календарной даты по юлианской. Глобальная переменная
JulFlag определяет используемый календарь: JulFlag<0 - юлианский,
JulFlag>0 - григорианский, при JulFlag=0 и jd>2299160 используется
григорианский календарь, иначе - юлианский.
Процедура дает корректный результат для любых юлианских дат,
включая отрицательные;
procedure caldat(jd,y,m,d,h,min,s);
value jd; real jd,s; integer y,m,d,h,min;
begin
real dt; integer j,c0,k1,k2,x1,x2,x3,r1,r2,r3,y2;
comment выделить дробную часть суток;
dt := jd + 0.5;
j := jd;
dt := (dt-j)×24.0;
if dt<0 then
begin
j := j - 1; dt := dt + 24.0
end;
h := entier(dt);
dt:= (dt-h)×60.0;
min := entier(dt);
s := (dt-min)×60.0;
if JulFlag<0 ∨ JulFlag=0 ∧ j ≤ 2299160 then
begin
comment юлианский календарь;
y2 := j - 1721118;
k2 := 4×y2 + 3;
k1 := 5×div(mod(k2,1461),4) + 2;
x1 := div(k1,153);
c0 := div(x1+2,12);
y := div(k2,1461) + c0;
m := x1 - 12×c0 + 3;
d := div(mod(k1,153),5) + 1
end
else
begin
comment григорианский календарь;
divmod(4×j-6884477,146097,x3,r3);
divmod(100×div(r3,4)+99,36525,x2,r2);
divmod(5×div(r2,100)+2,153,x1,r1);
d := div(r1,5) + 1;
c0 := div(x1+2,12);
y := 100×x3 + x2 + c0;
m := x1 - 12×c0 + 3
end
end caldat;
procedure test(y,m,d,h,min,s);
value y,m,d,h,min,s; integer y,m,d,h,min; real s;
begin
real jd;
jd := julday(y,m,d,h,min,s);
outreal(1,jd);
caldat(jd,y,m,d,h,min,s);
outinteger(1,y); outinteger(1,m); outinteger(1,d);
outinteger(1,h); outinteger(1,min); outreal(1,s);
outstring(1,‘‘NL’’)
end test;
JulFlag := 0;
test(1900,1,0,12,0,0);
test(2000,1,1,12,0,0);
test(-5000,7,12,16,6,32.87);
test(5000,7,12,16,6,32.87)
end
Почему же имеющиеся алгоритмя врут в таких случаях? Исследование показывает, что главная причина связана с понятиями целочисленного частного и остатка. Дональд Кнут в своем фундаментальном труде "Искусство программирования" определял целочисленное частное чисел x и y как наибольшее целое, не превосходящее x/y. Иными словами, частное чисел 10 и 3 ожидаемо равно 3, но вот для чисел -10 и 3 результат должен быть -4, а вовсе не -3, потому что -4 - как раз наибольшее целое, не превосходящее -10/3=-3.333... Сообразно с этим, остаток определяется как неотрицательное число r, такое, что x=q·y+r, где q - определенное выше целочисленное частное.
Проблема в том, что популярные языки программирования этим определениям не следуют. Как правило, при вычислении частного дробная часть просто отбрасывается, и получается -10÷3=-3. Остаток тоже сплошь и рядом получается отрицательным. Из известных мне языков только Модула-2 реализует эти операции в том же смысле, как их понимает Кнут. Но стоит только корректно выразить эти операции, как отыскание календарных формул, годных для любых дат, становится простым.
Ниже приводится алгоритм вычисления юлианской даты по календарной и обратно, пригодный для любых дат. Формулы и их вывод можно найти здесь. Алгоритм написан на ныне мертвом, но очень популярном в прошлом языке Алгол-60. Для любителей живых даю две реализации: на языках Паскаль (ISO 7185) и Фортран. Для языка Паскаль нужно обратить внимание, что стандарт ISO 7185 определяет операцию mod точно в смысле Кнута, чего не скажешь про операцию div. Если ваш компилятор не следует стандарту и реализует mod иначе - выкручивайтесь самостоятельно. Любители всяких Сей и Яв, тоже пишите сами!
Вывод для тестового примера такой:
2415020.000 1899 12 31 12 0 0.000000 2451545.000 2000 1 1 12 0 0.000000 -104998.829 -5000 7 12 16 6 32.870000 3547465.171 5000 7 12 16 6 32.869994begin
integer JulFlag;
comment Вычисление частного q и остатка r двух целых чисел с условием r≥0;
integer procedure mod(x,y);
value x,y; integer x,y;
begin
integer q,r;
if y<0 then
begin
outstring(1,‘mod: отрицательный делитель: ’); outinteger(1,y);
outstring(1,‘‘NL’’);
r := 0
end
else
begin
q := x÷y;
r := x - y×q;
if x<0 ∧ r≠0 then
r := r + y
end;
mod := r
end mod;
integer procedure div(x,y);
value x,y; integer x,y;
begin
integer q,r;
if y<0 then
begin
outstring(1,‘div: отрицательный делитель: ’); outinteger(1,y);
outstring(1,‘‘NL’’);
r := 0
end
else
begin
q := x÷y;
r := x - y×q;
if x<0 ∧ r≠0 then
q := q - 1
end;
div := q
end div;
procedure divmod(x,y,q,r);
value x,y; integer x,y,q,r;
begin
if y<0 then
begin
outstring(1,‘divmod: отрицательный делитель: ’); outinteger(1,y);
outstring(1,‘‘NL’’);
q := 0;
r := 0
end
else
begin
q := x÷y;
r := x - y×q;
if x<0 ∧ r≠0 then
begin
q := q - 1;
r := r + y
end
end
end divmod;
comment Вычисление юлианской даты по календарной. Глобальная переменная
JulFlag определяет используемый календарь: JulFlag<0 - юлианский,
JulFlag>0 - григорианский, при JulFlag=0 с 15.10.1582 года
используется григорианский календарь, иначе - юлианский.
Функция дает корректный результат для любых дат;
real procedure julday(j,m,d,h,min,s);
value j,m,d,h,min,s; integer j,m,d,h,min; real s;
begin
integer c0,j0,j1,j2,x1,x2,x3,x4,jd;
if JulFlag<0 ∨ JulFlag=0 ∧ (j<1582 ∨ (j=1582 ∧ m<10) ∨ (j=1582 ∧ m=10 ∧ d<15)) then
begin
comment юлианский календарь;
j0 := 1721117;
c0 := div(m-3,12);
j1 := div(1461×(j+c0),4);
j2 := div(153×m-1836×c0-457,5);
jd := j1 + j2 + d + j0
end
else
begin
comment григорианский календарь;
c0 := div(m-3,12);
x4 := j + c0;
divmod(x4,100,x3,x2);
x1 := m - 12×c0 - 3;
jd := div(146097×x3,4) + div(36525×x2,100) + div(153×x1+2,5) + d + 1721119
end;
julday := jd + (h+(min+s/60)/60)/24 - 0.5
end julday;
comment Вычисление календарной даты по юлианской. Глобальная переменная
JulFlag определяет используемый календарь: JulFlag<0 - юлианский,
JulFlag>0 - григорианский, при JulFlag=0 и jd>2299160 используется
григорианский календарь, иначе - юлианский.
Процедура дает корректный результат для любых юлианских дат,
включая отрицательные;
procedure caldat(jd,y,m,d,h,min,s);
value jd; real jd,s; integer y,m,d,h,min;
begin
real dt; integer j,c0,k1,k2,x1,x2,x3,r1,r2,r3,y2;
comment выделить дробную часть суток;
dt := jd + 0.5;
j := jd;
dt := (dt-j)×24.0;
if dt<0 then
begin
j := j - 1; dt := dt + 24.0
end;
h := entier(dt);
dt:= (dt-h)×60.0;
min := entier(dt);
s := (dt-min)×60.0;
if JulFlag<0 ∨ JulFlag=0 ∧ j ≤ 2299160 then
begin
comment юлианский календарь;
y2 := j - 1721118;
k2 := 4×y2 + 3;
k1 := 5×div(mod(k2,1461),4) + 2;
x1 := div(k1,153);
c0 := div(x1+2,12);
y := div(k2,1461) + c0;
m := x1 - 12×c0 + 3;
d := div(mod(k1,153),5) + 1
end
else
begin
comment григорианский календарь;
divmod(4×j-6884477,146097,x3,r3);
divmod(100×div(r3,4)+99,36525,x2,r2);
divmod(5×div(r2,100)+2,153,x1,r1);
d := div(r1,5) + 1;
c0 := div(x1+2,12);
y := 100×x3 + x2 + c0;
m := x1 - 12×c0 + 3
end
end caldat;
procedure test(y,m,d,h,min,s);
value y,m,d,h,min,s; integer y,m,d,h,min; real s;
begin
real jd;
jd := julday(y,m,d,h,min,s);
outreal(1,jd);
caldat(jd,y,m,d,h,min,s);
outinteger(1,y); outinteger(1,m); outinteger(1,d);
outinteger(1,h); outinteger(1,min); outreal(1,s);
outstring(1,‘‘NL’’)
end test;
JulFlag := 0;
test(1900,1,0,12,0,0);
test(2000,1,1,12,0,0);
test(-5000,7,12,16,6,32.87);
test(5000,7,12,16,6,32.87)
end
Tags: