接上一篇用ephem库计算24节气时间 ，用ephem库可以求得任意时间太阳的视黄经，因而可以求得24节气的时间；同理我们也可以求得任意时间几大行星的视黄经，有了和时间对应的行星黄经值，我们就可以计算一些和行星相关的天象的发生时间，像地内行星(水星、金星)的上合、下合、东大距、西大距，地外行星(火星、木星、土星、天王星、海王星)的上合、大冲、东方照、西方照，月亮的满月/望月、新月/朔月、上弦月、下弦月时间；几大行星的逆行时间，我们都可以迭代计算出来，因为这些天象都是和行星与太阳的相对位置有关，而它们的相对位置可以通过黄道坐标来计算。conjunction

地内行星地外行星大距

1.首先计算一下上面提到的和行星在一个会合周期内的四个特殊位置所构成的天文现象的发生时间

一个会合周期内，从地球观察，行星或者月亮相对于太阳的几个特殊位置：

地内行星：水星、金星

上合：行星和太阳黄经差绝对值减小到0°，即二者黄经相等，行星从太阳西边跑到太阳东边

东大距：行星在太阳东边，昏星，和太阳夹角开始从变大改为变小，即二者夹角达到极大

下合：行星和太阳黄经差绝对值减小到0°，即二者黄经相等，行星从太阳东边跑到太阳西边

西大距：行星在太阳西边，晨星，和太阳夹角开始从变大改为变小，即二者夹角达到极大

变化顺序：上合-东大距-下合-西大距-上合

地外行星：火星、木星、土星、天王星、海王星

上合：行星和太阳黄经差绝对值减小到0°，即二者黄经相等，行星从太阳西边跑到太阳东边

东方照：行星和太阳黄经差绝对值为90°，行星在太阳东边

大冲：行星和太阳黄经差绝对值为180°，行星从太阳东边跑到太阳西边

西方照：行星和太阳黄经差绝对值为90°，行星在太阳西边

变化顺序：上合-西方照-大冲-东方照-上合

月亮：

上弦月：月亮和太阳黄经差绝对值90°，月亮在太阳东边

满月/望：月亮和太阳黄经差绝对值为180°，月亮从太阳东边跑到太阳西边

下弦月：月亮和太阳黄经差绝对值90°，月亮在太阳西边

无月/朔：月亮和太阳黄经差绝对值为0°，月亮从太阳东边跑到太阳西边

变化顺序：上弦月-满月-下弦月-无月-上弦月

地外行星的位置变化顺序看上去是在自东向西转，其实是因为相对地球而言被反超。

所以下面代码主要有判断行星在太阳东或西的函数、上合/下合/大冲函数、大距函数、黄经差函数、夹角函数、迭代函数。

from ephem import *

import math

import datetime

event=[]#不同天体不同状态对应的天象名称

event.append({1:"西大距",2:"上　合",3:"东大距",4:"下　合"})#地内行星

event.append({1:"西方照",4:"上　合",2:"冲　日",3:"东方照"})#地外行星，是太阳在追行星，即太阳的黄经慢慢靠近行星

event.append({1:"下弦月",2:"新　月",3:"满　月",4:"下弦月"})#月亮

#设置数字对应的天体名称

body={0:"月亮",1:"水星",2:"金星",4:"火星",5:"木星",6:"土星",7:"天王星",8:"海王星"}

#迭代求天象时间

def iteration(jd,n,sta):#jd：要求的开始时间，n：上面不同天体的代号，sta：不同的状态函数

s1=sta(jd,n)#初始状态(状态可以是处于什么位置，或者是否处于逆行)

s0=s1

dt=1.0#初始时间改变量设为1天

while True:

jd+=dt

s=sta(jd,n)

if s0!=s:

s0=s

dt=-dt/2#使时间改变量折半减小

if abs(dt)<0.0000001 and s==s1:#s==s1是为了让求得的时间在天象发生之前

break

return jd

#计算黄经

def ecliptic_coor(jd_utc,n):

if n==0:

p=Moon(jd_utc)#月球

elif n==1:

p=Mercury(jd_utc)#水星

elif n==2:

p=Venus(jd_utc)#金星

elif n==4:

p=Mars(jd_utc)#火星

elif n==5:

p=Jupiter(jd_utc)#木星

elif n==6:

p=Saturn(jd_utc)#土星

elif n==7:

p=Uranus(jd_utc)#天王星

elif n==8:

p=Neptune(jd_utc)#海王星

else:

p=Sun(jd_utc)#太阳

eq=Equatorial(p.ra,p.dec,epoch=jd_utc)#求天体的视赤经视赤纬(epoch设为所求时间就是视赤经视赤纬)

ec=Ecliptic(eq)#赤经赤纬转到黄经黄纬

return ec.lon,ec.lat#返回黄经

#return eq.ra,eq.dec

#给定两个黄经a1和a2，判断a2在a1的哪一侧，东(设为0)或者西(设为1)

def east_west_angle(a1,a2):

if abs(a1-a2)

if a2

return 1

else:

return 0

else:

if a2

return 0

else:

return 1

#下面四个函数一般设m为太阳，也就是求相对于太阳的参数

#求两个天体之间的黄经差

def ecliptic_angle(jd,m,n):

a1=ecliptic_coor(jd,m)

a2=ecliptic_coor(jd,n)

if abs(a1[0]-a2[0])

return abs(a1[0]-a2[0])

else:

return 2*math.pi-abs(a1[0]-a2[0])

#计算两个天体之间的夹角

def included_angle(jd_utc,m,n):

a1=ecliptic_coor(jd_utc,m)

a2=ecliptic_coor(jd_utc,n)

a=math.acos(math.sin(a1[1])*math.sin(a2[1])+math.cos(a1[1])*math.cos(a2[1])*math.cos(a2[0]-a1[0]))

return a

#判断两个天体之间夹角在增大还是减小

def elongation(jd_utc,n):#大距判断函数

a1=included_angle(jd_utc,-1,n)

a2=included_angle(jd_utc+1/86400,-1,n)

if a1

return 1

else:#夹角减小

return 0

#两个天体之间黄经之差是否大于90度

def quadrature(jd_utc,n):#方照判断函数

a1=ecliptic_coor(jd_utc,-1)

a2=ecliptic_coor(jd_utc,n)

d=abs(a1[0]-a2[0])#黄经之差

if d>math.pi:#夹角要小于180°

d=2*math.pi-d

if d>math.pi/2:#地外行星或者月亮和太阳的黄经之差绝对值大于90度

return 1

else:#地外行星或者月亮和太阳的黄经之差绝对值小于90度

return 0

#根据黄经判断天体n在太阳的哪一侧，东(设为0)或者西(设为1)

def east_west(jd_utc,n):#上合、下合、大冲判断函数

a1=ecliptic_coor(jd_utc,-1)

a2=ecliptic_coor(jd_utc,n)

return east_west_angle(a1[0],a2[0])

#判断并输出计算的天象名称

def print_time(jd,n):

if n==1 or n==2:

i=0

elif n==0:

i=2

else:

i=1

e=ecliptic_angle(jd,-1,n)*180/math.pi#判断黄经差值

#这里是为了把太阳和行星的合/冲时的黄经差(接近0°)跟大距/方照是的黄经差区分开

if abs(180-e)<1:#为了把大冲时的黄经差(接近180°)转为0°，这样太阳和行星的合/冲时的黄经差判断标准统一为0°，1°已足够区分合/冲与大距/方

e=abs(180-e)#大冲

d=east_west(jd,n)#判断在东在西

if d==1 and e>1:

s=1

if d==1 and e<1:

s=2

if d==0 and e>1:

s=3

if d==0 and e<1:

s=4

if i==0 and (s==1 or s==3):#地内行星东大距或者西大距时输出和太阳夹角

a=included_angle(jd,-1,n)

print("{0}{1}：{2},夹角：{3}".format(body[n], event[i][s],Date(jd+1/3),a*180/math.pi))

else:

#if s==1:

print("{0}{1}：{2}".format(body[n], event[i][s],Date(jd+1/3)))

#计算天体n在jd时间时的下个大距时间

def next_elongation(jd,n):

s=elongation(jd,n)

if s==0:#跳过求夹角极小值的时间

jd=iteration(jd,n,elongation)+1

return iteration(jd,n,elongation)

#计算天体n在jd时间时的下个方照时间

def next_quadrature(jd,n):

return iteration(jd,n,quadrature)

#计算天体n在jd时间时的下个合/冲/满月/新月时间

def next_conjunction(jd,n):

return iteration(jd,n,east_west)

#求未来多次天体处于特殊几何位置的时间

def event_time(jd,n,num):#分别输入时间(ephem儒略日)，天体，要计算的接下来要发生的天象个数

jd0=jd

k=0

while True:

jd1=next_conjunction(jd0,n)#地内行星和地外行星的合/大冲统一使用next_conjunction来求

k+=1

if k==num:

print_time(jd1,n)

break

if n==1 or n==2:

jd2=next_elongation(jd0,n)#地内行星求大距

else:

jd2=next_quadrature(jd0,n)#地外行星求方照

k+=1

if jd1

t1=jd1

t2=jd2

else:

t1=jd2

t2=jd1

print_time(t1,n)

print_time(t2,n)

if k==num:

break

jd0=t2+2#加上2天以防止天体还是处在上一个状态

jd=now()

event_time(jd,2,8)

event_time(jd,4,8)

金星下　合：2020/6/4 01:43:28

金星西大距：2020/8/13 08:13:33,夹角：45.79094800265633

金星上　合：2021/3/26 14:57:32

金星东大距：2021/10/30 04:51:18,夹角：47.04505287489652

金星下　合：2022/1/9 08:47:30

金星西大距：2022/3/20 17:25:26,夹角：46.58633650874031

金星上　合：2022/10/23 05:17:06

金星东大距：2023/6/4 19:00:37,夹角：45.39928405369522

火星西方照：2020/6/7 03:10:53

火星大　冲：2020/10/14 07:25:46

火星东方照：2021/2/1 18:33:15

火星上　合：2021/10/8 12:00:51

火星西方照：2022/8/27 13:28:07

火星大　冲：2022/12/8 13:41:25

火星东方照：2023/3/17 02:09:40

火星上　合：2023/11/18 13:41:22

2.ephem库有现成的函数可以输出月亮的满月/望月、新月/朔月、上弦月、下弦月时间，上面的函数同样也可以计算，我们可以对比一下二者结果来看一下计算精度。

jd=now()

d=next_full_moon(jd)

print(Date(d+1/3))

d=next_last_quarter_moon(jd)

print(Date(d+1/3))

d=next_new_moon(jd)

print(Date(d+1/3))

d=next_first_quarter_moon(jd)

print(Date(d+1/3))

event_time(jd,0,4)

计算结果：

2020/1/11 03:21:17

2020/1/17 20:58:24

2020/1/25 05:41:59

2020/2/2 09:41:40

月亮满　月：2020/1/11 03:21:17

月亮下弦月：2020/1/17 20:58:24

月亮新　月：2020/1/25 05:41:59

月亮下弦月：2020/2/2 09:41:40

3.其次来计算一下行星逆行的时间，网上经常有人说“水逆”，水星的2020年的具体逆行时间可以计算如下：

#判断天体是否处于逆行状态,根据黄经判断天体n1秒后在之前位置的哪一侧，东(顺行，设为0)或者西(逆行，设为1)

def retrograde(jd,n):

a1=ecliptic_coor(jd,n)

a2=ecliptic_coor(jd+1/86400,n)

return east_west_angle(a1[0],a2[0])

#求下次逆行或者顺行的开始时间

def next_retrograde(jd,n):

return iteration(jd,n,retrograde)

#求未来多次逆行或者顺行开始时间

def retrograde_time(jd,n,num):

s=retrograde(jd,n)#初始状态

if s==0:

print("此刻{0}{1}在顺行".format(Date(jd+1/3),body[n]))

else:

print("此刻{0}{1}在逆行".format(Date(jd+1/3),body[n]))

jd0=jd

for i in range(num):

s=retrograde(jd0,n)

jd=next_retrograde(jd0,n)

if s==0:

print("{0}逆行开始时间：{1}".format(body[n], Date(jd+1/3)))

else:

print("{0}逆行结束时间：{1}".format(body[n], Date(jd+1/3)))

jd0=jd+2

jd=now()

retrograde_time(jd,1,10)

计算结果如下：

此刻2020/1/8 22:14:14水星在顺行

水星逆行开始时间：2020/2/17 08:53:55

水星逆行结束时间：2020/3/10 11:48:35

水星逆行开始时间：2020/6/18 12:58:49

水星逆行结束时间：2020/7/12 16:26:07

水星逆行开始时间：2020/10/14 09:04:59

水星逆行结束时间：2020/11/4 01:49:37

水星逆行开始时间：2021/1/30 23:51:40

水星逆行结束时间：2021/2/21 08:51:59

水星逆行开始时间：2021/5/30 06:33:42

水星逆行结束时间：2021/6/23 05:59:56

4.行星合/伴月时间

和上面计算行星的上合时间类似，只是把参照天体由太阳换成了月亮，判断的标准把黄经改为赤经即可，而且无需判断行星在月亮的东边还是西边。

from ephem import *

import math

import datetime

#设置数字对应的天体名称

body={0:"月亮",1:"水星",2:"金星",4:"火星",5:"木星",6:"土星",7:"天王星",8:"海王星"}

#迭代求天象时间

def iteration(jd,n,sta):#jd：要求的开始时间，n：上面不同天体的代号，sta：不同的状态函数

s1=sta(jd,n)#初始状态(状态可以是处于什么位置，或者是否处于逆行)

s0=s1

dt=1.0#初始时间改变量设为1天

while True:

jd+=dt

s=sta(jd,n)

if s0!=s:

s0=s

dt=-dt/2#使时间改变量折半减小

if abs(dt)<0.0000001 and s==s1:#s==s1是为了让求得的时间在天象发生之前

break

return jd

#计算黄经

def ecliptic_coor(jd_utc,n):

if n==0:

p=Moon(jd_utc)#月球

elif n==1:

p=Mercury(jd_utc)#水星

elif n==2:

p=Venus(jd_utc)#金星

elif n==4:

p=Mars(jd_utc)#火星

elif n==5:

p=Jupiter(jd_utc)#木星

elif n==6:

p=Saturn(jd_utc)#土星

elif n==7:

p=Uranus(jd_utc)#天王星

elif n==8:

p=Neptune(jd_utc)#海王星

else:

p=Sun(jd_utc)#太阳

return p.ra,p.dec#返回赤经赤纬

#给定两个赤经a1和a2，判断a2在a1的哪一侧，东(设为0)或者西(设为1)

def east_west_angle(a1,a2):

if abs(a1-a2)

if a2

return 1

else:

return 0

else:

if a2

return 0

else:

return 1

#求两个天体之间的黄经差(绝对值)

def ecliptic_angle(jd,m,n):

a1=ecliptic_coor(jd,m)

a2=ecliptic_coor(jd,n)

if abs(a1[0]-a2[0])

return abs(a1[0]-a2[0])

else:

return 2*math.pi-abs(a1[0]-a2[0])

#根据黄经判断天体n在天体m的哪一侧，东(设为0)或者西(设为1)

def east_west(jd_utc,n):#上合、下合、大冲判断函数

a1=ecliptic_coor(jd_utc,0)

a2=ecliptic_coor(jd_utc,n)

return east_west_angle(a1[0],a2[0])

#计算天体n在jd时间时的下个合/冲/满月/新月时间

def next_conjunction(jd,n):

return iteration(jd,n,east_west)

def included_angle(jd_utc,m,n):

a1=ecliptic_coor(jd_utc,m)

a2=ecliptic_coor(jd_utc,n)

a=(math.acos(math.sin(a1[1])*math.sin(a2[1])+math.cos(a1[1])*math.cos(a2[1])*math.cos(a2[0]-a1[0])))*180.0/math.pi

return a

#求未来多次天体处于特殊几何位置的时间

def event_time(jd,num):#分别输入时间(ephem儒略日)，要计算的接下来要发生的天象个数

for i in [2,4,5,6]:#输入要求的行星代号，可同时求多个

jd1=jd

for j in range(num):

jd1=next_conjunction(jd1,i)

e=ecliptic_angle(jd1,0,i)*180/math.pi#判断行星和月亮黄经差值,应该为0或者180

if e<1:#1°足够判断其是0还是180

a=included_angle(jd1,0,i)

print("{0}合月时间：{1},夹角：{2}".format(body[i],Date(jd1+1/3),a))

else:

jd1=next_conjunction(jd1+1,i)#证明上次计算的是行星和月亮黄经差值为180，需要再算一次，就是伴月时间

a=included_angle(jd1,0,i)

print("{0}合月时间：{1},夹角：{2}".format(body[i],Date(jd1+1/3),a))

jd1=jd1+1#防止再次计算本次合月时间，确保计算的是下次

jd=now()

event_time(jd,2)

结果：

金星合月时间：2020/4/26 23:24:41,夹角：6.053646752648658

金星合月时间：2020/5/24 10:41:39,夹角：3.696193694618995

火星合月时间：2020/4/16 12:32:42,夹角：2.0007430696592174

火星合月时间：2020/5/15 10:01:29,夹角：2.7540416415394176

木星合月时间：2020/4/15 07:04:40,夹角：1.999709418220439

木星合月时间：2020/5/12 17:40:49,夹角：2.2534055018064207

土星合月时间：2020/4/15 17:17:57,夹角：2.4626117119043833

土星合月时间：2020/5/13 02:10:35,夹角：2.675694194064673