sunangle.py

import math
import datetime
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from timezonefinder import TimezoneFinder
from pytz import timezone
import numpy as np

def get_julian_day(dt):
    """Calculate Julian Day Number"""
    a = (14 - dt.month) // 12
    y = dt.year + 4800 - a
    m = dt.month + 12 * a - 3
    
    jdn = dt.day + (153 * m + 2) // 5 + 365 * y + y // 4 - y // 100 + y // 400 - 32045
    
    # Add fraction of day
    jdn += (dt.hour + dt.minute / 60.0 + dt.second / 3600.0) / 24.0 - 0.5
    
    return jdn

def calculate_sun_position(lat, lon, dt):
    """Calculate sun position (elevation and azimuth) for given location and time"""
    
    # Convert to radians
    lat_rad = math.radians(lat)
    
    # Calculate Julian Day
    jd = get_julian_day(dt)
    
    # Calculate number of days since J2000.0
    n = jd - 2451545.0
    
    # Mean longitude of the Sun
    L = math.radians((280.460 + 0.9856474 * n) % 360)
    
    # Mean anomaly
    g = math.radians((357.528 + 0.9856003 * n) % 360)
    
    # Ecliptic longitude
    lambda_sun = L + math.radians(1.915) * math.sin(g) + math.radians(0.020) * math.sin(2 * g)
    
    # Obliquity of the ecliptic
    epsilon = math.radians(23.439 - 0.0000004 * n)
    
    # Right ascension
    alpha = math.atan2(math.cos(epsilon) * math.sin(lambda_sun), math.cos(lambda_sun))
    
    # Declination
    delta = math.asin(math.sin(epsilon) * math.sin(lambda_sun))
    
    # Greenwich Mean Sidereal Time
    gmst = 18.697374558 + 24.06570982441908 * n
    gmst = gmst % 24
    
    # Local Sidereal Time
    lst = gmst + lon / 15.0
    lst = lst % 24
    
    # Hour angle
    h = math.radians(15 * (lst - alpha * 12 / math.pi))
    
    # Elevation (altitude)
    elevation = math.asin(math.sin(lat_rad) * math.sin(delta) + 
                         math.cos(lat_rad) * math.cos(delta) * math.cos(h))
    
    # Azimuth
    azimuth = math.atan2(math.sin(h), 
                        math.cos(h) * math.sin(lat_rad) - math.tan(delta) * math.cos(lat_rad))
    
    # Convert to degrees
    elevation_deg = math.degrees(elevation)
    azimuth_deg = math.degrees(azimuth) + 180  # Convert to 0-360 range
    
    return elevation_deg, azimuth_deg

def draw_sun_angle(lat, lon, elevation_angle, azimuth_angle, current_time, tz_name):
    """Draw the sun angle visualization"""
    fig, ax = plt.subplots(1, 1, figsize=(10, 8))
    
    # Set up the plot
    ax.set_xlim(-10, 10)
    ax.set_ylim(-2, 10)
    ax.set_aspect('equal')
    
    # Draw ground
    ground = patches.Rectangle((-10, -2), 20, 2, linewidth=1, 
                              edgecolor='black', facecolor='lightgreen')
    ax.add_patch(ground)
    
    # Draw house
    house_base = patches.Rectangle((-1.5, 0), 3, 2, linewidth=2, 
                                  edgecolor='black', facecolor='lightcoral')
    ax.add_patch(house_base)
    
    # Draw roof
    roof = patches.Polygon([[-2, 2], [2, 2], [0, 3.5]], closed=True,
                          linewidth=2, edgecolor='black', facecolor='brown')
    ax.add_patch(roof)
    
    # Draw coordinate marker
    ax.plot(0, 0, 'ro', markersize=10, label=f'Location: ({lat:.2f}°, {lon:.2f}°)')
    
    # Calculate sun position for drawing
    if elevation_angle > 0:  # Sun is above horizon
        # Convert angles to drawing coordinates
        # Use a radius of 7 units for sun distance
        sun_distance = 7
        
        # For visualization, we'll show the sun from a side view
        # X position will be based on azimuth (simplified)
        # Y position will be based on elevation
        
        # Simplified visualization: place sun based on elevation
        sun_x = sun_distance * math.cos(math.radians(elevation_angle))
        sun_y = sun_distance * math.sin(math.radians(elevation_angle))
        
        # Draw sun
        sun = plt.Circle((sun_x, sun_y), 0.5, color='yellow', ec='orange', linewidth=2)
        ax.add_patch(sun)
        
        # Draw sun rays
        for angle in range(0, 360, 30):
            ray_start_x = sun_x + 0.5 * math.cos(math.radians(angle))
            ray_start_y = sun_y + 0.5 * math.sin(math.radians(angle))
            ray_end_x = sun_x + 0.8 * math.cos(math.radians(angle))
            ray_end_y = sun_y + 0.8 * math.sin(math.radians(angle))
            ax.plot([ray_start_x, ray_end_x], [ray_start_y, ray_end_y], 
                   'orange', linewidth=2)
        
        # Draw angle line from coordinate to sun
        ax.plot([0, sun_x], [0, sun_y], 'b--', linewidth=2, alpha=0.7)
        
        # Draw angle arc
        arc_radius = 2
        angles = np.linspace(0, math.radians(elevation_angle), 50)
        arc_x = arc_radius * np.cos(angles)
        arc_y = arc_radius * np.sin(angles)
        ax.plot(arc_x, arc_y, 'b-', linewidth=2)
        
        # Add angle text
        text_x = 2.5 * math.cos(math.radians(elevation_angle/2))
        text_y = 2.5 * math.sin(math.radians(elevation_angle/2))
        ax.text(text_x, text_y, f'{elevation_angle:.1f}°', fontsize=14, 
               ha='left', va='center', color='blue', weight='bold')
        
        # Add sun information
        info_text = f'Sun Elevation: {elevation_angle:.1f}°\nSun Azimuth: {azimuth_angle:.1f}°'
        ax.text(0.02, 0.98, info_text, transform=ax.transAxes, fontsize=12,
               verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
    else:
        # Sun is below horizon
        ax.text(0, 5, 'Sun is below horizon', fontsize=16, ha='center', 
               color='darkblue', weight='bold')
        info_text = f'Sun Elevation: {elevation_angle:.1f}° (below horizon)\nSun Azimuth: {azimuth_angle:.1f}°'
        ax.text(0.02, 0.98, info_text, transform=ax.transAxes, fontsize=12,
               verticalalignment='top', bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))
    
    # Add title and labels
    title = f'Sun Angle at {lat:.2f}°, {lon:.2f}°\n'
    title += f'Time: {current_time.strftime("%Y-%m-%d %H:%M:%S")} {tz_name}'
    ax.set_title(title, fontsize=16, weight='bold')
    
    # Remove axes for cleaner look
    ax.set_xticks([])
    ax.set_yticks([])
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    ax.spines['left'].set_visible(False)
    
    # Add grid for reference
    ax.grid(True, alpha=0.3, linestyle='--')
    
    plt.tight_layout()
    plt.show()

def main():
    print("Sun Angle Calculator and Visualizer")
    print("=" * 40)
    
    # Get user input
    while True:
        try:
            lat_str = input("Enter latitude (e.g., 40.7128 for NYC): ")
            lat = float(lat_str)
            if -90 <= lat <= 90:
                break
            else:
                print("Latitude must be between -90 and 90 degrees.")
        except ValueError:
            print("Please enter a valid number.")
    
    while True:
        try:
            lon_str = input("Enter longitude (e.g., -74.0060 for NYC): ")
            lon = float(lon_str)
            if -180 <= lon <= 180:
                break
            else:
                print("Longitude must be between -180 and 180 degrees.")
        except ValueError:
            print("Please enter a valid number.")
    
    # Find timezone for coordinates
    tf = TimezoneFinder()
    tz_name = tf.timezone_at(lat=lat, lng=lon)
    
    if tz_name is None:
        print(f"Could not determine timezone for coordinates ({lat}, {lon})")
        print("Using UTC as default.")
        tz_name = 'UTC'
    
    # Get current time in that timezone
    tz = timezone(tz_name)
    current_time = datetime.datetime.now(tz)
    
    print(f"\nLocation: {lat:.4f}°, {lon:.4f}°")
    print(f"Timezone: {tz_name}")
    print(f"Current time: {current_time.strftime('%Y-%m-%d %H:%M:%S %Z')}")
    
    # Calculate sun position
    elevation, azimuth = calculate_sun_position(lat, lon, current_time)
    
    print(f"\nSun position:")
    print(f"  Elevation: {elevation:.2f}° {'(above horizon)' if elevation > 0 else '(below horizon)'}")
    print(f"  Azimuth: {azimuth:.2f}° (from North)")
    
    # Create visualization
    print("\nGenerating visualization...")
    draw_sun_angle(lat, lon, elevation, azimuth, current_time, tz_name)

if __name__ == "__main__":
    main()