import math

def calculate_full_irradiance(solar_elevation, solar_azimuth, panel_tilt, panel_azimuth):
    """Calculate full irradiance including bifacial effect"""
    
    dni_clear_sky = 900  # W/m²
    bifaciality_factor = 0.80
    ground_albedo = 0.15
    
    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)
    
    # Front side 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)
    )
    
    # Direct irradiance - only front side gets direct sunlight
    front_direct = dni_clear_sky * max(0, cos_incidence)
    back_direct = 0  # Back side never gets direct sunlight
    
    print(f"  Front direct: {front_direct:.1f} W/m² (cos_inc: {cos_incidence:.4f})")
    
    # Diffuse irradiance (sky radiation)
    diffuse_factor = 0.15  # 15% of total radiation is diffuse
    sky_diffuse = dni_clear_sky * diffuse_factor
    
    front_diffuse = sky_diffuse * (1 + math.cos(pt_rad)) / 2
    back_diffuse = sky_diffuse * (1 - math.cos(pt_rad)) / 2
    
    print(f"  Front diffuse: {front_diffuse:.1f} W/m²")
    print(f"  Back diffuse: {back_diffuse:.1f} W/m²")
    
    # Ground reflected irradiance
    ground_reflected = dni_clear_sky * math.sin(se_rad) * ground_albedo
    front_ground = ground_reflected * (1 - math.cos(pt_rad)) / 2
    back_ground = ground_reflected * (1 + math.cos(pt_rad)) / 2
    
    print(f"  Ground reflected: {ground_reflected:.1f} W/m²")
    print(f"  Front ground: {front_ground:.1f} W/m²")
    print(f"  Back ground: {back_ground:.1f} W/m²")
    
    # Bifacial effect - more vertical panels have better bifacial exposure
    bifacial_bonus = math.sin(pt_rad)  # 0.0 for horizontal, 1.0 for vertical
    effective_bifaciality = bifaciality_factor * (0.5 + 0.5 * bifacial_bonus)
    
    print(f"  Bifacial bonus: {bifacial_bonus:.3f}")
    print(f"  Effective bifaciality: {effective_bifaciality:.3f}")
    
    # Total irradiance
    front_irradiance = front_direct + front_diffuse + front_ground
    back_irradiance = (back_direct + back_diffuse + back_ground) * effective_bifaciality
    
    print(f"  Front total: {front_irradiance:.1f} W/m²")
    print(f"  Back total: {back_irradiance:.1f} W/m²")
    
    total_irradiance = front_irradiance + back_irradiance
    power = (total_irradiance / 1000) * 700  # 700W panel
    
    bifacial_contribution = (back_irradiance / total_irradiance * 100) if total_irradiance > 0 else 0
    
    print(f"  Total: {total_irradiance:.1f} W/m²")
    print(f"  Power: {power:.1f} W")
    print(f"  Bifacial contribution: {bifacial_contribution:.1f}%")
    
    return total_irradiance, power

# June 15, 10:00 - elevation 9.7°, azimuth 39°
solar_elev = 9.7
solar_az = 39.0

print(f"June 15, 10:00: Sun at {solar_elev}° elevation, {solar_az}° azimuth")
print()

# Compare 70° vs 90° (vertical) panels
for tilt in [70, 90]:
    print(f"=== Panel tilt {tilt}° (North-facing) ===")
    total_irr, power = calculate_full_irradiance(solar_elev, solar_az, tilt, 0)
    print()