geo/sun.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
"""
geo.sun (Sun)
=============

**Author:**

* Dirk Alders <sudo-dirk@mount-mockery.de>
* Formula from Dr. Roland Brodbeck, Calsky (http://lexikon.astronomie.info/zeitgleichung/neu.html)
* First implementation by Alexander Klupp 2014-01-14

**Description:**

    This module is a submodule of :mod:`geo` and supports functions and classes for sun position.

**Contentlist:**

* :func:`geo.sun.sunset`
* :func:`geo.sun.sunrise`

**Unittest:**

    See also the :download:`unittest <../../geo/_testresults_/unittest.pdf>` documentation.
"""


import calendar
import math
import time

europe_berlin = {'lat': 52.520008, 'lon': 13.404954}


def JulianischesDatum(Jahr, Monat, Tag, Stunde, Minuten, Sekunden):
    if (Monat <= 2):
        Monat = Monat + 12
        Jahr = Jahr - 1
    Gregor = (Jahr / 400) - (Jahr / 100) + (Jahr / 4)   # Gregorianischer Kalender
    return 2400000.5 + 365 * Jahr - 679004 + Gregor + math.floor(30.6001 * (Monat + 1)) + Tag + Stunde / 24 + Minuten / 1440 + Sekunden / 86400


def InPi(x):
    n = int(x / (2 * math.pi))
    x = x - n * 2 * math.pi
    if (x < 0):
        x += 2 * math.pi
    return x


def eps(T):
    # Neigung der Erdachse
    return math.pi / 180 * (23.43929111 + (-46.8150 * T - 0.00059 * T ** 2 + 0.001813 * T ** 3) / 3600)


def BerechneZeitgleichung(T):
    RA_Mittel = 18.71506921 + 2400.0513369 * T + (2.5862e-5 - 1.72e-9 * T) * T ** 2
    M = InPi(2 * math.pi * (0.993133 + 99.997361 * T))
    L = InPi(2 * math.pi * (0.7859453 + M / (2 * math.pi) + (6893 * math.sin(M) + 72 * math.sin(2 * M) + 6191.2 * T) / 1296e3))
    e = eps(T)
    RA = math.atan(math.tan(L) * math.cos(e))
    if (RA < 0):
        RA += math.pi
    if (L > math.pi):
        RA += math.pi
    RA = 24 * RA / (2 * math.pi)
    DK = math.asin(math.sin(e) * math.sin(L))
    # Damit 0 <= RA_Mittel < 24
    RA_Mittel = 24.0 * InPi(2 * math.pi * RA_Mittel / 24.0) / (2 * math.pi)
    dRA = RA_Mittel - RA
    if (dRA < -12.0):
        dRA += 24.0
    if (dRA > 12.0):
        dRA -= 24.0
    dRA = dRA * 1.0027379
    return dRA, DK


def Sonnenauf_untergang(JD, Zeitzone, coord):
    # Zeitzone = 0 #Weltzeit
    # Zeitzone = 1 #Winterzeit
    # Zeitzone = 2 #Sommerzeit
    # JD = JulianischesDatum

    B = math.radians(coord.get('lat'))          # geographische Breite Erkelenz
    GeographischeLaenge = coord.get('lon')      # geographische Laenge

    JD2000 = 2451545
    h = -50.0 / 60.0 * math.pi / 180
    T = (JD - JD2000) / 36525

    Zeitgleichung, DK = BerechneZeitgleichung(T)

    Zeitdifferenz = 12 * math.acos((math.sin(h) - math.sin(B) * math.sin(DK)) / (math.cos(B) * math.cos(DK))) / math.pi

    AufgangOrtszeit = 12 - Zeitdifferenz - Zeitgleichung
    UntergangOrtszeit = 12 + Zeitdifferenz - Zeitgleichung
    AufgangWeltzeit = AufgangOrtszeit - GeographischeLaenge / 15
    UntergangWeltzeit = UntergangOrtszeit - GeographischeLaenge / 15

    Aufgang = AufgangWeltzeit + Zeitzone
    if (Aufgang < 0):
        Aufgang += 24
    elif (Aufgang >= 24):
        Aufgang -= 24

    AM = round(Aufgang * 60) / 60       # minutengenau runden

    Untergang = UntergangWeltzeit + Zeitzone
    if (Untergang < 0):
        Untergang += 24
    elif (Untergang >= 24):
        Untergang -= 24

    UM = round(Untergang * 60) / 60     # minutengenau runden

    return AM, UM


def sunrise(coord=europe_berlin, date=None):
    """
    :param coord: Target coordinate or None (default is central europe).
    :type coord: geo.gps.coordinate
    :param date: The day to calculate with or None (only year, month and day are relevant; default ist today)
    :type date: time.struct_time
    :return: The date and time information for the sunrise
    :rtype: time.struct_time

    This Method calculates the time for sunrise for a given date and coordinate.

    .. code-block:: python

        >>> import geo

        >>> ...
    """
    date = date or time.localtime()

    year, month, day = date[0:3]
    dst = date[8]                             # Sommerzeit

    AM = Sonnenauf_untergang(JulianischesDatum(year, month, day, 12, 0, 0), dst + 1, coord)[0]

    tup_list = list(date)
    tup_list[3] = 0
    tup_list[4] = 0
    tup_list[5] = 0
    tm = time.mktime(time.struct_time(tup_list))
    return time.localtime(tm + AM * 3600)


def sunset(coord=europe_berlin, date=None):
    """
    :param coord: Target coordinate or None (default is central europe).
    :type coord: geo.gps.coordinate
    :param date: The day to calculate with or None (only year, month and day are relevant; default ist today)
    :type date: time.struct_time
    :return: The date and time information for the sunset
    :rtype: time.struct_time

    This Method calculates the time for sunrise for a given date and coordinate.

    .. code-block:: python

        >>> import geo

        >>> ...
    """
    date = date or time.localtime()

    year, month, day = date[0:3]
    dst = date[8]                             # Sommerzeit

    UM = Sonnenauf_untergang(JulianischesDatum(year, month, day, 12, 0, 0), dst + 1, coord)[1]

    tup_list = list(date)
    tup_list[3] = 0
    tup_list[4] = 0
    tup_list[5] = 0
    tm = time.mktime(time.struct_time(tup_list))
    return time.localtime(tm + UM * 3600)