#! /usr/bin/python3
import math
import svgwrite

def generate_sundial_with_analemma(latitude_deg, nodus_height_mm=100, filename="sonnenuhr_analemma.svg"):
    phi = math.radians(latitude_deg)
    radius = 200  # Radius des Zifferblatts in mm
    center_x, center_y = 250, 150  # Fußpunkt B (Achsursprung)
    
    dwg = svgwrite.Drawing(filename, size=('500mm', '600mm'), viewBox='0 0 500 600')
    
    # 1. Zifferblatt-Rahmen & Fußpunkt (B)
    dwg.add(dwg.circle(center=(center_x, center_y), r=3, fill='red'))
    
    # 2. Stundenlinien (6 bis 18 Uhr)
    for hour in range(6, 19):
        t = hour - 12
        hour_angle = math.radians(t * 15)
        
        # Vertikal-Südsonnenuhr: tan(gamma) = cos(phi) * tan(h)
        gamma = math.atan(math.cos(phi) * math.tan(hour_angle))
        
        x = center_x + radius * math.sin(gamma)
        y = center_y + radius * math.cos(gamma)
        
        dwg.add(dwg.line(start=(center_x, center_y), end=(x, y), stroke='#888888', stroke_width=1, stroke_dasharray="4,4" if hour != 12 else "none"))
        
        # Beschriftung
        text_x = center_x + (radius + 15) * math.sin(gamma)
        text_y = center_y + (radius + 15) * math.cos(gamma)
        dwg.add(dwg.text(f"{hour}", insert=(text_x - 5, text_y + 5), font_size="14px", font_family="Arial", font_weight="bold"))

    # 3. Analemma-Kurve für jeden Tag des Jahres berechnen
    analemma_points = []
    
    for day in range(1, 366):
        # Näherungsformeln für Deklination (delta) und Zeitgleichung (ZGL in Minuten)
        # B = Winkelposition der Erde auf ihrer Bahn
        B = math.radians((360 / 365) * (day - 81))
        
        # Deklination delta in Radian
        delta = math.radians(23.44 * math.sin(math.radians((360 / 365) * (day - 80))))
        
        # Zeitgleichung ZGL in Minuten
        eot_minutes = 9.87 * math.sin(2 * B) - 7.53 * math.cos(B) - 1.5 * math.sin(B)
        
        # Umrechnung ZGL in Stundenwinkel-Abweichung (1 Min = 0.25 Grad = 0.00436 Rad)
        h_eot = math.radians(eot_minutes * 0.25)
        
        # Projektion des Nodus-Schattens auf die vertikale Südwand:
        # y_shadow: Vertikaler Abstand vom Noduspunkt (Sonne höher -> Schatten weiter unten)
        # x_shadow: Horizontaler Versatz durch ZGL
        y_shadow = center_y + nodus_height_mm * math.tan(phi - delta)
        x_shadow = center_x + nodus_height_mm * (math.sin(h_eot) / math.cos(phi - delta))
        
        analemma_points.append((x_shadow, y_shadow))

    # 4. Analemma als geschlossener Pfad (Polyline) einzeichnen
    dwg.add(dwg.polyline(points=analemma_points, stroke='red', stroke_width=1.5, fill='none'))
    
    # Optional: Monats-Markierungen hervorheben (z. B. jeweils am 1. des Monats)
    month_days = [1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335]
    month_names = ["Jan", "Feb", "Mär", "Apr", "Mai", "Jun", "Jul", "Aug", "Sep", "Okt", "Nov", "Dez"]
    
    for idx, d in enumerate(month_days):
        pt = analemma_points[d - 1]
        dwg.add(dwg.circle(center=pt, r=2, fill='blue'))
        dwg.add(dwg.text(month_names[idx], insert=(pt[0] + 6, pt[1] + 3), font_size="10px", font_family="Arial", fill='blue'))

    dwg.save()
    print(f"SVG mit Analemma erfolgreich gespeichert unter: {filename}")

# Aufruf für z. B. München (Breitengrad 48.1° N), Nodus-Höhe 100mm
generate_sundial_with_analemma(latitude_deg=48.1, nodus_height_mm=100)