import math

def solar_position(latitude, day_of_year, hour):
    """Calculate solar position for given day and hour"""
    # Solar declination
    declination = 23.45 * math.sin(math.radians(360 * (284 + day_of_year) / 365))
    
    # Hour angle
    hour_angle = 15 * (hour - 12)
    
    # Solar elevation angle
    elevation = math.asin(
        math.sin(math.radians(declination)) * math.sin(math.radians(latitude)) +
        math.cos(math.radians(declination)) * math.cos(math.radians(latitude)) * 
        math.cos(math.radians(hour_angle))
    )
    
    # Solar azimuth angle (standard formula gives azimuth from South)
    azimuth_from_south = math.atan2(
        math.sin(math.radians(hour_angle)),
        math.cos(math.radians(hour_angle)) * math.sin(math.radians(latitude)) - 
        math.tan(math.radians(declination)) * math.cos(math.radians(latitude))
    )
    
    # Convert to azimuth from North (0° = North, 90° = East, 180° = South, 270° = West)
    azimuth = math.degrees(azimuth_from_south) + 180
    if azimuth >= 360:
        azimuth -= 360
    if azimuth < 0:
        azimuth += 360
    
    return math.degrees(elevation), azimuth

def calculate_irradiance_on_tilted_surface(ghi, solar_elevation, solar_azimuth, panel_tilt, panel_azimuth):
    """Calculate irradiance on tilted surface"""
    if solar_elevation <= 0:
        return 0, 0
    
    # Convert to radians
    se_rad = math.radians(solar_elevation)
    sa_rad = math.radians(solar_azimuth)
    pt_rad = math.radians(panel_tilt)
    pa_rad = math.radians(panel_azimuth)
    
    # Angle of incidence
    cos_incidence = (
        math.sin(se_rad) * math.cos(pt_rad) +
        math.cos(se_rad) * math.sin(pt_rad) * math.cos(sa_rad - pa_rad)
    )
    
    if cos_incidence < 0:
        cos_incidence = 0
    
    # Direct normal irradiance (simplified)
    dni = ghi * 0.8 if solar_elevation > 0 else 0
    
    # Direct irradiance on panel
    direct_irradiance = dni * cos_incidence
    
    print(f"  Solar: elev={solar_elevation:.1f}°, azim={solar_azimuth:.1f}°")
    print(f"  Panel: tilt={panel_tilt}°, azim={panel_azimuth}°")
    print(f"  cos_incidence: {cos_incidence:.3f}")
    print(f"  GHI: {ghi:.1f} W/m², DNI: {dni:.1f} W/m²")
    print(f"  Direct on panel: {direct_irradiance:.1f} W/m²")
    
    return direct_irradiance, 0

# Test pro červenec (day 196) v poledne
latitude = -47.3919449
day_of_year = 7 * 30 - 15  # July 15th
hour = 12

solar_elev, solar_az = solar_position(latitude, day_of_year, hour)

print(f"July 15th, noon at Lago Brown:")
print(f"Solar elevation: {solar_elev:.1f}°, azimuth: {solar_az:.1f}°")
print()

# Simulace GHI pro červenec (velmi nízké kvůli zimě)
monthly_ghi = 1.4  # kWh/m²/day
hour_factor = max(0, math.cos(math.pi * abs(hour + 0.5 - 12.5) / 8))
hourly_ghi = monthly_ghi * hour_factor * 1000 / 6  # Convert to W/m²

print(f"July GHI distribution:")
print(f"Monthly GHI: {monthly_ghi} kWh/m²/day")
print(f"Hour factor at noon: {hour_factor:.3f}")
print(f"Hourly GHI at noon: {hourly_ghi:.1f} W/m²")
print()

# Test panelů
panels = [
    ("Vertical N-S", 90, 0),
    ("Vertical E-W", 90, 90), 
    ("Tilt 70°", 70, 0),
    ("Tilt 20°", 20, 0)
]

for name, tilt, azimuth in panels:
    print(f"\n{name}:")
    direct_irr, _ = calculate_irradiance_on_tilted_surface(hourly_ghi, solar_elev, solar_az, tilt, azimuth)
    
    # Power calculation (700W panel)
    power = (direct_irr / 1000) * 700
    print(f"  Power output: {power:.1f} W")