天泣記

time.c の、find_time_t をひさしぶりにいじる。

これが何をやる関数かというと、年月日時分秒の値から time_t の値を得る関数である。

これは、夏時間と閏秒を考慮すると、そう簡単ではない。

たとえば、PST8PDT という、アメリカ西海岸のタイムゾーンの去年から今年にかけての夏時間の切り替わりと閏秒のところを図にすると以下のようになる。ちなみに、PST と PDT をあわせたものが PST8PDT である。(ついでに GMT と JST もちょっと描いてある)

ここで PST/PDT の関係にしたがって縦軸の値から横軸の値を求めるというのが find_time_t の機能である。(横軸は UTC な時刻で目盛をふってあるが、ここでは time_t だと考える。)

これをどうやるかというと、mktime をつかうというのも一つの方法であるが、それをやると (閏秒込みの) UTC/GMT に対応できないなどの問題があるので、localtime/gmtime を使う。

localtime/gmtime は図の横軸から縦軸の値を求める関数である。これの逆関数を作ればいい。

具体的にどうやるかと言えば、localtime/gmtime にいろんな値を食わせてみて、縦軸が与えられた値になるような横軸の値を探すわけである。

幸いにして localtime/gmtime は「だいたい」単調なので、二分探索でもなんとか動かせる。(秋に単調でない部分があるが、考えるとその場合でも二分探索でどちらかひとつは見つかることがわかる。)

さて、以前、最初に find_time_t を実装したときは二分探索を使っていたのだが、すぐに補間探索 (interpolation search) を実装して、長年それでやってきた。そうやったのは2001年だから、8年くらいになる。

それでとくに文句も出ていないのだが、補間探索というのはじつはよくないとは思っていた。

というのは、補間探索はうまくいくときには O(log log n) で速いのだが、うまくいかないときには O(n) になってしまうのである。そして、デバッグでステップ実行したときに、どうも長くかかる、という経験もあった。

で、思い立って測ってみると、春の切り替わりの周辺で 90回以上 localtime を呼び出しているケースが存在することがわかった。この環境の time_t は 32bit であり、素朴に二分探索しても 32回で終わるはずで、多すぎである。

で、考えた結果、この関数はほぼ直線である (そして傾きが一定である) ことを積極的に利用することにした。

二分探索と同じく、はさみうちにするのは変わらないが、端の値から、関数が直線であることを仮定して答えを推測し、その値で探索を進める。

端には下端と上端があるので両方でやる。下端から上端まで夏時間の切り替わりも閏秒もなければ一発で見つかるし、ひとつしかなければどちらかで見つかる。

あと、補間探索と同じくひどい目にあうかもしれないので、両端から試したら次は半分に分割して探索する。つまり、上端から推測し、下端から推測し、半分と推測する、というのを繰り返す。3回に 1回は探索範囲が半分になるので、最悪でも O(log n) で済む。

そうやって実装してみると、ずいぶんとマシになって、90回以上かかっていたものも 20回以下で済むようになった。

なお、そうやってもやはり春がいちばん手間がかかる。というのは、春は縦軸の時刻に対して横軸の時刻が存在しないことがあって、そうすると推測は確実に外れて、二分探索になってしまうからである。

time-system.asy:

import graph;

real gmt(real t) {
  // 2008-11-01 00:00:00 UTC -> 2008-11-01 00:00:00 GMT
  return t;
}

real jp(real t) {
  // 2008-11-01 00:00:00 UTC -> 2008-11-01 09:00:00 +0900 (JST)
  return t + 9*3600;
}

real nov1 = 0*24*3600;
real nov2 = 1*24*3600;
real nov2_1 = 1*24*3600+1*3600;
real nov2_2 = 1*24*3600+2*3600;
real nov2_3 = 1*24*3600+3*3600;
real nov2_4 = 1*24*3600+4*3600;
real nov2_5 = 1*24*3600+5*3600;
real nov2_6 = 1*24*3600+6*3600;
real nov2_7 = 1*24*3600+7*3600;
real nov2_8 = 1*24*3600+8*3600;
real nov2_9 = 1*24*3600+9*3600;
real nov2_10 = 1*24*3600+10*3600;
real nov2_11 = 1*24*3600+11*3600;
real nov2_12 = 1*24*3600+12*3600;
real nov2_20 = 1*24*3600+20*3600;
real nov3 = 2*24*3600;

real dec31 = 0*3600;
real dec31_12 = 12*3600;
real dec31_13 = 13*3600;
real dec31_14 = 14*3600;
real dec31_15 = 15*3600;
real dec31_16 = 16*3600;
real dec31_17 = 17*3600;
real dec31_18 = 18*3600;
real dec31_19 = 19*3600;
real dec31_20 = 20*3600;
real dec31_21 = 21*3600;
real dec31_22 = 22*3600;
real dec31_23 = 23*3600;
real dec31_23_59 = 23*3600 + 59*60;
real dec31_23_59_50 = 23*3600 + 59*60 + 50;
real dec31_23_59_55 = 23*3600 + 59*60 + 55;
real dec31_23_59_56 = 23*3600 + 59*60 + 56;
real dec31_23_59_57 = 23*3600 + 59*60 + 57;
real dec31_23_59_58 = 23*3600 + 59*60 + 58;
real dec31_23_59_59 = 23*3600 + 59*60 + 59;
real dec31_23_59_60 = 23*3600 + 59*60 + 60;
real dec31_24 = 24*3600;
real jan1 = dec31_24 + 1;
real jan1_0_0_1 = jan1+1;
real jan1_0_0_2 = jan1+2;
real jan1_0_0_3 = jan1+3;
real jan1_0_0_4 = jan1+4;
real jan1_0_0_5 = jan1+5;
real jan1_0_0_10 = jan1+10;
real jan1_0_1 = jan1+60;
real jan1_1 = jan1+1*3600;
real jan1_2 = jan1+2*3600;
real jan1_3 = jan1+3*3600;
real jan1_4 = jan1+4*3600;
real jan1_5 = jan1+5*3600;
real jan1_6 = jan1+6*3600;
real jan1_7 = jan1+7*3600;
real jan1_8 = jan1+8*3600;
real jan1_9 = jan1+9*3600;
real jan1_10 = jan1+10*3600;
real jan1_11 = jan1+11*3600;
real jan1_12 = jan1+12*3600;

real mar8 = 0*24*3600;
real mar8_1 = 1*3600;
real mar8_2 = 2*3600;
real mar8_3 = 3*3600;
real mar8_4 = 4*3600;
real mar8_5 = 5*3600;
real mar8_6 = 6*3600;
real mar8_7 = 7*3600;
real mar8_8 = 8*3600;
real mar8_9 = 9*3600;
real mar8_10 = 10*3600;
real mar8_11 = 11*3600;
real mar8_12 = 12*3600;
real mar8_13 = 13*3600;
real mar8_14 = 14*3600;
real mar8_15 = 15*3600;
real mar8_20 = 20*3600;
real mar9 = 1*24*3600;

real nov_us_pacific(real t) {
  // 2008-11-01 00:00:00 UTC -> 2008-10-31 17:00:00 -0700 (PDT)
  // 2008-11-02 08:59:59 UTC -> 2008-11-02 01:59:59 -0800 (PDT)
  // 2008-11-02 09:00:00 UTC -> 2008-11-02 01:00:00 -0800 (PST)
  if (t < (24+9)*3600)
    return t - 7*3600;
  else
    return t - 8*3600;
}

real dec_us_pacific(real t) {
  // right/PST8PDT  Wed Dec 31 23:59:60 2008 UTC = Wed Dec 31 15:59:60 2008 PST isdst=0 gmtoff=-28800
  // right/PST8PDT  Thu Jan  1 00:00:00 2009 UTC = Wed Dec 31 16:00:00 2008 PST isdst=0 gmtoff=-28800
  return t - 8*3600;
}

real mar_us_pacific(real t) {
  // PST8PDT  Sun Mar  8 09:59:59 2009 UTC = Sun Mar  8 01:59:59 2009 PST isdst=0 gmtoff=-28800
  // PST8PDT  Sun Mar  8 10:00:00 2009 UTC = Sun Mar  8 03:00:00 2009 PDT isdst=1 gmtoff=-25200
  //
  // 2009-03-08 00:00:00 UTC -> 2008-03-07 16:00:00 -0800 (PST)
  // 2009-03-08 09:59:59 UTC -> 2009-03-08 01:59:59 -0800 (PST)
  // 2009-03-08 10:00:00 UTC -> 2009-03-08 03:00:00 -0700 (PDT)
  if (t < 10*3600)
    return t - 8*3600;
  else
    return t - 7*3600;
}

string nov_datefmt(real t) {
  int x = floor(t);
  int s = x % 60;
  int m = quotient(x, 60) % 60;
  int h = quotient(x, 3600) % 24;
  int d = quotient(x, 86400)+1;
  return format("2008-11-%02d", d) +
         format(" %02d", h) + 
         format(":%02d", m) + 
         format(":%02d", s);
}

string nov_xticklabel(real x) {
  if (nov1 <= x)
    return nov_datefmt(x) + " UTC";
  else
    return format("$%g$", x);
}

string nov_yticklabel(real y) {
  if (nov1 <= y)
    return nov_datefmt(y);
  else
    return format("$%g$", y);
}

string dec_utc_datefmt(real t) {
  int x = floor(t);
  bool before_leap_second = x < jan1;
  if (x == dec31_24)
    return "2008-12-31 23:59:60";
  else if (!before_leap_second)
    x = x-1;
  int s = x % 60;
  int m = quotient(x, 60) % 60;
  int h = quotient(x, 3600) % 24;
  int day = before_leap_second ? 31 : 1;
  int mon = before_leap_second ? 12 : 1;
  int year = before_leap_second ? 2008 : 2009;
  return format("%d", year) +
         format("-%02d", mon) + 
         format("-%02d", day) + 
         format(" %02d", h) + 
         format(":%02d", m) + 
         format(":%02d", s);
}

string dec_us_pacific_datefmt(real t) {
  int x = floor(t);
  bool before_leap_second = x < dec_us_pacific(jan1);
  if (x == dec_us_pacific(dec31_24))
    return "2008-12-31 15:59:60";
  else if (!before_leap_second)
    x = x-1;
  int s = x % 60;
  int m = quotient(x, 60) % 60;
  int h = quotient(x, 3600) % 24;
  int day = 31;
  int mon = 12;
  int year = 2008;
  return format("%d", year) +
         format("-%02d", mon) + 
         format("-%02d", day) + 
         format(" %02d", h) + 
         format(":%02d", m) + 
         format(":%02d", s);
}

string dec_xticklabel(real x) {
  return dec_utc_datefmt(x) + " UTC";
}

string dec_utc_yticklabel(real y) {
  return dec_utc_datefmt(y);
}

string dec_us_pacific_yticklabel(real y) {
  return dec_us_pacific_datefmt(y);
}

string mar_datefmt(real t) {
  int x = floor(t);
  int s = x % 60;
  int m = quotient(x, 60) % 60;
  int h = quotient(x, 3600) % 24;
  int d = quotient(x, 86400)+8;
  return format("2009-03-%02d", d) +
         format(" %02d", h) + 
         format(":%02d", m) + 
         format(":%02d", s);
}

string mar_xticklabel(real x) {
  return mar_datefmt(x) + " UTC";
}

string mar_yticklabel(real y) {
  return mar_datefmt(y);
}

Label xlabel = shift(13mm,-3mm)*rotate(-20)*Label("$%g$");
Label ylabel = shift(0mm,0mm)*rotate(0)*Label("$%g$");

label("Autumn DST Switch", (50,170));
label("Leap Second", (350,170));
label("Spring DST Switch", (650,170));

draw((30,10)--(0,-115), palecyan);
draw((49,10)--(110,-115), palecyan);

draw((351,17)--(300,-90), palered);
draw((353,17)--(410,-90), palered);

draw((634,20)--(600,-60), palecyan);
draw((653,20)--(710,-60), palecyan);

picture nov_pic1;
size(nov_pic1, 10cm);
filldraw(nov_pic1, shift(nov2_7, nov2)*scale(3600*4,3600*3)*unitsquare, palecyan, black);
draw(nov_pic1, graph(gmt, nov2, nov3));
draw(nov_pic1, graph(nov_us_pacific, nov2, nov2+9*3600-1));
draw(nov_pic1, graph(nov_us_pacific, nov2+9*3600, nov3));
draw(nov_pic1, graph(jp, nov2, nov3));
real[] xmajorticks = { nov2, nov2_12, nov3 };
real[] xminorticks;
for (int i = 1; i <= 23; ++i) xminorticks[xminorticks.length] = nov2+i*3600;
ticks xticks = Ticks(xlabel, nov_xticklabel, Ticks=xmajorticks, ticks=xminorticks);
xaxis(nov_pic1, YEquals(nov2, false), xticks);
real[] ymajorticks = { nov2, nov2_12, nov3 };
real[] yminorticks;
for (int i = -11; i <= 35; ++i) yminorticks[yminorticks.length] = nov2+i*3600;
ticks yticks = Ticks(ylabel, nov_yticklabel, Ticks=ymajorticks, ticks=yminorticks);
yaxis(nov_pic1, XEquals(nov2), ymin=nov2-11*3600, yticks);
label(nov_pic1, "JST", (nov2_20, jp(nov2_20)), E);
label(nov_pic1, "GMT", (nov2_20, gmt(nov2_20)), E);
label(nov_pic1, "PST", (nov2_20, nov_us_pacific(nov2_20)), E);
label(nov_pic1, "PDT", (nov2_5, nov_us_pacific(nov2_5)));
nov_pic1 = shift(-nov2,-nov2)*nov_pic1;
add(nov_pic1.scale(300,300),(0,0));

picture nov_pic2;
size(nov_pic2, 10cm);
draw(nov_pic2, graph(nov_us_pacific, nov2_7, nov2_9-1));
draw(nov_pic2, graph(nov_us_pacific, nov2_9, nov2_11));
fill(nov_pic2, shift(nov2_9, nov2_1)*scale(200)*unitcircle);
filldraw(nov_pic2, shift(nov2_9, nov2_2)*scale(200)*unitcircle, white, black);
real[] xmajorticks = { nov2_6, nov2_7, nov2_8, nov2_9, nov2_10, nov2_11 };
ticks xticks = Ticks(xlabel, nov_xticklabel, Ticks=xmajorticks);
xaxis(nov_pic2, YEquals(nov2, false), xticks);
real[] ymajorticks = { nov2, nov2_1, nov2_2, nov2_3 };
ticks yticks = Ticks(ylabel, nov_yticklabel, Ticks=ymajorticks);
yaxis(nov_pic2, XEquals(nov2_7), ymin=nov2, yticks);
label(nov_pic2, "PDT", (nov2_8+1800, nov_us_pacific(nov2_8+1800)), NW);
label(nov_pic2, "PST", (nov2_10, nov_us_pacific(nov2_10)), SE);
nov_pic2 = shift(-nov2_7,-nov_us_pacific(nov2_7))*nov_pic2;
add(nov_pic2.scale(300,300),(0,-200));

picture dec_pic1;
size(dec_pic1, 10cm);
draw(dec_pic1, graph(gmt, dec31_12, jan1_12));
draw(dec_pic1, graph(jp, dec31_12, jan1_12));
draw(dec_pic1, graph(dec_us_pacific, dec31_12, jan1_12));
real[] xmajorticks = { dec31_12, jan1, jan1_12};
real[] xminorticks;
for (int i = 1; i <= 23; ++i) xminorticks[xminorticks.length] = dec31_12+i*3600;
ticks xticks = Ticks(xlabel, dec_xticklabel, Ticks=xmajorticks, ticks=xminorticks);
xaxis(dec_pic1, YEquals(dec31_12, false), xticks);
real[] ymajorticks = { dec31_12, jan1, jan1_12};
real[] yminorticks;
for (int i = -11; i <= 35; ++i) yminorticks[yminorticks.length] = dec31_12+i*3600;
ticks yticks = Ticks(ylabel, dec_utc_yticklabel, Ticks=ymajorticks, ticks=yminorticks);
yaxis(dec_pic1, XEquals(dec31_12), ymin=dec31_12-11*3600, yticks);
label(dec_pic1, "JST", (jan1_9, jp(jan1_9)), E);
label(dec_pic1, "GMT", (jan1_9, gmt(jan1_9)), E);
label(dec_pic1, "PST", (jan1_9, dec_us_pacific(jan1_9)), E);
draw(dec_pic1, (dec31_24,dec31_12-11*3600)--(dec31_24,jan1+21*3600), palered);
draw(dec_pic1, (dec31_12,gmt(dec31_24))--(jan1_12,gmt(dec31_24)), palered);
draw(dec_pic1, (dec31_12,jp(dec31_24))--(jan1_12,jp(dec31_24)), palered);
draw(dec_pic1, (dec31_12,dec_us_pacific(dec31_24))--(jan1_12,dec_us_pacific(dec31_24)), palered);
dot(dec_pic1, (dec31_24, jp(dec31_24)), red);
dot(dec_pic1, (dec31_24, gmt(dec31_24)), red);
dot(dec_pic1, (dec31_24, dec_us_pacific(dec31_24)), red);
dec_pic1 = shift(-dec31_12,-dec31_12)*dec_pic1;
add(dec_pic1.scale(300,300),(300,0));

picture dec_pic2;
size(dec_pic2, 10cm);
fill(dec_pic2, shift(dec31_23_59_60, dec_us_pacific(dec31_23_59_58))*scale(1, jan1_0_0_2-dec31_23_59_58)*unitsquare, palered);
fill(dec_pic2, shift(dec31_23_59_58, dec_us_pacific(dec31_23_59_60))*scale(jan1_0_0_2-dec31_23_59_58, 1)*unitsquare, palered);
draw(dec_pic2, graph(dec_us_pacific, dec31_23_59_58, jan1_0_0_2));
draw(dec_pic2, graph(dec_us_pacific, dec31_23_59_60, jan1), red+2);
real[] xmajorticks = { dec31_23_59_58, dec31_23_59_60, jan1, jan1_0_0_2 };
real[] xminorticks;
for (int i = -2; i <= 2; ++i) xminorticks[xminorticks.length] = dec31_23_59_60+i;
ticks xticks = Ticks(xlabel, dec_xticklabel, Ticks=xmajorticks, ticks=xminorticks);
xaxis(dec_pic2, YEquals(dec_us_pacific(dec31_23_59_58), false), xticks);
real[] ymajorticks = { dec_us_pacific(dec31_23_59_58), dec_us_pacific(dec31_23_59_60), dec_us_pacific(jan1), dec_us_pacific(jan1_0_0_2) };
real[] yminorticks;
for (int i = -2; i <= 2; ++i) yminorticks[yminorticks.length] = dec_us_pacific(dec31_23_59_60)+i;
ticks yticks = Ticks(ylabel, dec_us_pacific_yticklabel, Ticks=ymajorticks, ticks=yminorticks);
yaxis(dec_pic2, XEquals(dec31_23_59_58), ymin=dec_us_pacific(dec31_23_59_58), yticks);
label(dec_pic2, "leap second", (dec31_23_59_60+0.5, dec_us_pacific(dec31_23_59_60+0.5)), E);
label(dec_pic2, "PST", (jan1_0_0_1, dec_us_pacific(jan1_0_0_1)), E);
dec_pic2 = shift(-dec31_23_59_58,-dec_us_pacific(dec31_23_59_58))*dec_pic2;
add(dec_pic2.scale(300,300),(300,-200));

picture mar_pic1;
size(mar_pic1, 10cm);
filldraw(mar_pic1, shift(mar8_8, mar8)*scale(3600*4,3600*5)*unitsquare, palecyan, black);
draw(mar_pic1, graph(gmt, mar8, mar9));
draw(mar_pic1, graph(jp, mar8, mar9));
draw(mar_pic1, graph(mar_us_pacific, mar8, mar8_10-1));
draw(mar_pic1, graph(mar_us_pacific, mar8_10, mar9));
real[] xmajorticks = { mar8, mar8_12, mar9 };
real[] xminorticks;
for (int i = 1; i <= 23; ++i) xminorticks[xminorticks.length] = mar8+i*3600;
ticks xticks = Ticks(xlabel, mar_xticklabel, Ticks=xmajorticks, ticks=xminorticks);
xaxis(mar_pic1, YEquals(mar8, false), xticks);
real[] ymajorticks = { mar8, mar8_12, mar9 };
real[] yminorticks;
for (int i = -11; i <= 35; ++i) yminorticks[yminorticks.length] = mar8+i*3600;
ticks yticks = Ticks(ylabel, mar_yticklabel, Ticks=ymajorticks, ticks=yminorticks);
yaxis(mar_pic1, XEquals(mar8), ymin=mar8-11*3600, yticks);
label(mar_pic1, "JST", (mar8_20, jp(mar8_20)), E);
label(mar_pic1, "GMT", (mar8_20, gmt(mar8_20)), E);
label(mar_pic1, "PDT", (mar8_20, mar_us_pacific(mar8_20)), E);
label(mar_pic1, "PST", (mar8_6, mar_us_pacific(mar8_6)));
mar_pic1 = shift(-mar8,-mar8)*mar_pic1;
add(mar_pic1.scale(300,300),(600,0));

picture mar_pic2;
size(mar_pic2, 10cm);
draw(mar_pic2, graph(mar_us_pacific, mar8_8, mar8_10-1));
draw(mar_pic2, graph(mar_us_pacific, mar8_10, mar8_12));
fill(mar_pic2, shift(mar8_10, mar8_3)*scale(200)*unitcircle);
filldraw(mar_pic2, shift(mar8_10, mar8_2)*scale(200)*unitcircle, white, black);
real[] xmajorticks = { mar8_8, mar8_9, mar8_10, mar8_11, mar8_12 };
ticks xticks = Ticks(xlabel, mar_xticklabel, Ticks=xmajorticks);
xaxis(mar_pic2, YEquals(mar8, false), xticks);
real[] ymajorticks = { mar8, mar8_1, mar8_2, mar8_3, mar8_4, mar8_5 };
ticks yticks = Ticks(ylabel, mar_yticklabel, Ticks=ymajorticks);
yaxis(mar_pic2, XEquals(mar8_8), ymin=mar8, yticks);
label(mar_pic2, "PST", (mar8_9+1800, mar_us_pacific(mar8_9+1800)), NW);
label(mar_pic2, "PDT", (mar8_11, mar_us_pacific(mar8_11)), SE);
mar_pic2 = shift(-mar8_8,-mar_us_pacific(mar8_8))*mar_pic2;
add(mar_pic2.scale(300,300),(600,-200));

Asymptote で複数のグラフを一枚に描いたのだが、ふたつのグラフの関係を書き込むのが難しい。

個々のグラフの座標系で点を定義して、その間に線を引きたいのだが、やりかたがよくわからない。

グローバルな座標系で指定すればもちろん引けるが、指定するのが大変である。

整理して描き直すとこうか。

縦軸を毎分61秒のスケールにした結果、グラフ全部で epoch を共通化でき、時刻のフォーマットもまとめられた。

time-system5.asy:

import graph;

int min = 60;
int hour = 3600;
int day = 24hour;

real feb1 = 31day;
real mar1 = feb1 + 29day;
real apr1 = mar1 + 31day;
real may1 = apr1 + 30day;
real jun1 = may1 + 31day;
real jul1 = jun1 + 30day;
real aug1 = jul1 + 31day;
real sep1 = aug1 + 31day;
real oct1 = sep1 + 30day;
real nov1 = oct1 + 31day;
real dec1 = nov1 + 30day;
real dec31_23_59_60 = dec1 + 31day; // the leap second
real jan1 = dec31_23_59_60 + 1;

real mar8 = mar1 + 7day;
real mar8_12 = mar8 + 12hour;
real mar8_13 = mar8 + 13hour;
real mar9 = mar1 + 8day;
real mar9_1 = mar9 + 1hour;
real mar9_2 = mar9 + 2hour;
real mar9_3 = mar9 + 3hour;
real mar9_4 = mar9 + 4hour;
real mar9_5 = mar9 + 5hour;
real mar9_6 = mar9 + 6hour;
real mar9_7 = mar9 + 7hour;
real mar9_8 = mar9 + 8hour;
real mar9_9 = mar9 + 9hour;
real mar9_10 = mar9 + 10hour;
real mar9_11 = mar9 + 11hour;
real mar9_12 = mar9 + 12hour;
real mar9_20 = mar9 + 20hour;
real mar10 = mar1 + 9day;
real mar10_12 = mar10 + 12hour;

real nov1_12 = nov1 + 12hour;
real nov1_13 = nov1 + 13hour;
real nov2 = nov1 + 1day;
real nov2_1 = nov2 + 1hour;
real nov2_2 = nov2 + 2hour;
real nov2_3 = nov2 + 3hour;
real nov2_4 = nov2 + 4hour;
real nov2_5 = nov2 + 5hour;
real nov2_6 = nov2 + 6hour;
real nov2_7 = nov2 + 7hour;
real nov2_8 = nov2 + 8hour;
real nov2_9 = nov2 + 9hour;
real nov2_10 = nov2 + 10hour;
real nov2_11 = nov2 + 11hour;
real nov2_12 = nov2 + 12hour;
real nov2_20 = nov2 + 20hour;
real nov3 = nov1 + 2day;
real nov3_12 = nov3 + 12hour;

real dec31 = dec1 + 30day;
real dec31_1 = dec31 + 1hour;
real dec31_12 = dec31 + 12hour;
real dec31_15 = dec31 + 15hour;
real dec31_15_59 = dec31_15 + 59min;
real dec31_15_59_58 = dec31_15_59 + 58;
real dec31_15_59_59 = dec31_15_59 + 59;
real dec31_16 = dec31 + 16hour;
real dec31_16_0_1 = dec31_16 + 1;
real dec31_16_0_2 = dec31_16 + 2;
real dec31_23_59_58 = dec31 + 23hour + 59min + 58;
real dec31_23_59_59 = dec31 + 23hour + 59min + 59;
real jan1_0_0_1 = jan1 + 1;
real jan1_0_0_2 = jan1 + 2;
real jan1_8 = jan1 + 8hour;
real jan1_12 = jan1 + 12hour;
real jan2 = jan1 + 1day;

int lmin = 61;
int lhour = 60*lmin;
int lday = 24*lhour;

string x_datefmt(int t) {
  if (t < dec31_23_59_60) {
    int s = t % min;
    int m = quotient(t, min) % 60;
    int h = quotient(t, hour) % 24;
    int d = quotient(t, day)+1;
    string hms = format(" %02d", h) + format(":%02d", m) + format(":%02d", s);
    if (d <= 31)
      return format("2008-01-%02d", d) + hms;
    d -= 31;
    if (d <= 29)
      return format("2008-02-%02d", d) + hms;
    d -= 29;
    if (d <= 31)
      return format("2008-03-%02d", d) + hms;
    d -= 31;
    if (d <= 30)
      return format("2008-04-%02d", d) + hms;
    d -= 30;
    if (d <= 31)
      return format("2008-05-%02d", d) + hms;
    d -= 31;
    if (d <= 30)
      return format("2008-06-%02d", d) + hms;
    d -= 30;
    if (d <= 31)
      return format("2008-07-%02d", d) + hms;
    d -= 31;
    if (d <= 31)
      return format("2008-08-%02d", d) + hms;
    d -= 31;
    if (d <= 30)
      return format("2008-09-%02d", d) + hms;
    d -= 30;
    if (d <= 31)
      return format("2008-10-%02d", d) + hms;
    d -= 31;
    if (d <= 30)
      return format("2008-11-%02d", d) + hms;
    d -= 30;
    return format("2008-12-%02d", d) + hms;
  }
  else if (t < jan1) {
    return "2008-12-31 23:59:60";
  }
  else {
    --t;
    int s = t % 60;
    int m = quotient(t, 60) % 60;
    int h = quotient(t, 3600) % 24;
    int d = quotient(t, 86400)+1-366;
    string hms = format(" %02d", h) + format(":%02d", m) + format(":%02d", s);
    return format("2009-01-%02d", d) + hms;
  }
}

string xticklabel(real x) {
  int t = floor(x);
  return x_datefmt(t) + " UTC";
}

string y_datefmt(int t) {
  int s = t % lmin;
  int m = quotient(t, lmin) % 60;
  int h = quotient(t, lhour) % 24;
  int d = quotient(t, lday)+1;
  string hms = format(" %02d", h) + format(":%02d", m) + format(":%02d", s);
  if (d <= 31)
    return format("2008-01-%02d", d) + hms;
  d -= 31;
  if (d <= 29)
    return format("2008-02-%02d", d) + hms;
  d -= 29;
  if (d <= 31)
    return format("2008-03-%02d", d) + hms;
  d -= 31;
  if (d <= 30)
    return format("2008-04-%02d", d) + hms;
  d -= 30;
  if (d <= 31)
    return format("2008-05-%02d", d) + hms;
  d -= 31;
  if (d <= 30)
    return format("2008-06-%02d", d) + hms;
  d -= 30;
  if (d <= 31)
    return format("2008-07-%02d", d) + hms;
  d -= 31;
  if (d <= 31)
    return format("2008-08-%02d", d) + hms;
  d -= 31;
  if (d <= 30)
    return format("2008-09-%02d", d) + hms;
  d -= 30;
  if (d <= 31)
    return format("2008-10-%02d", d) + hms;
  d -= 31;
  if (d <= 30)
    return format("2008-11-%02d", d) + hms;
  d -= 30;
  if (d <= 31)
    return format("2008-12-%02d", d) + hms;
  d -= 31;
  return format("2009-01-%02d", d) + hms;
}

string yticklabel(real y) {
  int t = floor(y);
  return y_datefmt(t);
}

real gmt(real x) {
  // 2008-01-01 00:00:00 UTC -> 2008-01-01 00:00:00 GMT
  int t = floor(x);
  real subsec = x-t;
  if (t < dec31_23_59_60)
    return t + quotient(t, 60) + subsec;
  if (t < jan1)
    return dec31_23_59_60 + quotient(floor(dec31_23_59_60), 60) + subsec;
  return (t-1) + quotient(t-1, 60) + subsec;
}

real jp(real x) {
  // 2008-11-01 00:00:00 UTC -> 2008-11-01 09:00:00 +0900 (JST)
  int t = floor(x);
  real subsec = x-t;
  if (t < dec31_23_59_60)
    return t + quotient(t, 60) + subsec + 9lhour;
  if (t < jan1)
    return dec31_23_59_60 + quotient(floor(dec31_23_59_60), 60) + subsec + 9lhour;
  return (t-1) + quotient(t-1, 60) + subsec + 9lhour;
}

real us_pacific(real x) {
  // 2008-11-01 00:00:00 UTC -> 2008-11-01 09:00:00 +0900 (JST)
  int t = floor(x);
  real subsec = x-t;
  if (t < mar9_10)
    return t + quotient(t, 60) + subsec - 8lhour;
  if (t < nov2_9)
    return t + quotient(t, 60) + subsec - 7lhour;
  if (t < dec31_23_59_60)
    return t + quotient(t, 60) + subsec - 8lhour;
  if (t < jan1)
    return dec31_23_59_60 + quotient(floor(dec31_23_59_60), 60) - 1 + subsec - 8lhour;
  return (t-1) + quotient(t-1, 60) + subsec - 8lhour;
}

real reg1_beg = mar9-1min;
real reg1_end = mar10+1min;
real reg2_beg = nov2-1min;
real reg2_end = nov3+1min;
real reg3_beg = dec31_12-1min;
real reg3_end = jan1_12+1min;

real aspect = 60/61;
real usize1x = 0.0010;
real usize2x = 0.007;
real usize3x = 22;

real usize1y = usize1x*aspect;
pair usize1 = (usize1x,usize1y);
real usize2y = usize2x*aspect;
pair usize2 = (usize2x,usize2y);
real usize3y = usize3x*aspect;
pair usize3 = (usize3x,usize3y);

pair step = (180,180);

pair mar_pos1 = (0,0);
pair nov_pos1 = mar_pos1 + step;
pair dec_pos1 = mar_pos1 + 2step;

pair mar_pos2 = (-170,240);
pair nov_pos2 = mar_pos2 + step - (usize2.x*hour,0);
pair dec_pos2 = mar_pos2 + 2step;
pair dec_pos3 = (mar_pos2.x,dec_pos2.y);

pair mar_orig1 = (mar9, gmt(mar8_13));
pair nov_orig1 = (nov2, gmt(nov1_13));
pair dec_orig1 = (dec31_12, gmt(dec31_1));

pair mar_orig2 = (mar9_8, gmt(mar9));
pair nov_orig2 = (nov2_7, gmt(nov2));
pair dec_orig2 = (dec31_23_59_58, gmt(dec31_15_59_58));
pair dec_orig3 = (dec31_15_59_58, gmt(dec31_15_59_58));

transform scalep(pair p) {
  return scale(p.x, p.y);
}

transform mar_trans1 = shift(mar_pos1)*scalep(usize1)*shift(-mar_orig1);
transform mar_trans2 = shift(mar_pos2)*scalep(usize2)*shift(-mar_orig2);

transform nov_trans1 = shift(nov_pos1)*scalep(usize1)*shift(-nov_orig1);
transform nov_trans2 = shift(nov_pos2)*scalep(usize2)*shift(-nov_orig2);

transform dec_trans1 = shift(dec_pos1)*scalep(usize1)*shift(-dec_orig1);
transform dec_trans2 = shift(dec_pos2)*scalep(usize3)*shift(-dec_orig2);
transform dec_trans3 = shift(dec_pos3)*scalep(usize3)*shift(-dec_orig3);

path rect(pair min, pair max) {
  return min--(max.x,min.y)--max--(min.x,max.y)--cycle;
}
fill(rect(dec_trans3*(dec31_15_59_58,gmt(dec31_16)-1), dec_trans2*(jan1_0_0_2,gmt(dec31_16))),lightred);

draw(mar_trans1*(mar9_8,gmt(mar9))--mar_trans2*(mar9_8,gmt(mar9)), lightblue);
draw(mar_trans1*(mar9_12,gmt(mar9_5))--mar_trans2*(mar9_12,gmt(mar9_5)), lightblue);

draw(nov_trans1*(nov2_7,gmt(nov2))--nov_trans2*(nov2_7,gmt(nov2)), lightblue);
draw(nov_trans1*(nov2_11,gmt(nov2_3))--nov_trans2*(nov2_11,gmt(nov2_3)), lightblue);

draw(dec_trans1*(dec31_23_59_58,gmt(dec31_15_59_58))--dec_trans2*(dec31_23_59_58,gmt(dec31_15_59_58)), lightblue);
draw(dec_trans1*(jan1_0_0_2,gmt(dec31_16_0_2))--dec_trans2*(jan1_0_0_2,gmt(dec31_16_0_2)), lightblue);

draw(dec_trans1*(dec31_15_59_58,gmt(dec31_15_59_58))--dec_trans3*(dec31_15_59_58,gmt(dec31_15_59_58)), lightblue);
draw(dec_trans1*(dec31_16_0_2,gmt(dec31_16_0_2))--dec_trans3*(dec31_16_0_2,gmt(dec31_16_0_2)), lightblue);

draw(mar_trans1*(mar10,gmt(mar10))--nov_trans1*(nov2,gmt(nov2)), Dotted);
draw(mar_trans1*(mar10,jp(mar10))--nov_trans1*(nov2,jp(nov2)), Dotted);
draw(mar_trans1*(mar10,us_pacific(mar10))--nov_trans1*(nov2,us_pacific(nov2)), Dotted);

draw(nov_trans1*(nov3,gmt(nov3))--dec_trans1*(dec31_12,gmt(dec31_12)), Dotted);
draw(nov_trans1*(nov3,jp(nov3))--dec_trans1*(dec31_12,jp(dec31_12)), Dotted);
draw(nov_trans1*(nov3,us_pacific(nov3))--dec_trans1*(dec31_12,us_pacific(dec31_12)), Dotted);

draw(mar_trans2*(mar9_12,gmt(mar9_5))--nov_trans2*(nov2_7, gmt(nov2)), Dotted);
draw(nov_trans2*(nov2_11,gmt(nov2_3))--dec_trans2*(dec31_23_59_58, us_pacific(dec31_23_59_58)), Dotted);

pair pos_spring_dst_switch = mar_trans1*(mar9_12,gmt(mar9-36hour));
pair pos_autumn_dst_switch = nov_trans1*(nov2_12,gmt(nov2-36hour));
pair pos_leap_second = dec_trans1*(jan1+0.5,gmt(dec31_12-36hour));
label("Spring DST Switch", pos_spring_dst_switch, RightSide, fontsize(18));
label("Autumn DST Switch", pos_autumn_dst_switch, RightSide, fontsize(18));
label("Leap Second", pos_leap_second, RightSide, fontsize(18));

label("PST8PDT in 2008", (pos_leap_second.x,pos_spring_dst_switch.y), fontsize(18));

Label xlabel = shift(8mm,-3mm)*rotate(-30)*Label();
Label ylabel = shift(0mm,0.3mm)*rotate(0)*Label();

picture mar_pic1;
filldraw(mar_pic1, shift(mar9_8, gmt(mar9))*scale(4hour,5lhour)*unitsquare, palecyan, black);
draw(mar_pic1, graph(gmt,mar9,mar10));
draw(mar_pic1, graph(jp,mar9,mar10));
draw(mar_pic1, graph(us_pacific,mar9,mar9_10-1));
draw(mar_pic1, graph(us_pacific,mar9_10,mar10));
real[] xmajor = {mar9,mar9_12,mar10};
real[] xminor; for (int i = 0; i < 24; ++i) xminor.push(mar9 + i * hour);
xaxis(mar_pic1, YEquals(gmt(mar8_13),false), LeftTicks(xlabel, xticklabel, Ticks=xmajor, ticks=xminor));
real[] ymajor = {gmt(mar9),gmt(mar9_12),gmt(mar10),gmt(mar10_12)};
real[] yminor; for (int i = -12; i < 36; ++i) yminor.push(gmt(mar9) + i * lhour);
yaxis(mar_pic1, XEquals(mar9,false), RightTicks(ylabel, yticklabel, Ticks=ymajor, ticks=yminor));
label(mar_pic1, "JST", (mar9_20, jp(mar9_20)), E);
label(mar_pic1, "GMT", (mar9_20, gmt(mar9_20)), E);
label(mar_pic1, "PDT", (mar9_20, us_pacific(mar9_20)), E);
label(mar_pic1, "PST", (mar9_6, us_pacific(mar9_6)), 3S);
unitsize(mar_pic1, usize1x, usize1y);
add((shift(-mar_orig1)*mar_pic1).fit(),mar_pos1);

picture mar_pic2;
draw(mar_pic2, graph(us_pacific,mar9_8,mar9_10-1));
draw(mar_pic2, graph(us_pacific,mar9_10,mar9_12));
fill(mar_pic2, shift(mar9_10, gmt(mar9_3))*scale(1mm/usize2x, 1mm/usize2y)*unitcircle);
filldraw(mar_pic2, shift(mar9_10, gmt(mar9_2))*scale(1mm/usize2x, 1mm/usize2y)*unitcircle, white, black);
real[] xmajor = {mar9_8,mar9_9,mar9_10,mar9_11,mar9_12};
xaxis(mar_pic2, YEquals(gmt(mar9),false), LeftTicks(xlabel, xticklabel, Ticks=xmajor));
real[] ymajor = {gmt(mar9),gmt(mar9_1),gmt(mar9_2),gmt(mar9_3),gmt(mar9_4),gmt(mar9_5)};
yaxis(mar_pic2, XEquals(mar9_8,false), RightTicks(ylabel, yticklabel, Ticks=ymajor));
label(mar_pic2, "PST", (mar9_9+1800, us_pacific(mar9_9+1800)), SE);
label(mar_pic2, "PDT", (mar9_11, us_pacific(mar9_11)), NW);
unitsize(mar_pic2, usize2x, usize2y);
add((shift(-mar_orig2)*mar_pic2).fit(),mar_pos2);

picture nov_pic1;
filldraw(nov_pic1, shift(nov2_7, gmt(nov2))*scale(4hour,3lhour)*unitsquare, palecyan, black);
draw(nov_pic1, graph(gmt,nov2,nov3));
draw(nov_pic1, graph(jp,nov2,nov3));
draw(nov_pic1, graph(us_pacific,nov2,nov2_9-1));
draw(nov_pic1, graph(us_pacific,nov2_9,nov3));
real[] xmajor = {nov2,nov2_12,nov3};
real[] xminor; for (int i = 0; i < 24; ++i) xminor.push(nov2 + i * hour);
xaxis(nov_pic1, YEquals(gmt(nov1_13),false), LeftTicks(xlabel, xticklabel, Ticks=xmajor, ticks=xminor));
real[] ymajor = {gmt(nov2),gmt(nov2_12),gmt(nov3),gmt(nov3_12)};
real[] yminor; for (int i = -12; i < 36; ++i) yminor.push(gmt(nov2) + i * lhour);
yaxis(nov_pic1, XEquals(nov2,false), RightTicks(ylabel, yticklabel, Ticks=ymajor, ticks=yminor));
label(nov_pic1, "JST", (nov2_20, jp(nov2_20)), E);
label(nov_pic1, "GMT", (nov2_20, gmt(nov2_20)), E);
label(nov_pic1, "PST", (nov2_20, us_pacific(nov2_20)), E);
label(nov_pic1, "PDT", (nov2_5, us_pacific(nov2_5)), 3S);
unitsize(nov_pic1, usize1x, usize1y);
add((shift(-nov_orig1)*nov_pic1).fit(),nov_pos1);

picture nov_pic2;
draw(nov_pic2, graph(us_pacific,nov2_7,nov2_9-1));
draw(nov_pic2, graph(us_pacific,nov2_9,nov2_11));
fill(nov_pic2, shift(nov2_9, gmt(nov2_1))*scale(1mm/usize2x, 1mm/usize2y)*unitcircle);
filldraw(nov_pic2, shift(nov2_9, gmt(nov2_2))*scale(1mm/usize2x, 1mm/usize2y)*unitcircle, white, black);
real[] xmajor = {nov2_7,nov2_8,nov2_9,nov2_10,nov2_11};
xaxis(nov_pic2, YEquals(gmt(nov2),false), LeftTicks(xlabel, xticklabel, Ticks=xmajor));
real[] ymajor = {gmt(nov2),gmt(nov2_1),gmt(nov2_2),gmt(nov2_3)};
yaxis(nov_pic2, XEquals(nov2_7,false), RightTicks(ylabel, yticklabel, Ticks=ymajor));
label(nov_pic2, "PDT", (nov2_8+1800, us_pacific(nov2_8+1800)), NW);
label(nov_pic2, "PST", (nov2_10, us_pacific(nov2_10)), SE);
unitsize(nov_pic2, usize2x, usize2y);
add((shift(-nov_orig2)*nov_pic2).fit(),nov_pos2);

picture dec_pic1;
draw(dec_pic1, graph(gmt,dec31_12,jan1_12));
draw(dec_pic1, graph(jp,dec31_12,jan1_12));
draw(dec_pic1, graph(us_pacific,dec31_12,jan1_12));
draw(dec_pic1, (dec31_23_59_60,gmt(dec31))--(dec31_23_59_60,gmt(jan1+21hour)), lightred);
draw(dec_pic1, (dec31_12,gmt(dec31_23_59_60))--(jan1_12,gmt(dec31_23_59_60)), lightred);
draw(dec_pic1, (dec31_12,jp(dec31_23_59_60))--(jan1_12,jp(dec31_23_59_60)), lightred);
draw(dec_pic1, (dec31_12,us_pacific(dec31_23_59_60))--(jan1_12,us_pacific(dec31_23_59_60)), lightred);
dot(dec_pic1, (dec31_23_59_60, gmt(dec31_23_59_60)), brown);
dot(dec_pic1, (dec31_23_59_60, jp(dec31_23_59_60)), brown);
dot(dec_pic1, (dec31_23_59_60, us_pacific(dec31_23_59_60)), brown);
real[] xmajor = {dec31_12,jan1,jan1_12};
real[] xminor; for (int i = 0; i < 24; ++i) xminor.push(dec31_12 + i * hour);
xaxis(dec_pic1, YEquals(gmt(dec31_1),false), LeftTicks(xlabel, xticklabel, Ticks=xmajor, ticks=xminor));
real[] ymajor = {gmt(dec31_12), gmt(jan1),gmt(jan1_12),gmt(jan2)};
real[] yminor; for (int i = -12; i < 36; ++i) yminor.push(gmt(dec31_12) + i * lhour);
yaxis(dec_pic1, XEquals(dec31_12,false), RightTicks(ylabel, yticklabel, Ticks=ymajor, ticks=yminor));
label(dec_pic1, "JST", (jan1_8, jp(jan1_8)), E);
label(dec_pic1, "GMT", (jan1_8, gmt(jan1_8)), E);
label(dec_pic1, "PST", (jan1_8, us_pacific(jan1_8)), E);
unitsize(dec_pic1, usize1x, usize1y);
add((shift(-dec_orig1)*dec_pic1).fit(),dec_pos1);

picture dec_pic2;
fill(dec_pic2, shift(dec31_23_59_60, us_pacific(dec31_23_59_58))*scale(1, jan1_0_0_2-dec31_23_59_58)*unitsquare, lightred);
draw(dec_pic2, graph(us_pacific,dec31_23_59_58,jan1_0_0_2));
draw(dec_pic2, graph(us_pacific,dec31_23_59_60,jan1),brown+2);
real[] xmajor = {dec31_23_59_58, dec31_23_59_59, dec31_23_59_60, jan1, jan1_0_0_1, jan1_0_0_2};
xaxis(dec_pic2, YEquals(gmt(dec31_15_59_58),false), LeftTicks(xlabel, xticklabel, Ticks=xmajor));
real[] ymajor = {gmt(dec31_15_59_58),gmt(dec31_15_59_59),gmt(dec31_15_59_59)+1,gmt(dec31_16),gmt(dec31_16_0_1),gmt(dec31_16_0_2)};
yaxis(dec_pic2, XEquals(dec31_23_59_58,false), RightTicks(ylabel, yticklabel, Ticks=ymajor));
label(dec_pic2, "leap second", (dec31_23_59_60+0.5, us_pacific(dec31_23_59_60+0.5)), E);
label(dec_pic2, "PST", (jan1_0_0_1, us_pacific(jan1_0_0_1)), SE);
unitsize(dec_pic2, usize3x, usize3y);
add((shift(-dec_orig2)*dec_pic2).fit(),dec_pos2);

picture dec_pic3;
;
draw(dec_pic3, graph(gmt,dec31_15_59_58,dec31_16-0.01));
draw(dec_pic3, graph(gmt,dec31_16,dec31_16_0_2));
fill(dec_pic3, shift(dec31_16, gmt(dec31_16))*scale(1mm/usize3x, 1mm/usize3y)*unitcircle);
filldraw(dec_pic3, shift(dec31_16, (gmt(dec31_16))-1)*scale(1mm/usize3x, 1mm/usize3y)*unitcircle, white, black);
real[] xmajor = {dec31_15_59_58, dec31_15_59_59, dec31_16, dec31_16_0_1, dec31_16_0_2};
xaxis(dec_pic3, YEquals(gmt(dec31_15_59_58),false), LeftTicks(xlabel, xticklabel, Ticks=xmajor));
real[] ymajor = {gmt(dec31_15_59_58),gmt(dec31_15_59_59),gmt(dec31_15_59_59)+1,gmt(dec31_16),gmt(dec31_16_0_1),gmt(dec31_16_0_2)};
yaxis(dec_pic3, XEquals(dec31_15_59_58,false), RightTicks(ylabel, yticklabel, Ticks=ymajor));
label(dec_pic3, "GMT", (dec31_16_0_1, gmt(dec31_16_0_1)), SE);
unitsize(dec_pic3, usize3x, usize3y);
add((shift(-dec_orig3)*dec_pic3).fit(),dec_pos3);

天泣記

2009-07-01 (Wed)

2009-07-02 (Thu)

2009-07-03 (Fri)

2009-07-04 (Sat)

2009-07-05 (Sun)

2009-07-06 (Mon)

2009-07-09 (Thu)

2009-07-18 (Sat)

2009-07-19 (Sun)

2009-07-20 (Mon)

2009-07-21 (Tue)

2009-07-22 (Wed)

2009-07-23 (Thu)

2009-07-24 (Fri)

2009-07-28 (Tue)

2009-07-29 (Wed)

2009-07-30 (Thu)