星座早見プログラムを少し変更して,赤緯45°以上の天の北極の歳差運動による変化を見られるようにしました。より短いプログラムですので,実行ファイル形式でも利用してもらえると思います。古代エジプトのピラミッド建設で,正確な東西南北方向を決める方法として紀元前25000年前には,こぐま座のコカブと北斗七星のミザールを結ぶ線が使われたという説を確かめたり出来ます。

実行ファイルは→こちらからダウンロード

ソースコードをパイソンで実行する場合には,同じフォルダーに恒星,星座名,星座線データのCSVファイルを入れて実行します。ファイルは→こちらにあります。

'''
******************************************************
Star chart program by Python(Archaeoastronomy version)
2019.08.18 Ver1.0 Archaeoastronomy version1.2
2020.03.13  from Ver.2.0 Execute file form
星座早見プログラムをもとに,天の北極の星々の変化を表示するプログラムへ
******************************************************
'''

import os
import tkinter as tk
import csv
import math

# Calculation to cure Julius day in the Gregorian calendar
def JDT(jd):
    Z = int(jd + 0.5)
    if Z >= 2299161:
        a = int((Z - 1867216.25) / 36524.25)
        A = Z + 1 + a - int(a / 4)
    else:
        A = Z
    B = A + 1524
    C = int((B - 122.1) / 365.25)
    K = int(365.25 * C)
    E = int((B - K)/30.6001)
    D = B - K - int(30.6001 * E) + (jd + 0.5) - int(jd + 0.5)
    if E < 13.5:
        M = E -1
    else:
        M = E -13
    if M > 2.5:
        Y = C - 4716
    else:
        Y = C - 4715
    if M >= 13:
        Y = Y + 1
        M = M -12
    if Y <= 0:
        Y = Y -1

    h = D - int(D)
    D = int(D)
    h = h*24
    lh = round(long/15)
    h = h +lh
    if h >= 24.0:
        h = h - 24
        D = D + 1
    if D >= 32:
        D = D - 31
        M = M + 1
    if M >= 13:
        Y = Y + 1
        M = M - 12
    h = round(h,1)
    return [Y,M,D,h]

def koseji(jd,long):    # Calculate sidereal time
    B = jd - 2415020.0
    R = 366.2422/365.2422
    ST = 18.6461 + 24*B*R + 3.24e-14*B*B + long/15
    ST = 24*(ST/24-int(ST/24))
    if ST < 0:
        ST = ST + 24
    return ST

def yogen_AD(arfa,drta):  # Calculation of the direction cosine
    a = math.radians(arfa)
    d = math.radians(drta)
    L = math.cos(a) * math.cos(d)
    M = math.cos(d) * math.sin(a)
    N = math.sin(d)
    return [L,M,N]

def horizon(ad):   # to horizontal coordinate(ST;sideral time LAT;latitude ad; direction cosine)
    sT = math.sin(math.radians(ST))
    cT = math.cos(math.radians(ST))
    # sL = math.sin(math.radians(LAT))  ** = 1
    # cL = math.cos(math.radians(LAT))  ** = 0
    L = cT*ad[0] + sT*ad[1] # - cL*ad[2]
    M = -sT*ad[0] + cT*ad[1]
    N = ad[2] # + cL*sT*ad[1] + cL*cT*ad[0]
    if L == 0:
        L = 0.01
    h = math.asin(N)
    h = math.degrees(h)
    A = math.atan(-M/L)
    A = math.degrees(A)
    if L < 0. :
        A = A + 180.
    return [h,A]

def proper_move(RA,DC,V1,V2):   # proper monve
    T  = Cent - 1
    RA = RA + V1*T/3600000/math.cos(math.radians(DC))
    DC = DC + V2*T/3600000
    return [RA , DC]

def saisa_hosei(ad):  # 歳差補正
    t = Cent - 1
    f = 0.640616 * t + 0.0000839* t*t + 0.000005*t**3
    z = 0.640616 * t + 0.000304 * t*t + 0.00000506*t**3
    s = 0.556753 * t - 0.000119 * t*t - 0.0000116*t**3
    sF = math.sin(math.radians(f))
    cF = math.cos(math.radians(f))
    sZ = math.sin(math.radians(z))
    cZ = math.cos(math.radians(z))
    sS = math.sin(math.radians(s))
    cS = math.cos(math.radians(s))
    L = (- sZ*sF + cZ*cS*cF)*ad[0] + (- sZ*cF - cZ*cS*sF)*ad[1] - cZ*sS*ad[2]
    M = (cZ*sF + sZ*cS*cF)*ad[0] + (cZ*cF - sZ*cS*sF)*ad[1] - sZ*sS*ad[2]
    N = sS*cF*ad[0] - sS*sF*ad[1] + cS*ad[2]
    return [L, M, N]

def dispXY(hh, AA):
    dot = 480  # 画面のパラメータ 中心座標(540,375)
    r = dot * math.sin(math.radians(90-hh))
    x = -r * math.sin(math.radians(AA)) + 540
    y = -r * math.cos(math.radians(AA)) + 375
    return [x, y]

def sekido():  # 赤経赤緯の描画
    xl = [0.0]*182
    yl = [0.0]*182
    rad = [0.0]*182
    for i in range(180):
        rad[i] = i*2

    for i in range(181):
        xy= dispXY(80,rad[i])
        xl[i] ,yl[i]= xy

    for i in range(180):
        app.co_line(xl[i], yl[i], xl[i+1], yl[i+1], '#710071')

    for i in range(181):
        xy= dispXY(60,rad[i])
        xl[i] ,yl[i] = xy

    for i in range(180):
        app.co_line(xl[i], yl[i], xl[i+1], yl[i+1], '#710071')

    xy = dispXY(45,ST)
    x1 ,y1 = xy
    xy = dispXY(45,ST+180)
    x2, y2 = xy
    app.co_line(x1, y1, x2, y2, '#710071')
    app.text_in(x1 + 5, y1, "0h", "white")
    app.text_in(x2 + 5, y2, "12h", "white")
    xy = dispXY(45, ST+90)
    x1, y1 = xy
    xy = dispXY(45, ST+270)
    x2, y2 = xy
    app.co_line(x1, y1, x2, y2, '#710071')
    app.text_in(x1 + 5, y1, "18h", "white")
    app.text_in(x2 + 5, y2, "6h", "white")

    app.co_line(540, 38, 540, 710, "#999900")


def star_color(CL):    # 恒星の色 CL:color index
    if CL < -0.16:
        c = "#a09eff"
    elif CL < 0.15:
        c = "#a0d7ff"
    elif CL < 0.45:
        c = "#d7e8ff"
    elif CL < 0.68:
        c = "#ffffff"
    elif CL < 1.15:
        c = "#ffffdc"
    elif CL < 1.6:
        c = "#ffe6aa"
    else:
        c = "#ffd7b1"
    return c

def magnitude(m):    # 等級を星の大きさに
    if m < 0.0:
       rd = 7
    elif m < 1.0:
       rd = 6
    elif m < 2.0:
       rd = 5
    elif m < 3.0:
       rd = 4
    elif m < 4.0:
        rd = 3
    else:
       rd = 2
    return rd

class Conline:  # 星座線クラス
    def __init__(self):
        self.lcnum = 0
        self.linRas = 0
        self.linDcs = 0
        self.linV1s = 0
        self.linV2s = 0
        self.linRae = 0
        self.linDce = 0
        self.linV1e = 0
        self.linV2e = 0

class Star:  # 星クラス
    def __init__(self):
        self.stnum = 0
        self.stV1 = 0
        self.stV2 = 0
        self.stRA = 0
        self.stDC = 0
        self.stMg = 0
        self.stCL = 0
class Constelation:  # 星座クラス
    def __init__(self):
        self.con_name = ""
        self.con_Ra = 0
        self.con_Dc = 0

class Planetarium:  # プラネタリウムクラス

    # 星表示のための変数
    XX = [0.0] * 1000
    YY = [0.0] * 1000
    hh = [0.0] * 1000
    AA = [0.0] * 1000
    MG = [0.0] * 1000
    CLs = [0.0] * 1000
    rd = [0] * 1000
    hLs = [0.0] * 500
    hLe = [0.0] * 500
    ALs = [0.0] * 500
    ALe = [0.0] * 500
    x1 = [0.0] * 500
    y1 = [0.0] * 500
    x2 = [0.0] * 500
    y2 = [0.0] * 500
    nh = [0.0] * 60
    nA = [0.0] * 60
    xn = [0.0] * 60
    yn = [0.0] * 60
    coname = [""] * 60
    star_counter = 0
    line_counter = 0
    con_counter = 0

    def __init__(self):
        # 星、星座線リストを作る
        self.star_list = []
        self.conline_list = []
        self.conste_list = []

        for i in range(1263):
            self.star_list.append(Star())

        for i in range(673):
            self.conline_list.append(Conline())

        for i in range(89):
            self.conste_list.append(Constelation())


        # データ読み込み
        with open("starData1263.csv", "r") as f:  # 恒星ファイル読み込み
            stDATA = csv.reader(f, delimiter=",")
            i = 0
            for row in stDATA:
                self.star_list[i].stnum = int(row[4])  # No.
                self.star_list[i].stV1 = float(row[1]) # 固有運動RA
                self.star_list[i].stV2 = float(row[2]) #        DC
                self.star_list[i].stRA = float(row[6]) # 赤経
                self.star_list[i].stDC = float(row[7]) # 赤緯
                self.star_list[i].stMg = float(row[5]) # 光度
                self.star_list[i].stCL = float(row[8]) # 色指数
                i += 1

        with open("cons_lineData.csv", "r") as f:  #星座線ファイル読み込み
            szL=csv.reader(f,delimiter=",")
            i = 0
            for row in szL:
                self.conline_list[i].lcnum = str(row[0]) #星座コード
                self.conline_list[i].linRas = float(row[1])  #始点
                self.conline_list[i].linDcs = float(row[2])
                self.conline_list[i].linV1s = float(row[3])
                self.conline_list[i].linV2s = float(row[4])
                self.conline_list[i].linRae = float(row[5])  #終点
                self.conline_list[i].linDce = float(row[6])
                self.conline_list[i].linV1e = float(row[7])
                self.conline_list[i].linV2e = float(row[8])
                i += 1
        with open("cons_nameData.csv", "r") as f:  # 星座名ファイル読み込み
            szM=csv.reader(f,delimiter=",")
            i = 0
            for row in szM:
                self.conste_list[i].con_name = str(row[0]) # 星座名
                self.conste_list[i].con_Ra = float(int(row[1])*15+int(row[2])*0.25)
                self.conste_list[i].con_Dc = float(row[3])
                i += 1

    def star_culc(self):       # 恒星表示メソッド

        for lin in self.conline_list:  #恒星座標 固有運動の補正
            ads = proper_move(lin.linRas,lin.linDcs,lin.linV1s,lin.linV2s)
            lin.linRas = ads[0]
            lin.linDcs = ads[1]
            ade = proper_move(lin.linRae,lin.linDce,lin.linV1e,lin.linV2e)
            lin.linRae = ade[0]
            lin.linDce = ade[1]

        m = 0      # 星座線の地平座標計算
        for line in self.conline_list:
            ad = yogen_AD(line.linRas,line.linDcs)
            ad = saisa_hosei(ad)
            hs = horizon(ad)
            if hs[0] < 45.0:
                continue            
            ad = yogen_AD(line.linRae,line.linDce)
            ad = saisa_hosei(ad)
            he = horizon(ad)
            if he[0] < 45.0:
                continue
            self.hLs[m],self.ALs[m] = hs
            self.hLe[m],self.ALe[m] = he
            m += 1
        self.line_counter = m

        i = 0
        for i in range( self.line_counter):   # 画面上の座標
            xy= dispXY(self.hLs[i],self.ALs[i])
            self.x1[i] ,self.y1[i] = xy
            xy= dispXY(self.hLe[i],self.ALe[i])
            self.x2[i] ,self.y2[i] = xy

            # 恒星の処理 ******

        for star in self.star_list:  # 固有運動の補正
            ad = proper_move(star.stRA,star.stDC,star.stV1,star.stV2)
            star.stRA = ad[0]
            star.stDC = ad[1]

        n = 0              # 恒星の地平座標計算
        for star in self.star_list:
            ad = yogen_AD(star.stRA,star.stDC)
            ad = saisa_hosei(ad)
            h = horizon(ad)
            if h[0] < 45.0:
                continue
            self.hh[n],self.AA[n] = h
            self.MG[n] = star.stMg
            self.CLs[n] = star.stCL
            n += 1
        self.star_counter = n

        for i in range(self.star_counter):
            xy= dispXY(self.hh[i],self.AA[i])
            self.XX[i] ,self.YY[i] = xy

            # 星座名の処理 *******

        n = 0    #星座名表示の地平座標計算
        for c in self.conste_list:
            ad = yogen_AD(c.con_Ra,c.con_Dc)
            ad = saisa_hosei(ad)
            h = horizon(ad)
            if h[0] < 45.0:
                continue
            self.nh[n] ,self.nA[n] = h
            self.coname[n] = c.con_name
            n += 1
        self.con_counter = n

        for i in range(self.con_counter):
            xy= dispXY(self.nh[i],self.nA[i])
            self.xn[i] ,self.yn[i] = xy


    def star_display(self):
        # 星座線を引く
        for i in range( self.line_counter):
            app.co_line(self.x1[i], self.y1[i], self.x2[i], self.y2[i], 'blue')

        # 星のプロット
        for i in range(self.star_counter):
            self.rd[i] = magnitude(self.MG[i])
            color = star_color(self.CLs[i])
            app.point_star(self.XX[i], self.YY[i], self.rd[i]/2, color)

        #星座名の表示
        for i in range(self.con_counter):
            app.text_in(self.xn[i], self.yn[i], self.coname[i], 'red')



class Time:
    def __init__(self,Y,D,lo):
        self.Ydate = Y
        self.Dtime = D
        self.LAT = 34.5
        self.JD = 0
        self.ST = 0
        self.Cent = 0
        self.YY = 0
        self.long = lo

    def Julian(self):  # (ユリウス日)の計算
        self.JD = 0
        if self.Ydate != abs(self.Ydate):
            SP1 = -1
            self.Ydate = abs(self.Ydate)
        else:
            SP1 = 1
        self.YY = int(self.Ydate/10000)
        MD = int(self.Ydate-10000*self.YY)
        MM = int(MD/100)
        DD = MD - 100 * MM
        HH = int(self.Dtime/100)
        MS = self.Dtime-100*HH
        if SP1 < 0:
            self.YY = self.YY * SP1  # +1  BC.でなく,BC1年を0年として-で入力する
        SP2 = self.YY + (MM-1)/12 + DD/365.25
        if MM <= 2 :
            MM = MM + 12
            self.YY = self.YY - 1
        if self.YY < 0:
            self.JD = math.floor(365.25*self.YY) + int(30.59*(MM-2)) + DD - self.long/360 + 1721086.5
        else:
            self.JD = int(365.25*self.YY) + int(30.59*(MM-2)) + DD - self.long/360 + 1721086.5
        if SP2 > 1582.78:  # グレゴリオ暦以降
            self.JD = self.JD + int(self.YY/400) - int(self.YY/100) + 2
        if MM > 12:
            MM = MM - 12
            self.YY = self.YY + 1
        self.JD = self.JD + HH/24 + MS/1440

        self.ST  = koseji(self.JD,self.long)*15      # ST 恒星時(度)
        self.Cent  = (self.JD - 2415021.0)/36525      # 元期1900年I.1 0.5


def Batch():
    app.refresh()
    plt.star_culc()
    plt.star_display()

def update(event):
    global LAT
    global long
    global ymd
    global ST
    global Cent
    global YY
    global JD

    gvl.getV()
    Dtime = gvl.D_time
    Dtime = float(Dtime)
    Ydate = gvl.Y_date
    Ydate = float(Ydate)
    long  = gvl.LONG
    long  = float(long)
    LAT = gvl.Lat
    LAT = float(LAT)
    tm = Time(Ydate, Dtime, long)
    tm.Julian()
    JD = tm.JD
    ST = tm.ST
    Cent = tm.Cent
    YY = tm.YY
    ymd = JDT(JD)

    Batch()

def Timeforward():
    global ST
    global JD
    global ymd
    ymd[3] = ymd[3] + 1
    ST = ST + 15
    JD = JD + (1/24)
    ymd = JDT(JD)

    Batch()

def Timeback():
    global ST
    global JD
    global ymd
    ST = ST - 15
    JD = JD - (1/24)
    ymd = JDT(JD)
    Batch()

def Dateforward():
    global ST
    global JD
    global ymd
    ST = ST + 1.002738
    JD = JD + 1
    ymd = JDT(JD)
    Batch()

def Dateback():
    global ST
    global JD
    global ymd
    ST = ST - 1.002738
    JD = JD -1
    ymd = JDT(JD)
    Batch()


class Application(tk.Frame):
    def __init__(self,master = None):
        super().__init__(master)
        master.title(u"北極星野(星図)の表示プログラム_Ver.1.1")
        master.geometry("1100x750")
        self.c0 = tk.Canvas(master=None, width = 1100, height = 750, bg="#001766")
        self.pack()
        self.create_widgets()

    def create_widgets(self):
        self.c0.create_text(955, 570, text="その日の時刻を1時間", font=('', 12), anchor='c', fill="white")
        self.c0.create_text(955, 635, text="日付を1日", font=('', 12), anchor='c', fill='white')
        self.c0.create_text(25, 667, text=" 紀元前はB.C.1年を0年として計算する。入力は-1で", font=('', 10), fill='white', anchor='w')
        self.lb1 = tk.Label(text=u"観測地の緯度・経度(南緯は-,東経は+)", fg='white', bg="#001766", anchor='w', font=('', 14)).place(x=730,
                                                                                      y=30)
        self.inputbox1 = tk.Entry(width=4, font=('', 14))
        self.inputbox1.place(x=820, y=60)
        self.inputbox1.insert(tk.END, "34.5")
        self.inputbox1.bind("<KeyPress-Return>", update)
        self.inputbox2 = tk.Entry(width=6, font=('', 14))
        self.inputbox2.place(x=900, y=60)
        self.inputbox2.insert(tk.END, "135.8")
        self.inputbox2.bind("<KeyPress-Return>", update)
        self.lb2 = tk.Label(text="日付(YYYYMMDD)  時刻(hhmm)", fg='white', bg="#001766", anchor='w', font=('', 14)).place(x=800,
                                                                                                                 y=95)
        self.inputbox3 = tk.Entry(width=10, font=('', 14))
        self.inputbox3.place(x=850, y=120)
        self.inputbox3.insert(tk.END, "6970822")
        self.inputbox3.bind("<KeyPress-Return>", update)
        self.inputbox4 = tk.Entry(width=6, font=('', 14))
        self.inputbox4.place(x=980, y=120)
        self.inputbox4.insert(tk.END, "2100")
        self.inputbox4.bind("<KeyPress-Return>", update)

        self.lb3 = tk.Label(text=u'入力→Enterで再表示します', width=25, font=('', 10), fg='#001766', bg='skyblue').place(x=870, y=165)
        self.Button1 = tk.Button(text=u'進む>>', command=Timeforward, width=10, font=('', 10)).place(x=960, y=590)
        self.Button2 = tk.Button(text=u'<<もどる', command=Timeback, width=10, font=('', 10)).place(x=860, y=590)
        self.Button3 = tk.Button(text=u'進む>', command=Dateforward, width=10, font=('', 10)).place(x=960, y=650)
        self.Button4 = tk.Button(text=u'<もどる', command=Dateback, width=10, font=('', 10)).place(x=860, y=650)

    def refresh(self):  # 再表示
        global img_moon   # create_imageに必要
        self.c0.create_oval(170, 10, 910, 745, fill="#001766", outline="#001766")
        self.c0.create_text(540, 20, text="S", font=('', 18), fill="yellow")
        self.c0.create_text(540, 730, text="N", font=('', 18), fill="yellow")
        self.c0.create_text(180, 375, text="W", font=('', 18), fill="yellow")
        self.c0.create_text(900, 375, text="E", font=('', 18), fill="yellow")
        self.c0.create_oval(200, 35, 880, 715, fill='black', outline='skyblue', width=1)
        self.c0.pack()

        if ymd[0] <= 0:
            self.lb4 = tk.Label(text=f"B.C. {-ymd[0]}年 {ymd[1]}月 {ymd[2]}日   {ymd[3]} 時の北極    ", fg='white', \
                           bg="#001766", anchor='w', font=('', 14))
        else:
            self.lb4 = tk.Label(text=f"A.D. {ymd[0]}年 {ymd[1]}月 {ymd[2]}日   {ymd[3]} 時の北極    ", fg='white', \
                       bg="#001766", anchor='w', font=('', 14))
        self.lb4.place(x=25, y=25)
        self.lb5 = tk.Label(text=f"緯度= {LAT}°経度= {long}°", fg='white', bg="#001766", anchor \
            ='w', font=('', 14))
        self.lb5.place(x=25, y=60)
        # self.c0.create_oval(185, 35, 1015, 865, fill='black', outline='skyblue', width=1)
        self.lb6 = tk.Label(text=f" ユリウス日は = {round(JD, 3)} 日     ", fg='white', bg="#001766", anchor \
            ='w', font=('', 10))
        self.lb6.place(x=25, y=600)
        self.lb7 = tk.Label(text=f" 恒星時は = {round(ST / 15, 1)}  時    ", fg='white', bg="#001766", anchor \
            ='w', font=('', 10))
        self.lb7.place(x=25, y=630)
        self.lb8 = tk.Label(text=f" 赤緯45°以上を示す 内側の円は赤緯80° 外側は60°", fg='white', bg="#001766", anchor \
                            ='w', font=('', 10))
        self.lb8.place(x=25, y=690)

        sekido()

    def point_star(self, X, Y, dir, color):
        self.c0.create_oval(X-dir, Y-dir, X+dir, Y+dir, fill = color)
        self.c0.pack()

    def co_line(self, X1, Y1, X2, Y2, color):
        self.c0.create_line(X1, Y1, X2, Y2, fill = color, smooth = "True")
        self.c0.pack()

    def text_in(self, X, Y, name, color):
        self.c0.create_text(X, Y, text = name, font = ('',10), fill = color)
        self.c0.pack()

    def point_eclip(self, X, Y, dir, color):
        self.c0.create_oval(X-dir, Y-dir, X+dir, Y+dir, outline = color, )
        self.c0.pack()

class Get_Value(Application):
    def __init__(self,master=None):
        super().__init__(master)
        self.Lat = 0
        self.LONG = 0
        self.Y_date = 0
        self.D_time = 0
    def getV(self):
        self.Lat = self.inputbox1.get()
        self.LONG = self.inputbox2.get()
        self.Y_date = self.inputbox3.get()
        self.D_time = self.inputbox4.get()


plt = Planetarium()
root = tk.Tk()
root.resizable(width=False, height=False)
app = Application(master=root)   # tkinter(GUI)
gvl = Get_Value(master=root)    # 設定変更値の取得

# 日時の計算
tm = Time(6970822, 2100, 135.8) # 初期設定(年月日,時刻,経度
tm.Julian()
JD = tm.JD
ST = tm.ST
Cent = tm.Cent
YY = tm.YY
long = tm.long  # 経度
LAT = tm.LAT  # 緯度
ymd = JDT(JD)
Batch()

if __name__ == '__main__':
    app.mainloop()