遅ればせながら、北口さん（あれ？南口さんのほうがいいのかな）が行っている月の分析をマネてみました。
くりちゃんさんの方式で天体を撮ったことがなかったので、このとき初めて意識して2分（結果的には2分2秒でした）おいて比較明合成してみました。
2013/6/4 3:03:22 → 3:05:24
それぞれの月を円とみなし、周縁部の座標を3点ずつ読み、中心座標を求めます。
右上の月から読み取った3点の座標
a=(1310, 127), b=(2022, 1315), c=(2145, 257)
3つの座標を通る円の中心の計算：公式があったのでパクリました＾＾
α=BC^2*(AC^2 + AB^2 - BC^2)
={(2022 - 2145)^2 + (1315 - 257)^2} * [{(1310 - 2145)^2 + (127 - 257)^2} + {..} - {..}]
=1699379754560
β=AC^2*(BC^2 + AB^2 - AC^2)
=1670092716000
γ=AB^2*(BC^2 + AC^2 - AB^2)
=-133647124960
x =(ax*α + bx*β + cx*γ)/(α+β+γ)
=(1310*α+2022*β+2145*γ)/3235825345600
=1642.994
y =(ay*α + by*β + cy*γ)/(α+β+γ)
=734.788
r =693.030
左下の月から読み取った3点の座標
(450, 1150), (1122, 2362), (1435, 2001)
（途中略）
(787.066, 1755.409)
r =692.917
画像上の月の半径は約693pixel、たいてい１４００ｘ１４００に収まりますがこの日も同様でした。
画像上で動く方向は、右に855.928(=1642.994-787.066)、上に1020.621。
傾きは atan(1020.621/855.928) = 50.0156°
この傾きが地球の自転軸に垂直方向のはずなので、右上の月を右に50度回転させると地球の自転軸を垂直とした場合の月の画像になります。
国立天文台によると、この日3:05のパラメータは、P=-23.84, B0=-3.13, L0=6.46 でした。
画像を回転するのは面倒なのでサボって、経緯線を書くプログラムのほうで補正をしてみます。（50.0156-23.84するだけです＾＾）
くりちゃんさんの方式で天体を撮ったことがなかったので、このとき初めて意識して2分（結果的には2分2秒でした）おいて比較明合成してみました。
2013/6/4 3:03:22 → 3:05:24
それぞれの月を円とみなし、周縁部の座標を3点ずつ読み、中心座標を求めます。
右上の月から読み取った3点の座標
a=(1310, 127), b=(2022, 1315), c=(2145, 257)
3つの座標を通る円の中心の計算：公式があったのでパクリました＾＾
α=BC^2*(AC^2 + AB^2 - BC^2)
={(2022 - 2145)^2 + (1315 - 257)^2} * [{(1310 - 2145)^2 + (127 - 257)^2} + {..} - {..}]
=1699379754560
β=AC^2*(BC^2 + AB^2 - AC^2)
=1670092716000
γ=AB^2*(BC^2 + AC^2 - AB^2)
=-133647124960
x =(ax*α + bx*β + cx*γ)/(α+β+γ)
=(1310*α+2022*β+2145*γ)/3235825345600
=1642.994
y =(ay*α + by*β + cy*γ)/(α+β+γ)
=734.788
r =693.030
左下の月から読み取った3点の座標
(450, 1150), (1122, 2362), (1435, 2001)
（途中略）
(787.066, 1755.409)
r =692.917
画像上の月の半径は約693pixel、たいてい１４００ｘ１４００に収まりますがこの日も同様でした。
画像上で動く方向は、右に855.928(=1642.994-787.066)、上に1020.621。
傾きは atan(1020.621/855.928) = 50.0156°
この傾きが地球の自転軸に垂直方向のはずなので、右上の月を右に50度回転させると地球の自転軸を垂直とした場合の月の画像になります。
国立天文台によると、この日3:05のパラメータは、P=-23.84, B0=-3.13, L0=6.46 でした。
画像を回転するのは面倒なのでサボって、経緯線を書くプログラムのほうで補正をしてみます。（50.0156-23.84するだけです＾＾）
この背景を透明にして、重ねてみると。
やってみてわかりましたが、輪郭円か、円の中心を示す＋（クロス）あたりが欲しいですね。
月の下側（南極側）が手前に来ていて北極より光っているのがわかりました。
これは1/4に縮小している画像なので、実際はもっと高解像です。
重ね合わせによってクレータの名前や動きがよくわかるようになると思います。
今回やったこと。まとめ。
・2分差の画像を比較明合成する。
・画像に合わせた経緯線を作成する。
・国立天文台から諸パラメータをひいてくる（将来的には観測でわかるようになる）。
・画像を重ね合わせる。
やっとスタートラインに立てたような気がします＾＾。
（追記＆修正版）
・緯度・経度0度をクロスで表示
・地心から月の中心を見た点を描画
これによってバグが見つかりました。
OpenGLでは軸による回転glRotateを呼ぶ順番に気をつけなければなりません。
今回、安易にL0分回転→B0分回転→P分回転の順に書いてしまいました。（ど忘れ；；）
最終的に変換されたところに描画をして戻ってくるので（なんといったらいいのか、変換を Push、Push、Push、していって最後に経緯線を書いてPop、Pop、Pop のイメージ）、自分が意図した変換をさせたい場合は逆に書くのが正解なのでした。
※いまUploadされている keisen.exe は、このバグ入りです。
正しい結果は以下です。
$ keisen.exe 6.46 -3.13 26.1756 347 で作成。
オーバーレイしたものがこれです。今度はガイドの円の中心があるので楽々重ね合わせです。
