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 angle_of_incidence(solar_elevation, solar_azimuth, panel_tilt, panel_azimuth):
    """Calculate angle of incidence on tilted surface"""
    se_rad = math.radians(solar_elevation)
    sa_rad = math.radians(solar_azimuth)
    pt_rad = math.radians(panel_tilt)
    pa_rad = math.radians(panel_azimuth)
    
    cos_incidence = (
        math.sin(se_rad) * math.cos(pt_rad) +
        math.cos(se_rad) * math.sin(pt_rad) * math.cos(sa_rad - pa_rad)
    )
    
    return max(0, cos_incidence)

# 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 ({latitude}°):")
print(f"Solar elevation: {solar_elev:.1f}°")
print(f"Solar azimuth: {solar_az:.1f}°")
print()

# Test úhlů dopadu pro různé panely
panels = [
    ("Vertical N-S", 90, 0),
    ("Vertical E-W", 90, 90), 
    ("Tilt 70°", 70, 0),
    ("Tilt 20°", 20, 0)
]

print("Angle of incidence cosines (higher = better):")
for name, tilt, azimuth in panels:
    cos_inc = angle_of_incidence(solar_elev, solar_az, tilt, azimuth)
    angle_inc = math.degrees(math.acos(cos_inc)) if cos_inc > 0 else 90
    print(f"{name:15}: cos = {cos_inc:.3f}, angle = {angle_inc:.1f}°")