Demo entry 1463097

sunrize_sunset

   

Submitted by anonymous on Mar 04, 2015 at 20:59
Language: Python. Code size: 7.2 kB.

#
# http://hhsprings.pinoko.jp/site-hhs/2015/02/%E6%97%A5%E5%87%BA%E6%97%A5%E6%B2%A1%E8%A8%88%E7%AE%97%E3%80%81%E3%82%84%E3%81%A3%E3%81%A6%E3%81%BF%E3%82%88%E3%81%86/
#
import numpy as np

# sin function using degree
def sind(d):
    return np.sin(np.radians(d))

# cos function using degree
def cosd(d):
    return np.cos(np.radians(d))

# tan function using degree
def tand(d):
    return np.tan(np.radians(d))

# calculate Julius year (year from 2000/1/1, for variable "t")
def jy(yy, mm, dd, h, m, s, i): # yy/mm/dd h:m:s, i: time difference
    yy -= 2000
    if (mm <= 2):
        mm += 12
        yy -= 1

    k = 365 * yy + 30 * mm + dd - 33.5 - i / 24.0 + np.floor(3 * (mm + 1) / 5.0) \
        + np.floor(yy / 4.0) - np.floor(yy / 100.0) + np.floor(yy / 400.0)
    k += ((s / 60.0 + m) / 60 + h) / 24.0 # plus time
    k += (65 + yy) / 86400.0 # plus delta T
    return k / 365.25

# solar position1 (celestial longitude, degree)
def spls(t): # t: Julius year
    L = 280.4603 + 360.00769 * t \
        + (1.9146 - 0.00005 * t) * sind(357.538 + 359.991 * t) \
        + 0.0200 * sind(355.05 +  719.981 * t) \
        + 0.0048 * sind(234.95 +   19.341 * t) \
        + 0.0020 * sind(247.1  +  329.640 * t) \
        + 0.0018 * sind(297.8  + 4452.67  * t) \
        + 0.0018 * sind(251.3  +    0.20  * t) \
        + 0.0015 * sind(343.2  +  450.37  * t) \
        + 0.0013 * sind( 81.4  +  225.18  * t) \
        + 0.0008 * sind(132.5  +  659.29  * t) \
        + 0.0007 * sind(153.3  +   90.38  * t) \
        + 0.0007 * sind(206.8  +   30.35  * t) \
        + 0.0006 * sind( 29.8  +  337.18  * t) \
        + 0.0005 * sind(207.4  +    1.50  * t) \
        + 0.0005 * sind(291.2  +   22.81  * t) \
        + 0.0004 * sind(234.9  +  315.56  * t) \
        + 0.0004 * sind(157.3  +  299.30  * t) \
        + 0.0004 * sind( 21.1  +  720.02  * t) \
        + 0.0003 * sind(352.5  + 1079.97  * t) \
        + 0.0003 * sind(329.7  +   44.43  * t)
    while (L >= 360):
        L -= 360
    while (L < 0):
        L += 360
    return L

# solar position2 (distance, AU)
def spds(t): # t: Julius year
    r = (0.007256 - 0.0000002 * t) * sind(267.54 + 359.991 * t) \
        + 0.000091 * sind(265.1 +  719.98 * t) \
        + 0.000030 * sind( 90.0) \
        + 0.000013 * sind( 27.8 + 4452.67 * t) \
        + 0.000007 * sind(254   +  450.4  * t) \
        + 0.000007 * sind(156   +  329.6  * t)
    r = np.power(10, r)
    return r

# solar position3 (declination, degree)
def spal(t): # t: Julius year
    ls = spls(t) ;
    ep = 23.439291 - 0.000130042 * t
    al = np.degrees(np.arctan(tand(ls) * cosd(ep)))
    if ((ls >= 0) and (ls < 180)):
        while (al < 0):
            al += 180
        while (al >= 180):
            al -= 180
    else:
        while (al < 180):
            al += 180
        while (al >= 360):
            al -= 180
    return al ;

# solar position4 (the right ascension, degree)
def spdl(t): # t: Julius year
    ls = spls(t)
    ep = 23.439291 - 0.000130042 * t
    dl = np.degrees(np.arcsin(sind(ls) * sind(ep)))
    return dl

# Calculate sidereal hour (degree)
def sh(t, h, m, s, l, i):
    # t: julius year, h: hour, m: minute, s: second,
    # l: longitude, i: time difference
    d = ((s / 60.0 + m) / 60.0 + h) / 24.0 # elapsed hour (from 0:00 a.m.)
    th = 100.4606 + 360.007700536 * t + 0.00000003879 * t * t - 15 * i
    th += l + 360 * d
    while (th >= 360):
        th -= 360
    while (th < 0):
        th += 360
    return th

# Calculating the seeming horizon altitude "sa"(degree)
def eandp(alt, ds): # subfunction for altitude and parallax
    e = 0.035333333 * np.sqrt(alt)
    p = 0.002442818 / ds
    return p - e

def sa(alt, ds): # alt: altitude (m), ds: solar distance (AU)
    s = 0.266994444 / ds
    r = 0.585555555
    k = eandp(alt,ds) - s - r
    return k

# Calculating solar alititude (degree) {
def soal(la, th, al, dl):
    # la: latitude, th: sidereal hour,
    # al: solar declination, dl: right ascension
    h = sind(dl) * sind(la) + cosd(dl) * cosd(la) * cosd(th - al)
    h = np.degrees(np.arcsin(h))
    return h

# Calculating solar direction (degree) {
def sodr(la, th, al, dl):
    # la: latitude, th: sidereal hour,
    # al: solar declination, dl: right ascension
    t = th - al
    dc = -cosd(dl) * sind(t)
    dm = sind(dl) * sind(la) - cosd(dl) * cosd(la) * cosd(t)
    if (dm == 0):
        st = sind(t)
        if (st > 0):
            dr = -90
        if (st == 0):
            dr = 9999
        if (st < 0):
            dr = 90
    else:
        dr = np.degrees(np.arctan(dc / dm))
        if (dm < 0):
            dr += 180

    if (dr < 0):
        dr += 360
    return dr

def calc(yy, mm, dd, i, lon, lat, alt=0, hh=-1, m=-1):
    """
    return list:
        item 0: astronomical twilight start
        item 1: voyage twilight start
        item 2: citizen twilight start
        item 3: sun rise
        item 4: meridian
        item 5: sun set
        item 6: citizen twilight end
        item 7: voyage twilight end
        item 8: astronomical twilight end

    each item has time, and solar direction (in degrees) if item is sun rise or sun set,
    solar altitude (in degrees) if item is meridian.
    """
    result = []

    t = jy(yy, mm, dd - 1, 23, 59, 0, i)
    th = sh(t, 23, 59, 0, lon, i)
    ds = spds(t)
    ls = spls(t)
    alp = spal(t)
    dlt = spdl(t)
    pht = soal(lat, th, alp, dlt)
    pdr = sodr(lat, th, alp, dlt)

    for hh in range(0, 24):
        for m in range(0, 60):
            t = jy(yy, mm, dd, hh, m, 0, i)
            th = sh(t, hh, m, 0, lon, i)
            ds = spds(t)
            ls = spls(t)
            alp = spal(t)
            dlt = spdl(t)
            ht = soal(lat, th, alp, dlt)
            dr = sodr(lat, th, alp, dlt)
            tt = eandp(alt, ds)
            t1 = tt - 18
            t2 = tt - 12
            t3 = tt - 6
            t4 = sa(alt, ds)

            key = "%d-%02d-%02d %02d:%02d:00+%02d" % (yy, mm, dd, hh, m, i)

            if ((pht < t1) and (ht > t1)):
                result.append((key, None))
            if ((pht < t2) and (ht > t2)):
                result.append((key, None))
            if ((pht < t3) and (ht > t3)):
                result.append((key, None))
            if ((pht < t4) and (ht > t4)):
                result.append((key, np.floor(dr)))
            if ((pdr < 180) and (dr > 180)):
                result.append((key, np.floor(ht)))
            if ((pht > t4) and (ht < t4)):
                result.append((key, np.floor(dr)))
            if ((pht > t3) and (ht < t3)):
                result.append((key, None))
            if ((pht > t2) and (ht < t2)):
                result.append((key, None))
            if ((pht > t1) and (ht < t1)):
                result.append((key, None))
            pht = ht
            pdr = dr

    #f.result.value = ans
    return result

if __name__ == '__main__':
    r = calc(2013, 4, 2, 9, 139.7909446414986, 35.54131474808071, 0.)
    for k, v in r:
        print("%s: %s" % (k, v))

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).