Sonnenzenitwinkel berechnen

[gkuhl]

Hallo allerseits,

ich möchte hier mal ein kleines Skript zeigen, dass den Sonnenzenitwinkel berechnet. Ich würde mich über eure Anregungen und Kritiken freuen.

Code: Alles auswählen

from math import acos, cos, degrees, radians, sin
from datetime import datetime

# define triogonometry with degrees
cos2 = lambda x: cos(radians(x))
sin2 = lambda x: sin(radians(x))
acos2 = lambda x: degrees(acos(x))

def sza(UTC = datetime.utcnow(), latitude = 52.37, longitude = 9.72):
    '''Returns the solar zenith angle (in degree)
    
    UTC         (as datetime.datetime Object)
    longitude   (in degree)
    latitude    (in degree)
    
    Default values: Hannover, Germany
    '''
      
    # parameter 
    day_of_year = UTC.timetuple().tm_yday
    leap_year_factor = (-0.375, 0.375, -0.125, 0.125)[UTC.year % 4]   
    UTC_min = UTC.hour * 60.0 + UTC.minute + UTC.second / 60.0
    J = 360.0 / 365.0 * (day_of_year + leap_year_factor + UTC_min / 1440.0)

    # hour angle (using the time equation)
    average_localtime = UTC_min + 4 * longitude
    true_solar_time = average_localtime + 0.0066 + 7.3525 * cos2(  J +  85.9) \
                                                 + 9.9359 * cos2(2*J + 108.9) \
                                                 + 0.3387 * cos2(3*J + 105.2)          
    hour_angle = 15.0 * (720 - true_solar_time) / 60.0
    
    # sun declination (using DIN 5034-2)
    declination = 0.3948 - 23.2559 * cos2(  J + 9.1 ) \
                         -  0.3915 * cos2(2*J + 5.4 ) \
                         -  0.1764 * cos2(3*J + 26.0)
    
    # solar zenith angle
    return acos2( sin2(latitude) * sin2(declination) + 
                  cos2(hour_angle) * cos2(latitude) * cos2(declination) 
                )   
    
if __name__ == '__main__':
    print sza()

Grüße
Gerrit

[mkesper]

Woher kommen die Konstanten bei true_solar_time und declination? Ist so pure Magie.
Die Backslashes kannst du dir sparen, da die Ausdrücke in Klammern stehen.

Code: Alles auswählen

from __future__ import division

würde ich noch empfehlen.

[wuf]

Hallo mkesper

Die Backslashes kannst du dir sparen, da die Ausdrücke in Klammern stehen

.

Hast du das Skript selber einmal ohne die Backslashes ausprobiert?

Gruss wuf

[gkuhl]

Woher kommen die Konstanten bei true_solar_time und declination?

Beide Gleichungen sind empirisch und in der DIN 5034-2 zu finden. Ohne die Backslash wäre das Programm nicht lauffähig, da es keine Klammern gibt. Warum ich das neue Verhalten von Division importieren soll, verstehe ich nicht. Wieso glaubst du, dass das sinnvoll wäre?

Grüße
Gerrit

[ms4py]

Warum ich das neue Verhalten von Division importieren soll, verstehe ich nicht. Wieso glaubst du, dass das sinnvoll wäre?

Dann kannst du auf die Float-Konvertierung àla 60.0 verzichten.

[gkuhl]

Dann kannst du auf die Float-Konvertierung àla 60.0 verzichten.

Das ist mir bekannt, allerdings hätte ich die .0 vermutlich auch so verwendet. Mir ist es lieber, wenn jeder sofort sieht, dass hier mit Floats gerechnet wird. Das neue Verhalten von ``/`` importiere ich meist nur, wenn ich Integer dividieren muss. Zum Beispiel beim Normieren von Verteilungen, wo der zusätzliche Aufruf von ``float()`` nervt:

Code: Alles auswählen

In [15]: y,x = histogram([1,2,3,1,2])

In [16]: y
Out[16]: array([2, 0, 0, 0, 0, 2, 0, 0, 0, 1])

In [17]: x
Out[17]: array([ 1. ,  1.2,  1.4,  1.6,  1.8,  2. ,  2.2,  2.4,  2.6,  2.8,  3. ])

In [18]: y / sum(y)
Out[18]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

In [19]: y / float(sum(y))
Out[19]: array([ 0.4,  0. ,  0. ,  0. ,  0. ,  0.4,  0. ,  0. ,  0. ,  0.2])

Ein anderer Grund ist, dass ich das Skript möglichst einfach halten wollte um Neulinge nicht unnötig zu verwirren. Ich hoffe noch darauf, dass in meinem Institut IDL langsam von Python verdrängt wird.

Grüße
Gerrit

[jens]

Man könnte das hier:

Code: Alles auswählen

    true_solar_time = average_localtime + 0.0066 + 7.3525 * cos2(  J +  85.9) \
                                                 + 9.9359 * cos2(2*J + 108.9) \
                                                 + 0.3387 * cos2(3*J + 105.2) 

auch so schreiben:

Code: Alles auswählen

    true_solar_time = average_localtime + 0.0066
    true_solar_time += 7.3525 * cos2(  J +  85.9)
    true_solar_time += 9.9359 * cos2(2*J + 108.9)
    true_solar_time += 0.3387 * cos2(3*J + 105.2)    

[mkesper]

Okay, sorry, hatte ich wohl Tomaten auf den Augen.
Die neue Division finde ich einfach logischer, weil sie der mathematischen entspricht.
Leider gibt es keine Onlinequellen der zitierten DIN ("Tageslicht in Innenräumen"), so dass man das einfach glauben muss.

[wuf]

...... Hier noch eine Variante um Backslashes zu umgehen:

Code: Alles auswählen

true_solar_time = (average_localtime + 0.0066 + 7.3525 * cos2(J + 85.9)
    + 9.9359 * cos2(2*J + 108.9) + 0.3387 * cos2(3*J + 105.2))

Gruss wuf