Jak výpočítat východ a západ slunce v pythonu

Ahoj,

tak jsem dostal za úkol vypočítat, kdy vychází slunce a kdy zapadá aby mohl ovládat relé které vypíná spotřebiče.

zde vidíte jak to vypadá. Program jsem našel na netu, ale byl napsaný v php, tak jsem to předělával do pythonu. Možná tam mám chybu, tak kdyby jste na něco přišli dejte vědet.

http://zhola.com/vychodySlunce/formular.php

Program:

#!/usr/bin/env python

# RASPBERRY PI VERSION

import time
import math

sirka = 50
minut1 = 05
vterin1 = 27
sirka2  = "severni"

delka  = 14 
minut2 = 24
vterin2 = 02
delka2 = "vychodni"

rok = 2016
mesic = 8
den = 26
zona = 1

if sirka2=="jizni": sirka = -sirka 
if delka2=="zapadni": delka = -delka

def SGN(VALUE):
  if (VALUE == 0):
   return 0 
  elif (VALUE > 0):
    return 1 
  else : return -1    

1
P1 = 3.14159265 
P2 = 2 * P1 
DR = P1 / 180 
K1 = 15 * DR * 1.0027379 

S_MESSAGE = "Zapad Slunce:" 
R_MESSAGE = "Vychod Slunce:" 
M1_MESSAGE = "Nedochazi k vychodu Slunce." 
M2_MESSAGE = "Nedochazi k zapadu Slunce." 
M3_MESSAGE = "Slunce nevychazi (polarni noc)" 
M4_MESSAGE = "Slunce sviti cely den (polarni den)." 
minut1 = minut1 + ((1/60)*vterin1) 
minut2 = minut2 + ((1/60)*vterin2) 
sirka = sirka + ((1/60)*minut1) 
delka = delka + ((1/60)*minut2) 
B5 = sirka 
L5 = delka 
H = float(-zona) 
L5 = L5 / 360 
Z0 = H / 24 
Y = rok 
M = mesic 
D = den 
G = 1 
if (Y < 1583) : G = 0 
D1 = D 
F = D - D1 - 0.5 
J = -int(7 * (int((M + 9) / 12) + Y) / 4) 
if G <> 0:
 S = SGN(M-9) 
 A = abs(M-9) 
 J3 = int(Y + S * int(A / 7)) 
 J3 = -int((int(J3 / 100) + 1) * 3 / 4) 

J = J + int(275 * M / 9) + D1 + G * J3 
J = J + 1721027 + 2 * G + 367 * Y 
if F < 0: 
 F = F + 1 
 J = J - 1 

T = (J - 2451545) + F 
TT = T / 36525 + 1 
T0 = T / 36525 
S = 24110.5 + 8640184.812999999 * T0 
S = S + 86636.6 * Z0 + 86400 * L5 
S = S / 86400; 
S = S - int(S) 
T0 = S * 360 * DR 
T = T + Z0 
L = 0.779072 + 0.00273790931 * T 
G = 0.993126 + 0.0027377785 * T 
L = L - int(L) 
G = G - int(G) 
L = L * P2 
G = G * P2 
V = 0.39785 *  math.sin(L) 
V = V - 0.01 *  math.sin(L - G) 
V = V + 0.00333 *  math.sin(L + G) 
V = V - 0.00021 * TT *  math.sin(L) 
U = 1 - 0.03349 *  math.cos(G) 
U = U - 0.00014 *  math.cos(2 * L) 
U = U + 0.00008 *  math.cos(L)
W = -0.0001 - 0.04129 *  math.sin(2 * L) 
W = W + 0.03211 *  math.sin(G) 
W = W + 0.00104 *  math.sin(2 * L - G) 
W = W - 0.00035 * math. sin(2 * L + G) 
W = W - 0.00008 * TT *  math.sin(G) 
S = W /  math.sqrt(U - V * V) 
A5 = L +  math.atan(S /  math.sqrt(1 - S * S)) 
S = V / math.sqrt(U) 
D5 =  math.atan(S /  math.sqrt(1 - S * S)) 
R5 = 1.00021 *  math.sqrt(U) 
A_1 = A5 
D_1 = D5 
T = T + 1 
L = 0.779072 + 0.00273790931 * T 
G = 0.993126 + 0.0027377785 * T 
L = L - int(L) 
G = G - int(G) 
L = L * P2 
G = G * P2 
V = 0.39785 *  math.sin(L) 
V = V - 0.01 *  math.sin(L - G) 
V = V + 0.00333 *  math.sin(L + G) 
V = V - 0.00021 * TT *  math.sin(L) 
U = 1 - 0.03349 *  math.cos(G) 
U = U - 0.00014 *  math.cos(2 * L) 
U = U + 0.00008 *  math.cos(L) 
W = -0.0001 - 0.04129 * math. sin(2 * L) 
W = W + 0.03211 *  math.sin(G) 
W = W + 0.00104 *  math.sin(2 * L - G) 
W = W - 0.00035 *  math.sin(2 * L + G) 
W = W - 0.00008 * TT * math.sin(G) 
S = W /  math.sqrt(U - V * V) 
A5 = L +  math.atan(S /  math.sqrt(1 - S * S)) 
S = V /  math.sqrt(U) 
D5 =  math.atan(S /  math.sqrt(1 - S * S)) 
R5 = 1.00021 *  math.sqrt(U) 
A_2 = A5 
D_2 = D5 
if A_2 < A_1: A_2 = A_2 + P2 
 
Z1 = DR * 90.833 
S =  math.sin(B5 * DR) 
C =  math.cos(B5 * DR) 
Z =  math.cos(Z1) 
M8 = 0 
W8 = 0 
A0 = A_1 
D0 = D_1 
DA = A_2 - A_1 
DD = D_2 - D_1 
for C0 in range (0 , 23):
 P = (C0 + 1) / 24
 A2 = A_1 + P * DA 
 D2 = D_1 + P * DD 
 L0 = T0 + C0 * K1 
 L2 = L0 + K1 
 H0 = L0 - A0 
 H2 = L2 - A2 
 H1 = (H2 + H0) / 2 
 D1 = (D2 + D0) / 2 
 if C0 <= 0: V0 = S *  math.sin(D0) + C *  math.cos(D0) *  math.cos(H0) - Z 
 
 V2 = S *  math.sin(D2) + C *  math.cos(D2) *  math.cos(H2) - Z 
 if (SGN(V0) <> SGN(V2)):
  V1 = S *  math.sin(D1) + C *  math.cos(D1) *  math.cos(H1) - Z 
  A = 2 * V2 - 4 * V1 + 2 * V0 
  B = 4 * V1 - 3 * V0 - V2 
  D = B * B - 4 * A * V0 
  if (D >= 0): 
   D =  math.sqrt(D) 
   if V0 < 0 and V2 > 0 : print R_MESSAGE 
   if V0 < 0 and V2 > 0 : M8 = 1 
   if V0 > 0 and V2 < 0 : print S_MESSAGE 
   if V0 > 0 and V2 < 0 : W8 = 1 
   E = (-B + D) / (2 * A) 
   if ((E > 1) or (E < 0)): E = (-B - D) / (2 * A) 
   T3 = C0 + E + 1 / 120 
   H3 = int(T3) 
   M3 = int((T3 - H3) * 60) 
   VYPIS_M3 = M3 
   if (VYPIS_M3 < 10) : VYPIS_M3 = '0'+ str(VYPIS_M3) 
   print " ", H3, ":",VYPIS_M3 
   H7 = H0 + E * (H2 - H0) 
   N7 = - math.cos(D1) *  math.sin(H7) 
   D7 = C *  math.sin(D1) - S *  math.cos(D1) *  math.cos(H7) 
   AZ = math.atan(N7 / D7) / DR 
   if (D7 < 0): AZ = AZ + 180 
   if (AZ < 0): AZ = AZ + 360 
   if (AZ > 360): AZ = AZ - 360 
   print "(azimut = ",int(AZ),"stupnu)" 
 A0 = A2; 
 D0 = D2; 
 V0 = V2; 

if ((M8 == 0) and (W8 == 0)): 
 if (V2 < 0): print M3_MESSAGE 
 if (V2 > 0): print M4_MESSAGE 
else: 
 if (M8 == 0): print M1_MESSAGE 
 if (W8 == 0): print M2_MESSAGE