import pandas as pd
import math

def wind_farm_analysis(device_name, area_km2, file_path, avg_temp, avg_wind_speed,
                       lateral_spacing_factor=5, longitudinal_spacing_factor=10, altitude=11,
                       hub_height=100, Cp=0.45, eta=0.8):
    """
    封装函数：分析风电场的风机数量及各项经济和技术指标，直接输入年平均气温和年平均风速
    参数：
        device_name (str): 风力发电机型号名称
        area_km2 (float): 风电场面积（平方公里）
        file_path (str): 包含风机参数的Excel文件路径
        avg_temp (float): 年平均气温（摄氏度）
        avg_wind_speed (float): 年平均风速（m/s）
        lateral_spacing_factor (float): 横向间距因子（默认为5倍叶片直径，5D）
        longitudinal_spacing_factor (float): 纵向间距因子（默认为10倍叶片直径，10D）
        altitude (float): 海拔高度（m），默认11m
        hub_height (float): 轮毂高度（m），默认100m
        Cp (float): 风能利用系数（功率系数），默认0.45，反映风能转换效率
        eta (float): 总系统效率（包括机械和电气效率），默认0.8
    返回：
        dict: 包含风电场分析结果的字典，包括装机容量、发电量、环境效益等
    """
    def estimate_wind_turbine_count(area_km2, blade_diameter):
        """
        估算风电场可容纳的风机数量，基于面积和风机间距
        参数：
            area_km2 (float): 风电场面积（平方公里）
            blade_diameter (float): 风机叶片直径（m）

        返回：
            int: 估算的风机数量
        """
        # 将面积从平方公里转换为平方米
        area_m2 = area_km2 * 1_000_000
        # 计算横向和纵向间距（以叶片直径为单位）
        lateral_spacing = lateral_spacing_factor * blade_diameter
        longitudinal_spacing = longitudinal_spacing_factor * blade_diameter
        # 单台风机占地面积 = 横向间距 * 纵向间距
        turbine_area = lateral_spacing * longitudinal_spacing
        # 风机数量 = 总面积 / 单台风机占地面积（取整数）
        turbine_count = int(area_m2 / turbine_area)
        print(f"单台风机占地面积: {turbine_area:,} 平方米 "
              f"(横向间距: {lateral_spacing} 米, 纵向间距: {longitudinal_spacing} 米)")
        print(f"估算风机数量: {turbine_count} 台")
        return turbine_count

    def get_wind_turbine_specs(device_name, file_path):
        """
        从Excel文件中获取指定风机的参数
        参数：
            device_name (str): 风机型号名称
            file_path (str): Excel文件路径
        返回：
            tuple: 额定功率（kW）、扫风面积（m²）、叶片直径（m）
        """
        try:
            # 读取Excel文件
            df = pd.read_excel(file_path)
            # 查找匹配的设备名称
            match = df[df.iloc[:, 0] == device_name]
            if not match.empty:
                rated_power = match.iloc[0, 1]  # 额定功率（kW）
                swept_area = match.iloc[0, 7]   # 扫风面积（m²）
                blade_diameter = match.iloc[0, 6]  # 叶片直径（m）
                print(f"找到设备 '{device_name}'，额定功率: {rated_power} KW, "
                      f"扫风面积: {swept_area} m², 叶片直径: {blade_diameter} 米")
                return rated_power, swept_area, blade_diameter
            else:
                raise ValueError(f"未找到设备名称: {device_name}")
        except FileNotFoundError:
            raise FileNotFoundError(f"文件未找到: {file_path}")
        except Exception as e:
            raise Exception(f"发生错误: {str(e)}")

    def air_density(altitude, hub_height, T0):
        """
        计算空气密度，考虑海拔和轮毂高度的影响
        参数：
            altitude (float): 海拔高度（m）
            hub_height (float): 轮毂高度（m）
            T0 (float): 地面平均气温（摄氏度）

        返回：
            float: 空气密度（kg/m³）
        公式：
            ρ = (353.05 / T) * exp(-0.034 * (z / T))
            其中 T = T0 - LR * z + 273.15（T0为地面温度，LR为温度递减率，z为总高度）
        """
        # 计算总高度（海拔 + 轮毂高度）
        z = altitude + hub_height
        # 温度递减率（lapse rate），每升高1米温度降低0.0065°C
        LR = 0.0065
        # 计算绝对温度（K），考虑高度引起的温度变化
        T = T0 - LR * z + 273.15
        # 计算空气密度
        return (353.05 / T) * math.exp(-0.034 * (z / T))

    def wind_power_density(density, velocity_avg):
        """
        计算风功率密度（单位面积的风能功率）
        参数：
            density (float): 空气密度（kg/m³）
            velocity_avg (float): 平均风速（m/s）
        返回：
            float: 风功率密度（W/m²）
        公式：
            P = 0.5 * ρ * v³（按照年平均风速来算）
        """
        return 0.5 * density * velocity_avg**3

    def estimated_wind_power(num_turbines, rated_power):
        """
        计算风电场总装机容量

        参数：
            num_turbines (int): 风机数量
            rated_power (float): 单台风机额定功率（kW）
        返回：
            float: 总装机容量（kW）
        """
        if not isinstance(num_turbines, int) or num_turbines < 0:
            raise ValueError("风机数量必须为非负整数")
        return rated_power * num_turbines

    def calculate_power_output(S, w, Cp, eta, num_turbines):
        """
        计算风电场年发电量

        参数：
            S (float): 扫风面积（m²）
            w (float): 风功率密度（W/m²）
            Cp (float): 风能利用系数
            eta (float): 系统效率
            num_turbines (int): 风机数量
        返回：
            float: 年发电量（Wh）
        公式：
            E = w * S * Cp * 8760 * η * N (N为风机个数）
            其中 8760 为一年小时数
        """
        return w * S * Cp * 8760 * eta * num_turbines

    def calculate_equivalent_hours(P, P_r):
        """
        计算等效满负荷小时数
        参数：
            P (float): 年发电量（Wh）
            P_r (float): 单台风机额定功率（kW）
        返回：
            float: 等效小时数（小时）
        """
        if P_r == 0:
            raise ValueError("额定功率不能为 0")
        return (P / 1000) / P_r

    def calculate_environmental_benefits(E_p_million_kwh):
        """
        计算环境效益（减排量）

        参数：
            E_p_million_kwh (float): 年发电量（万kWh）

        返回：
            dict: 包含标准煤、CO₂、SO₂、NOx减排量的字典

        假设：
            每万kWh可节约标准煤0.404吨，减排CO₂ 0.977吨，SO₂ 0.03吨，NOx 0.015吨
        """
        if E_p_million_kwh < 0:
            raise ValueError("年发电量需≥0")
        return {
            "coal_reduction": E_p_million_kwh * 0.404 * 10,  # kg
            "CO2_reduction": E_p_million_kwh * 0.977 * 10,   # kg
            "SO2_reduction": E_p_million_kwh * 0.03 * 10,    # kg
            "NOX_reduction": E_p_million_kwh * 0.015 * 10    # kg
        }

    # 获取风机参数
    rated_power, swept_area, blade_diameter = get_wind_turbine_specs(device_name, file_path)

    # 估算风机数量
    num_turbines = estimate_wind_turbine_count(area_km2, blade_diameter)

    # 计算空气密度
    avg_density = air_density(altitude, hub_height, avg_temp)

    # 计算风功率密度（W/m²）
    wpd = wind_power_density(avg_density, avg_wind_speed)

    # 计算总装机容量（kW）
    total_power = estimated_wind_power(num_turbines, rated_power)

    # 计算年发电量（Wh）
    P_test = calculate_power_output(swept_area, wpd, Cp, eta, num_turbines)

    # 计算等效满负荷小时数
    h = calculate_equivalent_hours(P_test, rated_power)

    # 转换为万kWh以计算环境效益
    E_p_million_kwh = P_test / 10000000
    env_benefits = calculate_environmental_benefits(E_p_million_kwh)

    # 返回结果字典
    return {
        "device": device_name,
        "rated_power": rated_power,
        "swept_area": swept_area,
        "blade_diameter": blade_diameter,
        "num_turbines": num_turbines,
        "avg_density": avg_density,
        "wpd": wpd,
        "total_power": total_power / 1000,  # 转换为MW
        "annual_power_output": P_test / 10000000,  # 转换为万kWh
        "equivalent_hours": h,
        "coal_reduction": env_benefits["coal_reduction"],
        "CO2_reduction": env_benefits["CO2_reduction"],
        "SO2_reduction": env_benefits["SO2_reduction"],
        "NOX_reduction": env_benefits["NOX_reduction"]
    }

# 主程序
if __name__ == "__main__":
    # 定义输入参数
    file_path = r".\wind_product.xlsx"  # 风机参数文件路径
    device_name = 'GW165-4.0'           # 风机型号
    area_km2 = 10                       # 风电场面积（平方公里）
    avg_temp = 13.0                     # 年平均气温（摄氏度）
    avg_wind_speed = 6                  # 年平均风速（m/s）

    # 调用风电场分析函数
    result = wind_farm_analysis(
        device_name=device_name,
        area_km2=area_km2,
        file_path=file_path,
        avg_temp=avg_temp,
        avg_wind_speed=avg_wind_speed
    )

    # 输出结果
    print(f"\n设备: {result['device']}")
    print(f"额定功率: {result['rated_power']:.2f} KW")
    print(f"扫风面积: {result['swept_area']:.2f} m^2")
    print(f"叶片直径: {result['blade_diameter']:.2f} m")
    print(f"风机数量: {result['num_turbines']} 台")
    print(f"平均空气密度: {result['avg_density']:.3f} kg/m^3")
    print(f"风功率密度: {result['wpd']:.2f} W/m^2")
    print(f"项目装机容量: {result['total_power']:.2f} MW")
    print(f"年发电量: {result['annual_power_output']:.3f} 万 kWh")
    print(f"等效小时数: {result['equivalent_hours']:.2f} 小时")
    print(f"标准煤减排量：{result['coal_reduction']:,.0f} kg")
    print(f"CO₂减排量：{result['CO2_reduction']:,.0f} kg")
    print(f"SO₂减排量：{result['SO2_reduction']:,.0f} kg")
    print(f"NOx减排量：{result['NOX_reduction']:,.0f} kg")