はじめに

2023年以降海上保安庁から天測暦が発行されないことが発表されています。
この手法では海上保安庁が無償で公開しているデータをもとに計算を行なっているため、2023年以降この手法が利用可能か不明です。

計算方法

大まかな計算方法は必要なデータとともに以下のページの「解説」というリンクにまとめられています。
https://www1.kaiho.mlit.go.jp/KOHO/
こちらに載っている数式は自分で確認して解釈してください。
このページでは地球時 - 世界時がYYYY年ではN秒と大まかな説明されていません。
そこで計算方法を調べたところNASAのホームページに記載されていました。
2005年から2050年までは以下の式で求めます。
$ y = year + \frac{month - 0.5}{12} $
$ t = y - 2000 $
$ ΔT = 62.92 + 0.32217 \times t + 0.005589 \times t^2 $
https://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html
これにより時角と赤緯を求めることができ、天測暦は不要になりました。

次に第一改正から第三改正までを行う必要があります。
第一改正は天測計算表、巻末のⅡページに記載されているとおりですが、大気差の求め方が示されていません。おそらく海上保安庁では実際に観測した時の値を参考にして計算してると考えられますが、式で求められないのはあまり美しくないため、一般式を探したところ以下のページに記載されていました。
https://eco.mtk.nao.ac.jp/koyomi/wiki/C2E7B5A42FB6FEC0DE.html
$ R(ha) = \frac{0.0167^\circ}{\tan\frac{ha + 7.31}{ha + 4.4}} $
これで求めることが出来ます。haは測高度から眼高差を引いた視高度という値を代入します。
第二改正は太陽の下辺観測とすると以下の式で求められます。

Corr2 = \left\{
\begin{array}{ll}
(month - 7) \times 0.1  & (month \geq 7) \\
(6 - month) \times 0.1  & (month \lt 7)
\end{array}
\right.

第三改正は海気温差に0.2をかけることで求めることが可能です。

コード

これらをPythonで計算するコードを作成しました。

Monring_sight.py

import math

Year, Month, Day = (int(x) for x in input("Date (UTC) YYYY:MM:DD >>> ").split(":"))
Hour, Minute, Second = (int(x) for x in input("Time (UTC) HH:MM:SS >>> ").split(":"))
lat = input("lat XX XX.X >>> ").split()
lat = float(lat[0]) + float(lat[1])/60
Long = input("Long XXX XX.X >>> ").split()
Long = float(Long[0]) + float(Long[1])/60
Height = int(input("Height >>> "))
Ao = input("Observe Altitude XX XX.X >>> ").split()
Ao = float(Ao[0]) + float(Ao[1])/60
Dtemp = float(input("Sea-Air temperature difference >>> "))

a = 120 #4/30から9/1まで
b = 244 #4/30から9/1まで

P = Month - 1
Q = int((Month + 7) / 10)
Y = int((Year / 4)-int(Year / 4) + 0.77)
S = int(P * 0.55 - 0.33)
F = Hour / 24 + Minute / 1440 + Second / 86400
Ty = Year + (Month - 0.5) / 12
Tt = Ty - 2000
Dt = 62.92 + 0.32217 * Tt + 0.005589 * Tt ** 2

t = 30 * P + Q * (S - Y) + P * (1 - Q) + Day + F + Dt/86400
Seta = math.degrees(math.acos((2 * t - (a + b))/(b - a)))

t1 = 30 * P + Q * (S - Y) + P * (1 - Q) + Day + F
Seta1 = math.degrees(math.acos((2 * t1 - (a + b))/(b - a)))

#4/30から9/1までのみ有効
RA = [6.617773, 4.135644, -0.041829, -0.040668, 0.00536, 0.003267, -0.000421, -0.00012, 0.000017, 0.000015, -0.000021, -0.000064, 0.000028, 0.000047, -0.000014, 0.000017, 0.000005, 0.000007]
Dec = [17.16569, -3.37392, -5.79007, 0.21379, 0.1598, -0.01083, -0.00504, 0.00078, 0.00005, 0.00003, 0.00016, -0.00001, 0.00011, -0.00006, -0.00014, 0.00004, 0.00004, -0.00001]
R = [18.601928, 4.07404, -0.000007, -0.000004, 0.000002, 0, -0.000003, 0.000001]

def get_value(Seta, Datas):
    value = 0
    for i, j in enumerate(Datas):
        value += j * math.cos(math.radians(i * Seta))
    return value

RA_value = get_value(Seta, RA)
Dec_value = get_value(Seta, Dec)
R_value = get_value(Seta1, R)

h = R_value - RA_value
LHA = (h + Hour + Minute / 60 + Second / 3600) * 15 + Long
print("LHA  {}°{}′{}″".format(int(LHA),int(LHA % 1 *60), int(LHA % 1 * 60 % 1 * 60)))
print("d  {}°{}′".format(int(Dec_value),int(Dec_value % 1 * 60)))

Ac = math.asin(math.sin(math.radians(lat)) * math.sin(math.radians(Dec_value)) + math.cos(math.radians(lat)) * math.cos(math.radians(Dec_value)) * math.cos(math.radians(LHA)))
print("Ac  {}°{}′".format(int(Ac), int(Ac % 1 * 600 + 0.5)/10))
Z = math.acos((math.sin(math.radians(Dec_value)) - math.sin(math.radians(lat)) * math.sin(math.radians(Ac)))/(math.cos(math.radians(lat)) * math.cos(math.radians(Ac))))
print("Z  {}°".format(int(Z * 100 + 0.5)/100))

Dip = 1.776 / 60 * math.sqrt(Height)
Appalt = Ao - Dip
Ref = 0.0167 / math.tan(math.radians(Appalt + (7.31 / (Appalt + 4.4))))
Par = 9 * math.cos(math.radians(Ao))
SD = 15/60 + 45 / 3600
Corr1 = -Dip -Ref + Par +SD
if(Month > 6):
    Corr2 = Month - 7
else:
    Corr2 = 6 - Month
Corr3 = Dtemp * 0.2

Ac1 = Ac - Corr1 - Corr2 - Corr3
print("Ac  {}°{}′".format(int(Ac1), int(Ac1 % 1 * 600 + 0.5)/100))
I = Ao - Ac
print("I  {}′".format(int(I % 1 * 600 + 0.5)/100))