#!/usr/bin/env python3
"""
SOx Emissions Analysis for Aviation Fuel Sulfur Content
Quantitative comparison of 10 ppm vs 300 ppm vs 3,000 ppm sulfur scenarios
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# Constants and Data Sources
GLOBAL_FUEL_CONSUMPTION_2024 = 299.8e6  # tonnes per year (IATA)
BASE_EMISSION_FACTOR = 1.90e-7  # kg SOx per kg fuel per ppm sulfur (EUROCONTROL)
MORTALITY_RATE_PER_TONNE_SOX = 0.014044  # deaths per tonne SOx
VALUE_OF_STATISTICAL_LIFE = 10e6  # USD per life (EPA)
ENVIRONMENTAL_DAMAGE_COST = 2000  # USD per tonne SOx

# Sulfur content scenarios (ppm)
scenarios = {
    'Current Maximum (3,000 ppm)': 3000,
    'ASTM Minimum (300 ppm)': 300,
    'dLCAF Ultra-Low (10 ppm)': 10
}

def calculate_sox_emissions(fuel_consumption_tonnes, sulfur_ppm):
    """Calculate annual SOx emissions in tonnes"""
    fuel_consumption_kg = fuel_consumption_tonnes * 1000
    sox_emissions_kg = fuel_consumption_kg * BASE_EMISSION_FACTOR * sulfur_ppm
    return sox_emissions_kg / 1000  # Convert back to tonnes

def calculate_health_impacts(sox_emissions_tonnes):
    """Calculate annual deaths from SOx emissions"""
    return sox_emissions_tonnes * MORTALITY_RATE_PER_TONNE_SOX

def calculate_economic_value(deaths_prevented):
    """Calculate economic value of lives saved"""
    return deaths_prevented * VALUE_OF_STATISTICAL_LIFE

def calculate_environmental_damage(sox_emissions_tonnes):
    """Calculate environmental damage costs"""
    return sox_emissions_tonnes * ENVIRONMENTAL_DAMAGE_COST

# Perform calculations for all scenarios
results = {}
for scenario_name, sulfur_ppm in scenarios.items():
    sox_emissions = calculate_sox_emissions(GLOBAL_FUEL_CONSUMPTION_2024, sulfur_ppm)
    deaths = calculate_health_impacts(sox_emissions)
    env_damage = calculate_environmental_damage(sox_emissions)
    
    results[scenario_name] = {
        'sulfur_ppm': sulfur_ppm,
        'sox_emissions_tonnes': sox_emissions,
        'annual_deaths': deaths,
        'environmental_damage_usd': env_damage
    }

# Calculate benefits of switching scenarios
def calculate_benefits(baseline_scenario, improved_scenario):
    baseline = results[baseline_scenario]
    improved = results[improved_scenario]
    
    sox_reduction = baseline['sox_emissions_tonnes'] - improved['sox_emissions_tonnes']
    lives_saved = baseline['annual_deaths'] - improved['annual_deaths']
    health_benefits = calculate_economic_value(lives_saved)
    env_benefits = baseline['environmental_damage_usd'] - improved['environmental_damage_usd']
    total_benefits = health_benefits + env_benefits
    
    return {
        'sox_reduction_tonnes': sox_reduction,
        'sox_reduction_percent': (sox_reduction / baseline['sox_emissions_tonnes']) * 100,
        'lives_saved': lives_saved,
        'health_benefits_usd': health_benefits,
        'environmental_benefits_usd': env_benefits,
        'total_annual_benefits_usd': total_benefits
    }

# Key comparisons
comparisons = {
    '10 ppm vs 3,000 ppm': calculate_benefits('Current Maximum (3,000 ppm)', 'dLCAF Ultra-Low (10 ppm)'),
    '300 ppm vs 3,000 ppm': calculate_benefits('Current Maximum (3,000 ppm)', 'ASTM Minimum (300 ppm)'),
    '10 ppm vs 300 ppm': calculate_benefits('ASTM Minimum (300 ppm)', 'dLCAF Ultra-Low (10 ppm)')
}

# Print results
print("=== SOx EMISSIONS ANALYSIS RESULTS ===\n")

print("SCENARIO RESULTS:")
for scenario, data in results.items():
    print(f"\n{scenario}:")
    print(f"  Sulfur Content: {data['sulfur_ppm']:,} ppm")
    print(f"  Annual SOx Emissions: {data['sox_emissions_tonnes']:,.0f} tonnes")
    print(f"  Annual Deaths: {data['annual_deaths']:,.0f}")
    print(f"  Environmental Damage: ${data['environmental_damage_usd']:,.0f}")

print("\n" + "="*50)
print("BENEFITS ANALYSIS:")

for comparison_name, benefits in comparisons.items():
    print(f"\n{comparison_name}:")
    print(f"  SOx Reduction: {benefits['sox_reduction_tonnes']:,.0f} tonnes ({benefits['sox_reduction_percent']:.1f}%)")
    print(f"  Lives Saved: {benefits['lives_saved']:,.0f} annually")
    print(f"  Health Benefits: ${benefits['health_benefits_usd']:,.0f}")
    print(f"  Environmental Benefits: ${benefits['environmental_benefits_usd']:,.0f}")
    print(f"  Total Annual Benefits: ${benefits['total_annual_benefits_usd']:,.0f}")

# Key findings summary
print("\n" + "="*50)
print("KEY FINDINGS:")
print(f"• dLCAF (10 ppm) reduces SOx emissions by {comparisons['10 ppm vs 3,000 ppm']['sox_reduction_percent']:.1f}%")
print(f"• Annual lives saved with dLCAF: {comparisons['10 ppm vs 3,000 ppm']['lives_saved']:,.0f}")
print(f"• Total annual benefits: ${comparisons['10 ppm vs 3,000 ppm']['total_annual_benefits_usd']/1e9:.1f} billion")
print(f"• Even ASTM minimum (300 ppm) saves {comparisons['300 ppm vs 3,000 ppm']['lives_saved']:,.0f} lives annually")

if __name__ == "__main__":
    # Create visualization
    plt.style.use('seaborn-v0_8')
    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
    
    # SOx Emissions Chart
    scenarios_list = list(scenarios.keys())
    emissions_data = [results[s]['sox_emissions_tonnes'] for s in scenarios_list]
    colors = ['#e74c3c', '#3498db', '#27ae60']
    
    ax1.bar(scenarios_list, emissions_data, color=colors, alpha=0.8)
    ax1.set_title('Annual SOx Emissions by Scenario', fontsize=14, fontweight='bold')
    ax1.set_ylabel('SOx Emissions (tonnes/year)')
    ax1.tick_params(axis='x', rotation=45)
    
    # Health Impact Chart
    deaths_data = [results[s]['annual_deaths'] for s in scenarios_list]
    ax2.bar(scenarios_list, deaths_data, color=colors, alpha=0.8)
    ax2.set_title('Annual Deaths by Scenario', fontsize=14, fontweight='bold')
    ax2.set_ylabel('Annual Deaths')
    ax2.tick_params(axis='x', rotation=45)
    
    # Benefits Comparison
    comparison_names = list(comparisons.keys())
    benefits_data = [comparisons[c]['total_annual_benefits_usd']/1e9 for c in comparison_names]
    ax3.bar(comparison_names, benefits_data, color=['#27ae60', '#3498db', '#f39c12'], alpha=0.8)
    ax3.set_title('Annual Benefits by Comparison', fontsize=14, fontweight='bold')
    ax3.set_ylabel('Annual Benefits (Billion USD)')
    ax3.tick_params(axis='x', rotation=45)
    
    # SOx Reduction Potential
    reduction_data = [comparisons[c]['sox_reduction_percent'] for c in comparison_names]
    ax4.bar(comparison_names, reduction_data, color=['#27ae60', '#3498db', '#f39c12'], alpha=0.8)
    ax4.set_title('SOx Emission Reduction Potential', fontsize=14, fontweight='bold')
    ax4.set_ylabel('Reduction (%)')
    ax4.tick_params(axis='x', rotation=45)
    
    plt.tight_layout()
    plt.savefig('sox_analysis_results.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    print(f"\nVisualization saved as 'sox_analysis_results.png'")