太阳方位角和高度角计算

太阳高度角和方位角的计算

代码思想

  1. 分析公式中需要的参数，分别为地理纬度、太阳赤纬、时角
  2. 先计算年积日，根据年积日计算太阳赤纬
  3. 将地理纬度换算成弧度
  4. 计算时角
  5. 根据上述参数传入方位角和高度角函数计算高度角和方位角
  6. 定义所要计算的时间地区
  7. 对结果进行输出

代码示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
import math

# 年积日的计算
def yearToYear(year, month, day, hour):
    # 儒略日(由常规通用时转换到儒略日)
    JD0 = int(365.25*(year-1))+int(30.6001*(1+13))+1+hour/24+1720981.5
    if month <= 2:
        JD2 = int(365.25*(year-1))+int(30.6001 *
                                       (month+13))+day+hour/24+1720981.5
    else:
        JD2 = int(365.25*year)+int(30.6001*(month+1))+day+hour/24+1720981.5
    # 年积日DOY
    DOY = JD2-JD0+1
    return DOY

# 太阳赤纬的计算
def EDCompute(DOY, year):
    N0 = 79.6764 + 0.2422*(year-1985) - int((year-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.758*math.cos(sitar) + 0.3656*math.cos(2*sitar) + 0.0201*math.cos(3*sitar)
    ED = ED1*math.pi/180  # ED本身有符号
    return ED

# 地理纬度的计算
def latCompute(lat):
    latitudeArc = lat*math.pi/180
    return latitudeArc

# 时角的计算
def dTimeCompute(lon, TimeZone, DOY, year):
    N0 = 79.6764 + 0.2422*(year-1985) - int((year-1985)/4.0)
    sitar = 2*math.pi*(DOY-N0)/365.2422
    if lon >= 0:
        if TimeZone == -13:
            dLon = lon - (math.floor((lon*10-75)/150)+1)*15.0
        else:
            dLon = lon - TimeZone*15.0  # 地球上某一点与其所在时区中心的经度差
    else:
        if TimeZone == -13:
            dLon = (math.floor((lon*10-75)/150)+1)*15.0 - lon
        else:
            dLon = TimeZone*15.0 - lon
    # 时差
    Et = 0.0028 - 1.9857*math.sin(sitar) + 9.9059*math.sin(2*sitar) - \
        7.0924*math.cos(sitar) - 0.6882*math.cos(2*sitar)
    gtdt1 = hour + min/60.0 + sec/3600.0 + dLon/15  # 地方时
    gtdt = gtdt1 + Et/60.0
    dTimeAngle1 = 15.0*(gtdt-12)
    dTimeAngle = dTimeAngle1*math.pi/180
    return dTimeAngle

# 高度角的计算
def heightAngle(latitudeArc, ED, dTimeAngle):
    HeightAngleArc = math.asin(math.sin(
        latitudeArc)*math.sin(ED)+math.cos(latitudeArc)*math.cos(ED)*math.cos(dTimeAngle))
    return HeightAngleArc

# 方位角的计算
def azimuthAngle(HeightAngleArc, latitudeArc, ED):
    CosAzimuthAngle = (math.sin(HeightAngleArc)*math.sin(latitudeArc) -
                       math.sin(ED))/math.cos(HeightAngleArc)/math.cos(latitudeArc)
    AzimuthAngleArc = math.acos(CosAzimuthAngle)
    AzimuthAngle1 = AzimuthAngleArc * 180/math.pi
    if dTimeAngle < 0:
        AzimuthAngle = 180 - AzimuthAngle1
    else:
        AzimuthAngle = 180 + AzimuthAngle1
    return AzimuthAngle

if __name__ == '__main__':
    # 参数设置
    year = 2022  # 年份
    month = 7  # 月份
    day = 31  # 日期
    hour = 9  # 时间
    min = 00  # 分钟
    sec = 00  # 秒
    lon = 120  # 经度
    lat = 30  # 纬度
    TimeZone = 8  # 时区
    DOY = yearToYear(year, month, day, hour)  # 年积日计算
    ED = EDCompute(DOY, year)  # 太阳赤纬计算
    latitudeArc = latCompute(lat)  # 地理纬度计算
    dTimeAngle = dTimeCompute(lon, TimeZone, DOY, year)  # 时角计算
    HeightAngle = heightAngle(latitudeArc, ED, dTimeAngle)*180/math.pi  # 高度角计算
    AzimuthAngle = azimuthAngle(heightAngle(
        latitudeArc, ED, dTimeAngle), latitudeArc, ED)  # 方位角计算
    print('%s年%s月%s日%s时%s分%s秒，地理位置位于经度%s，纬度%s的高度角(deg)为：%f，方位角(deg)为：%f ' %
          (year, month, day, hour, min, sec, lon, lat, HeightAngle, AzimuthAngle))

运行结果

image-20221107202928594