太陽高度の取得

2022/12/15

Python3.6.7

from math import sin,cos,asin,pi,degrees,radians

lat = radians(35.6544) #緯度
lon = radians(139.7447) #経度

days = 348 #12/15
hour = 12

# 今の時刻にする場合
#from datetime import datetime
#dt1 = datetime.now()
#dt2 = datetime(dt1.year,1,1)
#dt3 = dt1 - dt2
#days = dt3.days
#hour = dt1.hour

# 太陽赤緯
w = 2*pi/365
J = days + 0.5
sun = (  0.33281 - 22.984*cos(w*J) - 0.34990*cos(2*w*J)
       - 0.13980*cos(3*w*J) + 3.7872*sin(w*J)
       + 0.03250*sin(2*w*J) + 0.07187*sin(3*w*J) )
sun = radians(sun)
#print('sun',degrees(sun))

# 均時差
e = ( 0.0072*cos(w*J) - 0.0528*cos(2*w*J) - 0.0012*cos(3*w*J)
    - 0.1229*sin(w*J) - 0.1565*sin(2*w*J) - 0.0041*sin(3*w*J))
#print('e',e)

# 時角
Ts = hour
T = Ts + (degrees(lon) - 135)/15 + e
t = 15*T - 180
t = radians(t)
#print('t',degrees(t))

#太陽高度
sinh = sin(lat)*sin(sun) + cos(lat)*cos(sun)*cos(t)
h = asin(sinh)
print('h', degrees(h))

参照先：■太陽の高度と方位角（JavaScriptによる簡易版）
http://k-ichikawa.blog.enjoy.jp/etc/HP/js/sunShineAngle/ssa.html

違うサイトの掲載式も試してみました（微妙に違いますが、ほぼ同じです）。

from math import sin,cos,asin,pi,degrees,radians

lat = radians(35.6544) #緯度
lon = radians(139.7447) #経度

days = 348 #12/15
hour = 12

po = 2*pi*(days-1)/365
sun = (0.006918-0.399912*cos(po)+0.070257*sin(po)-0.006758*cos(2*po)+
       0.000907*sin(2*po)-0.002697*cos(3*po)+0.001480*sin(3*po) )
#print('sun',degrees(sun))

# 均時差
e = ( 0.000075+0.001868*cos(po)-0.032077*sin(po)
      - 0.014615*cos(2*po)-0.040849*sin(2*po) )
#print('e',e)

JST = hour
t = (JST-12)*pi/12+(lon-radians(135))+e
#print('t',degrees(t))
h = asin(sin(lat)*sin(sun)+cos(lat)*cos(sun)*cos(t))
print('h', degrees(h))

参照先:太陽方位、高度、大気外日射量の計算
http://www.es.ris.ac.jp/~nakagawa/met_cal/solar.html

答え合わせ：太陽高度（一日の変化）
https://keisan.casio.jp/exec/system/1185781259

解説：太陽高度・太陽方位角の計算式・符号・求め方について
https://environmental-engineering.work/archives/219