python太阳代码_利⽤python计算太阳天顶⾓、⽅位⾓、⾼度
1、编程语⾔改为Python
2、输⼊和输出⽂件格式均改为excel
3、以为太阳⾼度⾓与天顶⾓互余,这次也加进去了。
欢迎讨论,下⾯附上主要代码:
table = data.sheets()[0] # 打开第⼀张表
nrows = ws # 获取表的⾏数
ncols = ls
l_values(0)[1:]
l_values(1)[1:]
l_values(2)[1:]
l_values(3)[1:]
l_values(4)[1:]
l_values(5)[1:]
l_values(6)[1:]
l_values(7)[1:]
lat =l_values(8)[1:]
TimeZone =l_values(9)[1:]
wb = py(data)
ws = wb.get_sheet(0)
ws.write(0, 11, 'Day of Year')
ws.write(0, 12, 'Local Time') #实际是gtdt
ws.write(0, 13, 'Sun Angle')#(sitar)
ws.write(0, 14, 'Declination Angle')
ws.write(0, 15, 'Equation of Time')
style = xlwt.easyxf('pattern: pattern solid, fore_color yellow;')
ws.write(0, 16, 'ZenithAngle(deg)',style)
ws.write(0, 17, 'HeightAngle(deg)',style)
ws.write(0, 18, 'AzimuthAngle(deg)',style)
for n in range(1,nrows):
m=n-1
#儒略⽇ Julian day(由通⽤时转换到儒略⽇)
JD0 = int(365.25*(year[m]-1))+int(30.6001*(1+13))+1+hour[m]/24+1720981.5
if month[m]<=2:
JD2 = int(365.25*(year[m]-1))+int(30.6001*(month[m]+13))+day[m]+hour[m]/24+1720981.5
else:
JD2 = int(365.25*year[m])+int(30.6001*(month[m]+1))+day[m]+hour[m]/24+1720981.5
#年积⽇ Day of year
DOY = JD2-JD0+1
N0 sitar=θ
N0 = 79.6764 + 0.2422*(year[m]-1985) - int((year[m]-1985)/4.0)
sitar = 2*math.pi*(DOY-N0)/365.2422
ED1 = 0.3723 + 23.2567*math.sin(sitar) + 0.1149*math.sin(2*sitar) - 0.1712*math.sin(3*sitar)- 0.s(sitar) + 0.s(2*sitar) + 0.s(3*sitar)
ED = ED1*math.pi/180 #ED本⾝有符号
if lon[m] >= 0:
if TimeZone == -13:
dLon = lon[m] - (math.floor((lon[m]*10-75)/150)+1)*15.0
else:
dLon = lon[m] - TimeZone[m]*15.0 #地球上某⼀点与其所在时区中⼼的经度差
else:
if TimeZone[m] == -13:
dLon = (math.floor((lon[m]*10-75)/150)+1)*15.0- lon[m]
else:
dLon = TimeZone[m]*15.0- lon[m]
#时差
Et = 0.0028 - 1.9857*math.sin(sitar) + 9.9059*math.sin(2*sitar) - 7.s(sitar)- 0.s(2*sitar)
gtdt1 = hour[m] + min[m]/60.0 + sec[m]/3600.0 + dLon/15 #地⽅时
gtdt = gtdt1 + Et/60.0
dTimeAngle1 = 15.0*(gtdt-12)
dTimeAngle = dTimeAngle1*math.pi/180
latitudeArc = lat[m]*math.pi/180
⾼度⾓计算公式
HeightAngleArc =
math.asin(math.sin(latitudeArc)*math.sin(ED)+s(latitudeArc)*s(ED)*s(dTimeAngle))
python代码转换CosAzimuthAngle = (math.sin(HeightAngleArc)*math.sin(latitudeArc)-
math.sin(ED))/s(HeightAngleArc)/s(latitudeArc)
AzimuthAngleArc = math.acos(CosAzimuthAngle)
HeightAngle = HeightAngleArc*180/math.pi
ZenithAngle = 90-HeightAngle
AzimuthAngle1 = AzimuthAngleArc *180/math.pi
if dTimeAngle < 0:
AzimuthAngle = 180 - AzimuthAngle1
else:
AzimuthAngle = 180 + AzimuthAngle1
print('站位:'+station[m]+' 太阳天顶⾓(deg):%f ⾼度⾓(deg):%f ⽅位⾓(deg):%f ' % (ZenithAngle,HeightAngle,AzimuthAngle))
附上运⾏结果:
图: 天顶⾓、⽅位⾓、⾼度⾓运⾏结果
最后,感谢汪⾃军博⼠的思路和中国⽓象科学研究院王炳忠研究员编写的《太阳辐射计算讲座》。
2017年4⽉27⽇

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。