Hatena::ブログ(Diary)

理想のユーザ・インターフェイスを求めて このページをアンテナに追加 RSSフィード

2016-02-03

日の出・日の入りの時刻を求める計算

日の出と日の入りの時刻を求めるコードを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

グラフにすると下の通り。

f:id:Megumi221:20160203214955p:image

f:id:Megumi221:20160203214950p:image

横軸は月、縦軸は時刻を表す。上図は日の出の時刻、下図は日の入りの時刻。