日の出・日の入りの時刻を求める計算
日の出と日の入りの時刻を求めるコードをPythonで作成した。
計算方法は、参考文献（下の本）で紹介されている逐次近似法による反復計算。このコードを実行すると、2016年（変数yearで指定）の毎月3日（変数dayで指定）における東京での（緯度latitudeと経度longitudeで指定）日の出・日の入り時刻を算出する。綺麗に求められるので驚いた。この文献は一読の価値あり。
日の出・日の入りの計算―天体の出没時刻の求め方 | |
 | 長沢 工 地人書館 1999-12 売り上げランキング : 291817 Amazonで詳しく見る by G-Tools |
このコードで計算できるのは、日本近傍の日の出・日の入りに限定される。他の地域へ拡張するためには、時刻系の違い、位置座標の違い等を考慮した修正が必要になる。太陽が一日中出なかったり、沈まなかったりする地域でもこのままでは（収束解が得られず）使えない。
from math import * class Sun(object): def __init__(self, year, month, day, longitude, latitude): self.year = year self.month = month self.day = day self.longitude = longitude self.latitude = latitude self.hour = 12.0 self.dtr = pi/180.0 def julius_year(self): # obtain eqs.(3.13)(3.14)(3.15) # --- year = self.year month = self.month day = self.day hour = self.hour if month == 1 or month == 2: year -= 1 month = 12 + month k0 = 365*(year-2000) + 30*month + day -33.875\ + int(3*(month+1)/5) + int((year-2000)/4) dt = 64.0 + (year - 1999) k1 = k0 + hour/24.0 + dt/86400.0 timej = k1/365.25 return timej def distance(self): # obtain r in table 3.11 # --- dtr = self.dtr t = self.julius_year() q = (0.007256 - 0.0000002*t)*sin(dtr*(267.54 + 359.991*t)) q += 0.000091*sin(dtr*(265.1+ 719.98*t)) q += 0.000030*sin(dtr*(90.0)) q += 0.000013*sin(dtr*(27.8+ 4452.67*t)) q += 0.000007*sin(dtr*(254.0+ 450.4*t)) q += 0.000007*sin(dtr*(156.0+ 329.6*t)) dis = 10**q return dis def lambda_s(self): # obtain lambda_s in table 3.11 # --- dtr = self.dtr t = self.julius_year() val = 280.4603 + 360.00769*t val += (1.9146 - 0.00005*t)*sin(dtr*(357.538+359.991*t)) val += 0.0200*sin(dtr*(355.05+719.981*t)) val += 0.0048*sin(dtr*(234.95+ 19.341*t)) val += 0.0020*sin(dtr*(247.1+ 329.64*t)) val += 0.0018*sin(dtr*(297.8+ 4452.67*t)) val += 0.0018*sin(dtr*(251.3+ 0.20*t)) val += 0.0015*sin(dtr*(343.2+ 450.37*t)) val += 0.0013*sin(dtr*(81.4+ 225.18*t)) val += 0.0008*sin(dtr*(132.5+ 659.29*t)) val += 0.0007*sin(dtr*(153.3+ 90.38*t)) val += 0.0007*sin(dtr*(206.8+ 30.35*t)) val += 0.0006*sin(dtr*(29.8+ 337.18*t)) val += 0.0005*sin(dtr*(207.4+ 1.50*t)) val += 0.0005*sin(dtr*(291.2+ 22.81*t)) val += 0.0004*sin(dtr*(234.9+ 315.56*t)) val += 0.0004*sin(dtr*(157.3+ 299.30*t)) val += 0.0004*sin(dtr*(21.1+ 720.02*t)) val += 0.0003*sin(dtr*(352.5+ 1079.97*t)) val += 0.0003*sin(dtr*(329.7+ 44.43*t)) while val > 360.0 or val < 0.0: if val < 0.0: val += 360.0 if val > 360.0: val -= 360.0 return val def epsilon(self): # obtain eq.(3.17) # --- t = self.julius_year() eps = 23.439291 -0.000130042*t return eps def alpha(self): # obtain eq.(3.16) # --- dtr = self.dtr eps = self.epsilon() lam_s = self.lambda_s() tana = tan(dtr*lam_s)*cos(dtr*eps) val = atan(tana)/dtr if lam_s < 180.0: if val < 0.0: val += 180.0 else: if val > 0.0: val += 180.0 else: val += 360.0 return val def delta(self): # obtain eq.(3.16) # --- dtr = self.dtr eps = self.epsilon() lam_s = self.lambda_s() sind = sin(dtr*lam_s)*sin(dtr*eps) val = asin(sind)/dtr return val def theta(self): # obtain eq.(3.18) # --- ph = self.hour/24.0 elong = self.longitude t = self.julius_year() a1 = 325.4606 a2 = 360.007700536 a3 = 0.00000003879 val = a1 + a2*t + a3*t**2 + 360.0*ph + elong while val<0.0 or val>360.0: if val > 360.0: val -= 360 if val < 0.0: val += 360 return val def kodo(self): # obtain eq.(1.19) or p.37 # --- rl = self.distance() sv = 0.266994/rl rv = 0.585556 pv = 0.0024428/rl val = -sv - rv + pv return val def jikaku(self): # obtain eq.(2.17) # --- dtr = self.dtr kv = dtr*self.kodo() dl = dtr*self.delta() lant = dtr*self.latitude c = (sin(kv)-sin(dl)*sin(lant))/ (cos(dl)*cos(lant)) val = acos(c)/dtr return val def jikaku2(self): # obtain eq.(2.18) # --- return self.theta() - self.alpha() def getRiseTime(self): self.hour = 12.0 for i in xrange(10): tk = self.jikaku() if tk > 0.0: tk = -tk dd = (tk - self.jikaku2())/360.0 dd = dd - int(dd) d = self.hour/24.0 + dd if d < 0.0: d = d + 1.0 self.hour = d * 24.0 if abs(dd) < 0.00005: break dval = d*24.0 ohour = int(dval) ominu = int(round((dval - ohour)*60.0)) return (ohour, ominu) def getSetTime(self): self.hour = 12.0 for i in xrange(10): tk = self.jikaku() if tk < 0.0: tk = -tk dd = (tk - self.jikaku2())/360.0 dd = dd - int(dd) d = self.hour/24.0 + dd if d < 0.0: d = d + 1.0 self.hour = d * 24.0 if abs(dd) < 0.00005: break dval = d*24.0 ohour = int(dval) ominu = int(round((dval - ohour)*60.0)) return (ohour, ominu) if __name__ == '__main__': year = 2016 day = 3 longitude = 139.7447 # toukei latitude = 35.6544 # hokui for month in range(1,13): sun = Sun(year, month, day, longitude, latitude) trise= sun.getRiseTime() tset = sun.getSetTime() hr1 = str(trise[0]) mn1 = str(trise[1]) hr2 = str(tset[0]) mn2 = str(tset[1]) if len(hr1) == 1: hr1 = "0" + hr1 if len(hr2) == 1: hr2 = "0" + hr2 if len(mn1) == 1: mn1 = "0" + mn1 if len(mn2) == 1: mn2 = "0" + mn2 print '%s/%s/%s :[rise time]=%s:%s [set time]=%s:%s' % ( year,month,day,hr1,mn1,hr2,mn2)
実行することにより、以下の結果が得られる。
[work]$ python calculateSun.py 2016/1/3 :[rise time]=06:51 [set time]=16:40 2016/2/3 :[rise time]=06:40 [set time]=17:10 2016/3/3 :[rise time]=06:08 [set time]=17:38 2016/4/3 :[rise time]=05:25 [set time]=18:04 2016/5/3 :[rise time]=04:47 [set time]=18:29 2016/6/3 :[rise time]=04:26 [set time]=18:53 2016/7/3 :[rise time]=04:30 [set time]=19:01 2016/8/3 :[rise time]=04:50 [set time]=18:43 2016/9/3 :[rise time]=05:14 [set time]=18:06 2016/10/3 :[rise time]=05:37 [set time]=17:22 2016/11/3 :[rise time]=06:05 [set time]=16:44 2016/12/3 :[rise time]=06:34 [set time]=16:28
グラフにすると下の通り。
横軸は月、縦軸は時刻を表す。上図は日の出の時刻、下図は日の入りの時刻。
