如何计算由tmin/tmax生成的正弦波与tbase/tlimit相交的阴影面积？

一、正弦波温度曲线定义 #

根据需求，单日温度变化的正弦波需满足：

• 起始点x=0时，y=tmin（最低温）
• 峰值点x=12时，y=tmax（最高温）
• 周期24小时（符合单日温度变化规律）

推导得到曲线表达式：

import math

def temp_wave(x, tmin, tmax):
    # x取值范围：0 ≤ x ≤24
    return (tmax + tmin)/2 - (tmax - tmin)/2 * math.cos(math.pi * x / 12)

二、全场景阴影面积计算逻辑 #

阴影面积为正弦波与tbase（下限阈值）、tlimit（上限阈值）围成区域的总面积，即对每个x，计算min(y(x), tlimit) - max(y(x), tbase)在0-24区间的积分（结果为负时取0）。

以下分8种场景处理，同时通过参数钳制避免arccos产生NaN：

场景1：tlimit ≥ tmax 且 tbase ≤ tmin #

正弦波完全落在[tbase, tlimit]区间内，面积为正弦波积分减去tbase的积分：

def case1(tmin, tmax, tbase):
    return 12 * (tmax + tmin) - 24 * tbase

场景2：tlimit ≥ tmax 且 tmin < tbase < tmax #

正弦波部分高于tbase，需先求y(x)=tbase的交点，再计算交点间的积分：

def case2(tmin, tmax, tbase):
    C = (tmax + tmin - 2 * tbase) / (tmax - tmin)
    C = max(min(C, 1.0), -1.0)  # 钳制参数避免NaN
    x1 = (12 / math.pi) * math.acos(C)
    x2 = 24 - x1
    A = (tmax + tmin)/2 - tbase
    sin_term = math.sin(math.pi * x1 / 12)
    term1 = A * (x2 - x1)
    term2 = (12 * (tmax - tmin) / math.pi) * sin_term
    return term1 + term2

场景3：tlimit ≥ tmax 且 tbase ≥ tmax #

正弦波完全低于tbase，面积为0。

场景4：tbase ≤ tmin 且 tmin < tlimit < tmax #

正弦波部分高于tlimit，总面积为场景1面积减去超出tlimit部分的积分：

def case4(tmin, tmax, tbase, tlimit):
    full_area = 12 * (tmax + tmin) - 24 * tbase
    D = (tmax + tmin - 2 * tlimit) / (tmax - tmin)
    D = max(min(D, 1.0), -1.0)
    x1 = (12 / math.pi) * math.acos(D)
    x2 = 24 - x1
    A = (tmax + tmin)/2 - tlimit
    sin_term = math.sin(math.pi * x1 / 12)
    excess_area = A * (x2 - x1) + (12 * (tmax - tmin) / math.pi) * sin_term
    return full_area - excess_area

场景5：tmin < tbase < tlimit < tmax #

需同时求y(x)=tbase和y(x)=tlimit的交点，分三段计算面积：

def case5(tmin, tmax, tbase, tlimit):
    # 求y=tbase的交点
    C = (tmax + tmin - 2*tbase)/(tmax - tmin)
    C = max(min(C, 1.0), -1.0)
    x_a = (12/math.pi)*math.acos(C)
    x_d = 24 - x_a
    # 求y=tlimit的交点
    D = (tmax + tmin - 2*tlimit)/(tmax - tmin)
    D = max(min(D, 1.0), -1.0)
    x_b = (12/math.pi)*math.acos(D)
    x_c = 24 - x_b
    
    # 第一段：x_a到x_b的积分
    A = (tmax+tmin)/2 - tbase
    sin_a = math.sin(math.pi * x_a /12)
    sin_b = math.sin(math.pi * x_b /12)
    integral1 = A*(x_b - x_a) - ((tmax-tmin)/2)*(12/math.pi)*(sin_b - sin_a)
    # 第二段：x_b到x_c的矩形面积
    integral2 = (tlimit - tbase)*(x_c - x_b)
    # 第三段：x_c到x_d的积分（与第一段对称）
    sin_d = math.sin(math.pi * x_d /12)
    sin_c = math.sin(math.pi * x_c /12)
    integral3 = A*(x_d - x_c) - ((tmax-tmin)/2)*(12/math.pi)*(sin_d - sin_c)
    
    return integral1 + integral2 + integral3

场景6：tlimit < tmax 且 tbase ≥ tlimit #

若tbase ≥ tmax，面积为0；否则计算y(x)=tbase交点间的积分（同场景2）。

场景7：tlimit ≤ tmin 且 tbase ≤ tlimit #

正弦波完全高于tlimit，面积为矩形面积：

def case7(tlimit, tbase):
    return 24 * (tlimit - tbase)

场景8：tlimit ≤ tmin 且 tbase > tlimit #

无有效阴影区域，面积为0。

三、整合所有场景的计算函数 #

def calculate_gdd_area(tmin, tmax, tbase, tlimit):
    if tlimit <= tbase:
        return 0.0
    
    # 场景1
    if tlimit >= tmax and tbase <= tmin:
        return 12 * (tmax + tmin) - 24 * tbase
    # 场景2
    elif tlimit >= tmax and tmin < tbase < tmax:
        C = (tmax + tmin - 2 * tbase) / (tmax - tmin)
        C = max(min(C, 1.0), -1.0)
        x1 = (12 / math.pi) * math.acos(C)
        x2 = 24 - x1
        A = (tmax + tmin)/2 - tbase
        sin_term = math.sin(math.pi * x1 / 12)
        term1 = A * (x2 - x1)
        term2 = (12 * (tmax - tmin) / math.pi) * sin_term
        return term1 + term2
    # 场景3
    elif tlimit >= tmax and tbase >= tmax:
        return 0.0
    # 场景4
    elif tbase <= tmin and tmin < tlimit < tmax:
        full_area = 12 * (tmax + tmin) - 24 * tbase
        D = (tmax + tmin - 2 * tlimit) / (tmax - tmin)
        D = max(min(D, 1.0), -1.0)
        x1 = (12 / math.pi) * math.acos(D)
        x2 = 24 - x1
        A = (tmax + tmin)/2 - tlimit
        sin_term = math.sin(math.pi * x1 / 12)
        excess_area = A * (x2 - x1) + (12 * (tmax - tmin) / math.pi) * sin_term
        return full_area - excess_area
    # 场景5
    elif tmin < tbase < tlimit < tmax:
        C = (tmax + tmin - 2*tbase)/(tmax - tmin)
        C = max(min(C, 1.0), -1.0)
        x_a = (12/math.pi)*math.acos(C)
        x_d = 24 - x_a
        
        D = (tmax + tmin - 2*tlimit)/(tmax - tmin)
        D = max(min(D, 1.0), -1.0)
        x_b = (12/math.pi)*math.acos(D)
        x_c = 24 - x_b
        
        A = (tmax+tmin)/2 - tbase
        sin_a = math.sin(math.pi * x_a /12)
        sin_b = math.sin(math.pi * x_b /12)
        integral1 = A*(x_b - x_a) - ((tmax-tmin)/2)*(12/math.pi)*(sin_b - sin_a)
        
        integral2 = (tlimit - tbase)*(x_c - x_b)
        
        sin_d = math.sin(math.pi * x_d /12)
        sin_c = math.sin(math.pi * x_c /12)
        integral3 = A*(x_d - x_c) - ((tmax-tmin)/2)*(12/math.pi)*(sin_d - sin_c)
        
        return integral1 + integral2 + integral3
    # 场景6
    elif tlimit < tmax and tbase >= tlimit:
        if tbase >= tmax:
            return 0.0
        else:
            C = (tmax + tmin - 2 * tbase) / (tmax - tmin)
            C = max(min(C, 1.0), -1.0)
            x1 = (12 / math.pi) * math.acos(C)
            x2 = 24 - x1
            A = (tmax + tmin)/2 - tbase
            sin_term = math.sin(math.pi * x1 / 12)
            term1 = A * (x2 - x1)
            term2 = (12 * (tmax - tmin) / math.pi) * sin_term
            return term1 + term2
    # 场景7
    elif tlimit <= tmin and tbase <= tlimit:
        return 24 * (tlimit - tbase)
    # 场景8
    elif tlimit <= tmin and tbase > tlimit:
        return 0.0
    # 兜底
    return 0.0

四、关键问题解决说明 #

1. NaN错误处理：在计算arccos前，将参数钳制在[-1,1]区间，避免因浮点精度问题导致参数超出定义域。
2. 全场景覆盖：通过分支判断覆盖所有可能的tbase、tlimit与tmin、tmax的大小关系，确保无遗漏场景。
3. 对称简化计算：利用正弦波的对称性（x与24-x处的函数值对称），减少重复计算，提升效率。

内容的提问来源于stack exchange，提问作者X3R0