#由方位角高度角计算赤经赤纬
from ephem import *
import math
#(1)构造观测者位置变量
wpf=Observer()
wpf.lat='34.179569'#纬度
wpf.lon='115.93635'#经度
wpf.elevation=40#海拔
wpf.pressure=0#压强
wpf.date=now()#取当前时间
#wpf.date="2020/4/10 10:44:53"#时间根据需要修改
p=Venus()#金星
p.compute(wpf)
#世界时转为儒略日
def utc2jd(t):s=(Date(t)).tuple()#世界时的y=s[0]#年m=s[1]#月d=s[2]#日h=s[3]#时dh=h+s[4]/60.0+s[5]/3600.0if m <= 2:y = y - 1m = m + 12day_num=int(365.25*y) + int(30.6001*(m + 1))+ d + int(dh / 24.0 + 1720981.5)hour_num=(h + 12) % 24#儒略日从12点取jd=day_num+hour_num/24+s[4]/1440+s[5]/86400return jd
#由世界时求恒星时，t是格林尼治时间，返回值是角时
def sidereal_t(t):jd=utc2jd(t)mjd=jd-2451545.0#约化儒略日Tu=mjd/36525s0=6.697374558+2400.05134*Tu+0.000025862333*Tu*Tu-0.000000001722222*Tu*Tu*Tu#格林尼治时0时的恒星时s=(Date(t)).tuple()M=s[3]+s[4]/60.0+s[5]/3600.0#格林尼治时的小时数s=s0+M*366.2422/365.2422#格林尼治恒星时时加上地方经度（单位化为小时）等于地方恒星时s=s%24return s
#给出观测者，天体，时间信息，求天体方位角高度角
def ra_dec(obs,a,h):s=sidereal_t(obs.date)#格里尼治恒星时ls=wpf.sidereal_time()ls=s+obs.lon*180/math.pi/15#地方恒星时ls=ls*15*math.pi/180#角时转为弧度ls1=wpf.sidereal_time()sint=-math.cos(h)*math.sin(a)cost=math.cos(obs.lat)*math.sin(h)-math.cos(h)*math.sin(obs.lat)*math.cos(a)t=math.atan2(sint,cost)if t<0:t=t+2*math.pira=ls-tra1=ls1-tdec=math.asin(math.sin(obs.lat)*math.sin(h)+math.cos(obs.lat)*math.cos(h)*math.cos(a))return ra,ra1,dec
#A,h= az_alt(wpf,s,t)
a=p.az
h=p.alt
r,r1,d=ra_dec(wpf,a,h)
print(degrees(p.ra),degrees(p.dec))
print(degrees(r1),degrees(d))
print(degrees(r),degrees(d))