二十八宿距星位置计算

该程序适用于其他恒星位置的计算，需要录入基本数据，包括J2000.0的恒星位置、恒星自行等。 

# 恒星位置计算
import math
import ephem


# 二十八宿距星数据结构
class Star:
	def __init__(self, ChineseName, Designation, RightAscension, Declination, ProperMotion_ra,
	             ProperMotion_dec):
		# 需初始化的数据项：中文名、英文名、J2000.0赤经赤纬、周年视差、恒星自行
		self.name = ChineseName
		self.design = Designation
		self.__RA = RightAscension  # J2000.0赤经（单位：时分秒），不可更改
		self.__Dec = Declination  # J2000.0赤纬（单位：度分秒），不可更改
		self.pm_ra = ProperMotion_ra  # 赤经方向的自行（单位：毫角秒/年）
		self.pm_dec = ProperMotion_dec  # 赤纬方向的自行（单位：毫角秒/年）
		self.ra = 0  # 赤经计算用值
		self.dec = 0  # 赤纬计算用值

	# 读取赤经赤纬并转为计算用的度数
	def ra2deg(self):
		h, m, s = self.__RA.split(' ')
		deg = (int(h) + int(m) / 60 + float(s) / 3600) / 24 * 360
		return deg

	def dec2deg(self):
		d, m, s = self.__Dec.split(' ')
		deg = int(d) + int(m) / 60 + float(s) / 3600
		return deg


##########    格式转换工具    ##########
gdzh = 360 / 365.25  # 每古度的今度数。今度÷该值＝古度，古度×该值＝今度
rad = 180 / math.pi  # 每弧度的度数。  弧度×该值＝度数，度数÷该值＝弧度


def deg2dms(deg):  # 度数转度分秒，仅用于输出表示
	a = deg - int(deg)  # 获取度数的小数部分
	b = a * 3600  # 小数部分全部转为角秒表示
	m = int(b / 60)  # 整数部分为分
	s = b % 60  # 余数部分为秒
	return "{:02d}°{:02d}′{:05.2f}″".format(int(deg), m, round(s, 2))


def deg2hms(deg):  # 度数转时分秒，仅用于赤经输出表示
	a = deg / 15 - int(deg / 15)
	b = a * 3600
	m = int(b / 60)
	s = b % 60
	return "{:02d}h {:02d}m {:05.2f}s ".format(int(deg / 15), m, round(s, 2))


##########    视位置修正项计算函数   ##########
def ProperMotion(star, T):  # 计算自行
	star.ra += star.pm_ra / 1000 / 3600 * T * 100 / rad
	star.dec += star.pm_dec / 1000 / 3600 * T * 100 / rad
	return star.ra, star.dec  # Date点的平坐标


def Precession(star, T):  # 历元J2000.0至T时刻的岁差
	zeta = (2306.2181 * T + 0.30188 * T ** 2 + 0.017998 * T ** 3) / 3600 / rad
	z = (2306.2181 * T + 1.09468 * T ** 2 + 0.018203 * T ** 3) / 3600 / rad
	theta = (2004.3109 * T - 0.42665 * T ** 2 - 0.041833 * T ** 3) / 3600 / rad
	A = math.cos(star.dec) * math.sin(star.ra + zeta)
	B = math.cos(theta) * math.cos(star.dec) * math.cos(star.ra + zeta) - math.sin(theta) * math.sin(star.dec)
	C = math.sin(theta) * math.cos(star.dec) * math.cos(star.ra + zeta) + math.cos(theta) * math.sin(star.dec)
	star.ra = math.atan2(A, B) + z
	star.dec = math.asin(C)
	return star.ra, star.dec  # Date点的平坐标


def Nutation(star, T):  # T时刻的章动（IAU 2000B）
	D = ((1072260.70369 + 1602961601.2090 * T - 6.3706 * pow(T, 2) + 0.006593 * pow(T, 3) - 0.00003169 * pow(T,4)) / 3600.0) / rad
	L = ((485868.249036 + 1717915923.2178 * T + 31.8792 * pow(T, 2) + 0.051635 * pow(T, 3) - 0.00024470 * pow(T,4)) / 3600.0) / rad
	Lp = ((1287104.79305 + 129596581.0481 * T - 0.5532 * pow(T, 2) + 0.000136 * pow(T, 3) - 0.00001149 * pow(T,4)) / 3600.0) / rad
	F = ((335779.526232 + 1739527262.8478 * T - 12.7512 * pow(T, 2) - 0.001037 * pow(T, 3) + 0.00000417 * pow(T,4)) / 3600.0) / rad
	Om = ((450160.398036 - 6962890.5431 * T + 7.4722 * pow(T, 2) + 0.007702 * pow(T, 3) - 0.00005939 * pow(T,4)) / 3600.0) / rad
	# 章动计算的周期序列
	cl = [0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, -1, -1, 1, -1, -1, 1, -2, 0, 0, 0, -2, 2, 1, -1, 2, 0, 0, -1, 0, 0, 1, 0, -1, 0, 1, -2, 0, 0, 0, 0, 1, 2, -2, 2, 0, 0, -1, 2, 1, 0, 1, -2, 3, 0, 1, 0, -1, -1, 0, -2, 1, 2, -1, 1, 1, -1, 1, -1, 0, -1, -1, 0, 1, -2, -1, 1]
	clp = [0, 0, 0, 0, 1, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 1, 0, 2, 0, 0, -1, 0, 2, 0, 0, 1, 0, -1, 0, 0, 0, 0, 0, -1, 0, -1, 0, 0, 1, -1, 0, 0, -1, -1, 0, -1, 0, -1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, -2, 0, 0, 0, 1]
	cf = [0, 2, 2, 0, 0, 2, 0, 2, 2, 2, 2, 2, 0, 0, 0, 2, 2, 2, 0, 2, 2, 0, 2, 2, 2, 0, 2, 0, 0, 2, -2, 0, 0, 2, 0, 2, 2, 2, 2, 2, 0, 2, 2, 0, 2, 2, 0, 0, 0, 0, 2, 0, 2, 2, 0, 2, 0, 2, 2, 2, 0, 2, 0, 0, 0, 2, 2, 0, 0, 2, 2, 0, 2, 2, 2, 0, 2]
	cd = [0, -2, 0, 0, 0, -2, 0, 0, 0, -2, -2, 0, 2, 0, 0, 2, 0, 0, 2, 2, -2, 2, 0, -2, 0, 0, 0, 0, 2, -2, 2, -2, 0, 2, 0, 2, 0, 0, 2, 0, 2, -2, -2, 2, 0, -2, -2, 2, -2, 2, -2, 0, 0, 0, 2, 0, 1, 2, 0, 2, 0, 0, 0, 1, 0, 0, -2, 0, 1, 1, 4, 1, -2, 2, 2, 0, -2]
	co = [1, 2, 2, 2, 0, 2, 0, 1, 2, 2, 1, 2, 0, 1, 1, 2, 1, 1, 0, 2, 2, 0, 2, 2, 1, 0, 0, 1, 1, 2, 0, 1, 1, 1, 0, 2, 0, 2, 1, 2, 1, 1, 2, 1, 1, 1, 1, 0, 1, 0, 1, 0, 2, 2, 0, 2, 0, 2, 0, 2, 1, 2, 1, 0, 0, 0, 1, 2, 0, 2, 2, 1, 1, 1, 2, 2, 2]
	pa = [-172064161, -13170906, -2276413, 2074554, 1475877, -516821, 711159, -387298, -301461, 215829, 128227, 123457, 156994, 63110, -57976, -59641, -51613, 45893, 63384, -38571, 32481, -47722, -31046, 28593, 20441, 29243, 25887, -14053, 15164, -15794, 21783, -12873, -12654, -10204, 16707, -7691, -11024, 7566, -6637, -7141, -6302, 5800, 6443, -5774, -5350, -4752, -4940, 7350, 4065, 6579, 3579, 4725, -3075, -2904, 4348, -2878, -4230, -2819, -4056, -2647, -2294, 2481, 2179, 3276, -3389, 3339, -1987, -1981, 4026, 1660, -1521, 1314, -1283, -1331, 1383, 1405, 1290]
	pb = [-174666, -1675, -234, 207, -3633, 1226, 73, -367, -36, -494, 137, 11, 10, 63, -63, -11, -42, 50, 11, -1, 0, 0, -1, 0, 21, 0, 0, -25, 10, 72, 0, -10, 11, 0, -85, 0, 0, -21, -11, 21, -11, 10, 0, -11, 0, -11, -11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
	pc = [33386, -13696, 2796, -698, 11817, -524, -872, 380, 816, 111, 181, 19, -168, 27, -189, 149, 129, 31, -150, 158, 0, -18, 131, -1, 10, -74, -66, 79, 11, -16, 13, -37, 63, 25, -10, 44, -14, -11, 25, 8, 2, 2, -7, -15, 21, -3, -21, -8, 6, -24, 5, -6, -2, 15, -10, 8, 5, 7, 5, 11, -10, -7, -2, 1, 5, -13, -6, 0, -353, -5, 9, 0, 0, 8, -2, 4, 0]
	ea = [92052331, 5730336, 978459, -897492, 73871, 224386, -6750, 200728, 129025, -95929, -68982, -53311, -1235, -33228, 31429, 25543, 26366, -24236, -1220, 16452, -13870, 477, 13238, -12338, -10758, -609, -550, 8551, -8001, 6850, -167, 6953, 6415, 5222, 168, 3268, 104, -3250, 3353, 3070, 3272, -3045, -2768, 3041, 2695, 2719, 2720, -51, -2206, -199, -1900, -41, 1313, 1233, -81, 1232, -20, 1207, 40, 1129, 1266, -1062, -1129, -9, 35, -107, 1073, 854, -553, -710, 647, -700, 672, 663, -594, -610, -556]
	eb = [9086, -3015, -485, 470, -184, -677, 0, 18, -63, 299, -9, 32, 0, 0, 0, -11, 0, -10, 0, -11, 0, 0, -11, 10, 0, 0, 0, -2, 0, -42, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
	ec = [15377, -4587, 1374, -291, -1924, -174, 358, 318, 367, 132, 39, -4, 82, -9, -75, 66, 78, 20, 29, 68, 0, -25, 59, -3, -3, 13, 11, -45, -1, -5, 13, -14, 26, 15, 10, 19, 2, -5, 14, 4, 4, -1, -4, -5, 12, -3, -9, 4, 1, 2, 1, 3, -1, 7, 2, 4, -2, 3, -2, 5, -4, -3, -2, 0, -2, 1, -2, 0, -139, -2, 4, 0, 0, 4, -2, 2, 0]
	dp = 0
	de = 0
	for i in range(len(pa)):
		arg = cl[i] * L + clp[i] * Lp + cf[i] * F + cd[i] * D + co[i] * Om
		dp += (pa[i] + pb[i] * T) * math.sin(arg) + pc[i] * math.cos(arg)  # 黄经章动
		de += (ea[i] + eb[i] * T) * math.cos(arg) + ec[i] * math.sin(arg)  # 交角章动
	dp /= 10000000  # Δψ 单位：角秒
	de /= 10000000  # Δε 单位：角秒
	# 黄赤交角修正
	U = T / 100
	e0 = 23 * 3600 + 26 * 60 + 21.448 - 4680.93 * U - 1.55 * pow(U, 2) + 1999.25 * pow(U, 3) - 51.38 * pow(U,4) - 249.67 * pow(U, 5) - 39.05 * pow(U, 6) + 7.12 * pow(U, 7) + 27.87 * pow(U, 8) + 5.79 * pow(U, 9) + 2.45 * pow(U, 10) # 平黄赤交角，单位角秒
	e = (e0 + de) / 3600 / rad
	# 转为赤经、赤纬章动（不适用于赤纬近于±90°的计算，宜由黄道坐标转为赤道坐标，但对二十八宿和日月五星都是有效的）
	nut_ra = (math.cos(e) + math.sin(e) * math.sin(star.ra) * math.tan(star.dec)) * dp - math.cos(star.ra) * math.tan(star.dec) * de
	nut_dec = math.sin(e) * math.cos(star.ra) * dp + math.sin(star.ra) * de
	return nut_ra, nut_dec, e # 赤经赤纬章动（单位：角秒），e单位弧度


def Aberration(eps, star, T):  # 赤经赤纬的光行差（黄赤交角ε、恒星数据项、儒略世纪数）
	JD = T * 36525 + 2451545.0
	date = ephem.Date(JD - 2415020)
    s = ephem.Sun(date)
	L = (s.hlon / ephem.degree + 180) / rad  # Date点的太阳地心坐标几何黄经（单位：弧度）
	e = (0.016708617 - 0.000042037 * T - 0.0000001236 * pow(T, 2))  # 地球轨道离心率
	pi = (102.93735 + 1.71953 * T + 0.00046 * T ** 2) / rad  # 轨道近日点经度（单位：弧度）
	k = 20.49552  # 光行差常数
	abr_ra = -k * (math.cos(star.ra) * math.cos(L) * math.cos(eps) + math.sin(star.ra) * math.sin(L)) / math.cos(
		star.dec) + e * k * (math.cos(star.ra) * math.cos(pi) * math.cos(eps) + math.sin(star.ra) * math.sin(
		pi)) / math.cos(star.dec)
	abr_dec = -k * (math.cos(L) * math.cos(eps) * (
				math.tan(eps) * math.cos(star.dec) - math.sin(star.ra) * math.sin(star.dec)) + math.cos(
		star.ra) * math.sin(star.dec) * math.sin(L)) + e * k * (math.cos(pi) * math.cos(eps) * (
				math.tan(eps) * math.cos(star.dec) - math.sin(star.ra) * math.sin(star.dec)) + math.cos(
		star.ra) * math.sin(star.dec) * math.sin(pi))
	return abr_ra, abr_dec  # 光行差修正值（单位：角秒）


##########    主函数    ##########
def ApparentPositon(star, T):  # 由J2000.0恒星位置计算任意时刻视位置
	# 获得度数表示的J2000.0恒星位置
	star.ra = star.ra2deg() / rad
	star.dec = star.dec2deg() / rad
	# 获得T时刻平位置
	ProperMotion(star, T)  # 计算自行
	Precession(star, T)  # 计算岁差
	# 获得T时刻视位置
	nut_ra, nut_dec, eps = Nutation(star, T)  # 计算章动和黄赤交角
	abr_ra, abr_dec = Aberration(eps, star, T)  # 计算光行差
	star.ra = (star.ra * rad + nut_ra / 3600 + abr_ra / 3600) % 360  # 修正章动和光行差
	star.dec = (star.dec * rad + nut_dec / 3600 + abr_dec / 3600)
	return star.ra, star.dec  # 度数表示的视赤经、视赤纬


# 数据来源：依巴谷星表
# 数据格式（中文名，英文名，赤经h m s，赤纬d m s，赤经自行，赤纬自行
jiao = Star('角', 'α Vir', '13 25 11.5793', "-11 9 11.5793", -42.35, -30.67)  # Spica
kang = Star('亢', 'κ Vir', '14 12 53.74109', "-10 16 26.5579", 7.25, 139.88)
di = Star('氐', 'α Lib', '14 50 52.77724', "-16 2 29.7970", -105.68, -68.4)  # Zubenelgenubi
fang = Star('房', 'π Sco', '15 58 51.12066', "-26 6 50.5538", -11.42, -26.83)
xin = Star('心', 'σ Sco', '16 21 11.32257', "-25 35 33.9090", -10.6, -16.28)  # 一为α Sco
wei1 = Star('尾', 'μ¹ Sco', '16 51 52.23895', "-38 2 50.3764", -10.58, -22.06)
ji = Star('箕', 'γ Sgr', '18 5 48.52457', "-30 25 25.1406", -53.92, -180.9)
dou = Star('斗', 'φ Sgr', '18 45 39.35297', "-26 59 26.8051", 50.61, 1.22)
niu = Star('牛', 'β Cap', '20 21 0.64616', "-14 46 53.0436", 44.92, 7.38)
nv = Star('女', 'ε Aqr', '20 47 40.53251', "-9 29 44.4835", 33.98, -34.77)
xu = Star('虚', 'β Aqr', '21 31 33.52071', "-5 34 16.1602", 18.77, -8.21)
wei2 = Star('危', 'α Aqr', '22 5 47.03593', "-0 19 11.4568", 18.25, -9.39)  # Sadalmelik
shi = Star('室', 'α Peg', '23 4 45.61693', "15 12 19.3231", 61.1, -42.56)  # Markab
bi = Star('壁', 'γ Peg', '0 13 14.15003', "15 11 1.0180", 1.98, -9.28)  # Algenib
kui = Star('奎', 'ζ And', '0 47 20.39021', "24 16 2.5563", -101.17, -81.77)
# kui = Star('奎', 'η And', '0 57 12.42782', "23 25 3.9358", -43.72, -46.06) # 宋以后改
lou = Star('娄', 'β Ari', '1 54 38.34937', "20 48 29.8794", 98.74, -110.41)
wei3 = Star('胃', '35 Ari', '2 43 27.11185', "27 42 25.7233", 2.06, -10.37)
mao = Star('昴', '17 Tau', '3 44 52.52356', "24 6 48.4142", 21.55, -44.92)  # Electra
bi3 = Star('毕', 'ε Tau', '4 28 36.93323', "19 10 49.8757", 106.19, -37.84)
zi = Star('觜', 'φ1 Ori', '5 34 49.23788', "9 29 22.5076", 0.27, -2.26)
# zi = Star('觜', 'λ Ori', '5 35 8.27781', "9 56 2.9869",  -0.34,  -2.94) # 宋以后改
shen = Star('参', 'δ Ori', '5 32 0.39972', "-0 17 56.7364", 0.64, -0.69)  # 一为ζ Ori # Mintaka
jing = Star('井', 'μ Gem', '6 22 57.59125', "22 30 49.8606", 56.39, -110.03)
gui = Star('鬼', 'θ Cnc', '8 31 35.72962', "18 5 39.9160", -60.85, -56.14)
liu = Star('柳', 'δ Hya', '8 37 39.40742', "5 42 13.6749", -70.19, -7.9)
xing = Star('星', 'α Hya', '9 27 35.25169', "-8 39 31.2591", -15.23, 34.37)  # Alphard
zhang = Star('张', 'υ¹ Hya', '9 51 28.68244', "-14 50 47.5798", 18.88, -21.85)
yi = Star('翼', 'α Cra', '10 59 46.74887', "-18 17 56.7503", -462.26, 129.49)
zhen = Star('轸', 'γ Crv', '12 15 48.46784', "-17 32 31.1409", -158.61, 21.86)  # Gienah
xin2 = Star('心宿二', 'α Sco', '16 29 24.46759', "-26 25 55.0055", -12.11, -23.3)  # Antares
esbx = [jiao, kang, di, fang, xin, wei1, ji, dou, niu, nv, xu, wei2, shi, bi, kui, lou, wei3, mao, bi3, zi, shen, jing,
        gui, liu, xing, zhang, yi, zhen]  # 二十八宿
szzx = [xing, xin2, xu, mao]  # 四仲中星

# 计算示例
date = '2000/1/1 00:00:00'  # UT时间，近似用于TD
JD = ephem.julian_date(date)
t = (JD - 2451545) / 36525
for star in esbx:
	ApparentPositon(star, t)
	RA = deg2hms(star.ra)
	DEC = deg2dms(star.dec)
	print(star.name, RA, DEC)

参考资料

基本算法：Jean Meeus《Astronomical Algorithms》2nd

依巴谷星表数据引用：https://www.heavens-above.com、https://www.wikipedia.org

二十八宿距星：陈遵妫《中国天文学史》第二卷