#! /usr/bin/python3
import math
import matplotlib.pyplot as plt
import numpy as np

def generate_sundial_png_and_table(latitude_deg, nodus_height_mm=100, png_filename="sonnenuhr.png", txt_filename="zgl_tabelle.txt"):
    phi = math.radians(latitude_deg)
    
    # ---------------------------------------------------------
    # 1. ZGL-Tabelle erzeugen und als TXT speichern
    # ---------------------------------------------------------
    days_in_months = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    month_names = ["Januar", "Februar", "März", "April", "Mai", "Juni", 
                   "Juli", "August", "September", "Oktober", "November", "Dezember"]
    
    zgl_day_values = []  # Speichert (Tag_des_Jahres, Datum_Text, ZGL_Minuten)
    
    day_counter = 1
    with open(txt_filename, "w", encoding="utf-8") as f:
        f.write("=========================================================\n")
        f.write(f" MINUTENGENAUE ZEITGLEICHUNG (ZGL) UND TAGESWERTE\n")
        f.write(f" Standort-Breite: {latitude_deg}° N\n")
        f.write("=========================================================\n\n")
        f.write("Datum       | ZGL (Minuten) | Korrektur-Hinweis\n")
        f.write("------------+---------------+-----------------------------------\n")
        
        for m_idx, num_days in enumerate(days_in_months):
            for day_of_month in range(1, num_days + 1):
                # Astronomische Näherung für ZGL
                B = math.radians((360 / 365) * (day_counter - 81))
                eot_minutes = 9.87 * math.sin(2 * B) - 7.53 * math.cos(B) - 1.5 * math.sin(B)
                
                # Runden auf ganze Minuten
                eot_rounded = round(eot_minutes)
                
                date_str = f"{day_of_month:02d}. {month_names[m_idx]:<9}"
                
                if eot_rounded > 0:
                    sign_str = f"+{eot_rounded:2d} min"
                    note = "Sonnenuhr geht VOR  (ZGL abziehen)"
                elif eot_rounded < 0:
                    sign_str = f"{eot_rounded:3d} min"
                    note = "Sonnenuhr geht NACH (ZGL addieren)"
                else:
                    sign_str = "  0 min"
                    note = "Exakt synchron"
                
                f.write(f"{date_str} |    {sign_str}    | {note}\n")
                
                zgl_day_values.append((day_counter, day_of_month, m_idx, eot_minutes))
                day_counter += 1
            f.write("------------+---------------+-----------------------------------\n")
            
    print(f"Tabelle gespeichert unter: {txt_filename}")

    # ---------------------------------------------------------
    # 2. PNG-Grafik mit Matplotlib zeichnen
    # ---------------------------------------------------------
    fig, ax = plt.subplots(figsize=(10, 12), dpi=300)
    ax.set_aspect('equal')
    
    center_x, center_y = 0, 0  # Fußpunkt B
    radius = 200  # Radius für Stundenlinien
    
    # A. Stundenlinien (6 bis 18 Uhr)
    for hour in range(6, 19):
        t = hour - 12
        hour_angle = math.radians(t * 15)
        
        # Formel Vertikal-Südsonnenuhr: tan(gamma) = cos(phi) * tan(h)
        gamma = math.atan(math.cos(phi) * math.tan(hour_angle))
        
        x = radius * math.sin(gamma)
        y = -radius * math.cos(gamma)  # Nach unten abtragen
        
        line_style = '-' if hour in [6, 12, 18] else '--'
        line_color = 'black' if hour == 12 else 'gray'
        
        ax.plot([center_x, x], [center_y, y], color=line_color, linestyle=line_style, linewidth=1)
        
        # Beschriftung
        text_x = (radius + 15) * math.sin(gamma)
        text_y = -(radius + 15) * math.cos(gamma)
        ax.text(text_x, text_y, f"{hour}", fontsize=11, fontweight='bold', ha='center', va='center')

    # B. Analemma-Kurve auf der 12-Uhr-Linie berechnen & zeichnen
    analemma_x = []
    analemma_y = []
    
    for day_counter, day_m, m_idx, eot_min in zgl_day_values:
        # Deklination delta
        delta = math.radians(23.44 * math.sin(math.radians((360 / 365) * (day_counter - 80))))
        h_eot = math.radians(eot_min * 0.25)
        
        # Schattenpunkt des Nodus
        y_shadow = -nodus_height_mm * math.tan(phi - delta)
        x_shadow = nodus_height_mm * (math.sin(h_eot) / math.cos(phi - delta))
        
        analemma_x.append(x_shadow)
        analemma_y.append(y_shadow)

    # Analemma rot einzeichnen
    ax.plot(analemma_x, analemma_y, color='red', linewidth=1.5, label='Analemma (ZGL-Kurve)')
    
    # C. Monatsmarkierungen (Jeweils der 1. des Monats)
    month_days_indices = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]
    for idx, day_idx in enumerate(month_days_indices):
        mx = analemma_x[day_idx]
        my = analemma_y[day_idx]
        ax.plot(mx, my, 'bo', markersize=3)
        ax.text(mx + 4, my, month_names[idx][:3], fontsize=8, color='blue', va='center')

    # Fußpunkt B hervorheben
    ax.plot(0, 0, 'ro', markersize=6, label='Fußpunkt B (Polstab)')
    
    # Layout-Anpassungen
    ax.set_title(f"Vertikal-Südsonnenuhr ({latitude_deg}° N)\nNodus-Höhe: {nodus_height_mm} mm", fontsize=14, pad=20)
    ax.grid(True, linestyle=':', alpha=0.5)
    ax.legend(loc='upper right', fontsize=9)
    plt.axis('off')  # Achsenskalen ausblenden
    
    plt.tight_layout()
    plt.savefig(png_filename, format='png', dpi=300)
    plt.close()
    
    print(f"PNG-Grafik gespeichert unter: {png_filename}")

# Skript ausführen (z. B. für München 48.173333° N, Nodus-Höhe 100 mm)
generate_sundial_png_and_table(latitude_deg=48.173333, nodus_height_mm=100)