###  Y/M/D 2025/02/28_Beta_hqw
###  This function is to calculate Si module of soil water balance.   # 根系吸水量。

import datetime
import json
import math
from dataclasses import dataclass
from typing import Any, Dict, List, Tuple

import numpy as np
import pandas as pd
import uvicorn
from fastapi import FastAPI, File, HTTPException, UploadFile
from pydantic import BaseModel


class DateUtils:
    @staticmethod
    def cal_Doy(date: str) -> int:
        """计算日序"""
        try:
            month, day, year = map(int, date.split('/'))
            date_obj = datetime.date(year, month, day)
            return int(date_obj.strftime('%j'))
        except (ValueError, IndexError):
            raise ValueError(f"日期格式错误: {date}")

    @staticmethod
    def cal_DayLength(date: str, latitude: float) -> float:
        """计算日长"""
        day_of_year = DateUtils.cal_Doy(date)
        RAD = math.pi / 180.0
        DEC = -math.asin(math.sin(23.45 * RAD) * math.cos(2.0 * math.pi * (day_of_year + 10.0) / 365.0))
        dSINLD = math.sin(RAD * latitude) * math.sin(DEC)
        dCOSLD = math.cos(RAD * latitude) * math.cos(DEC)
        dCOSLD = dCOSLD if dCOSLD != 0.0 else 0.001
        SSCC = dSINLD / dCOSLD
        return 12.0 * (1.0 + 2.0 * math.asin(SSCC) / math.pi)

class TemperatureUtils:
    @staticmethod
    def get_hourly_temps(min_temperature: float, max_temperature: float) -> list:
        """获取逐时温度"""
        if min_temperature > max_temperature:
            raise ValueError("温度参数错误")
        return [(max_temperature + min_temperature)/2 +
                (max_temperature - min_temperature)*math.cos(0.2618*(i-14))/2
                for i in range(1, 25)]

    @staticmethod
    def get_day_temps(min_temperature: float, max_temperature: float) -> float:
        """获取白天平均温度"""
        if min_temperature > max_temperature:
            raise ValueError("温度参数错误")
        else:
            T5 = []
            for i in range(1, 9):
                # print(i)
                if i == 4 or i == 5 or i == 6:
                    continue
                else:
                    Tfac = 0.9310 + 0.1140 * i - 0.0703 * i ** 2 + 0.0053 * i ** 3
                    # print(format(min_temperature + Tfac * (max_temperature - min_temperature), '.4f'))
                    T5.append(format(min_temperature + Tfac * (max_temperature - min_temperature), '.4f'))
                    T5 = [float(x) for x in T5]
            average_Tday = np.mean(T5)
        return average_Tday

    @staticmethod
    def dailyGDD(Tmax: float, Tmin: float, Tb: float, Tm: float) -> float:
        """计算日有效积温"""
        Tmean = (Tmax + Tmin) / 2
        Tmean = max(min(Tmean, Tm), Tb)
        return Tmean - Tb

class VernalizationCalculator:
    @staticmethod
    def cal_DVE(PVT: float, VD: float, min_temp: float, max_temp: float) -> float:
        """春化效应计算"""
        T24 = TemperatureUtils.get_hourly_temps(min_temp, max_temp)
        Tbv = -PVT / 12
        Tol = 5 - PVT / 12
        Tou = 10 - PVT / 12
        Tmv = 20 - PVT / 6
        vef = 1 / (2 - 0.0167 * PVT)

        sum_RVE = sum_DDVE = 0
        for Ti in T24:
            # 春化效应计算
            if Ti <= Tbv or Ti >= Tmv:
                RVEi = 0
            elif Tbv < Ti < Tol:
                RVEi = math.sin((Ti-Tbv)/(Tol-Tbv)*math.pi/2)**0.5
            elif Tol <= Ti <= Tou:
                RVEi = 1
            else:
                RVEi = math.sin((Ti-Tou)/(Tmv-Tou)*math.pi/2)**vef
            sum_RVE += RVEi

            # 反春化效应计算
            sum_DDVE += 0.5*(Ti-27) if Ti > 27 else 0

        DVE = sum_RVE / 24
        DDVE = sum_DDVE / 24

        if VD < 0.3*PVT:
            return DVE - DDVE
        elif 0.3*PVT <= VD <= PVT:
            return DVE
        return 0


class PhotoperiodCalculator:
    @staticmethod
    def cal_DPE(date: str, latitude: float, P: float, PS: float) -> float:
        """光周期效应计算"""
        DL = DateUtils.cal_DayLength(date, latitude)
        return max(0, 2/(1 + math.exp(PS*(P - DL))))

class StressCalculator:
    def __init__(self, plant_params: dict):
        self.plant_params = plant_params

    def _init_parameters(self):
        """参数初始化"""
        # 物候阶段阈值
        self.PDT_EM = self.plant_params['1.1']
        self.PDT_BT = self.plant_params['1.5']
        self.PDT_SF = self.plant_params['2.2']
        self.PDT_EF = self.plant_params['2.4']
        self.PDT_MT = self.plant_params['3.0']


    @staticmethod
    def temperature_stress(Tmin: float, Tmax: float) -> float:
        """温度胁迫"""
        Tmn = 0
        Tpl = 20
        Tpu = 25
        Tmx = 35

        Tmean = (Tmax + Tmin) / 2
        if Tmn <= Tmean <= Tpl:
            TE = math.sin((Tmean - Tmn) / (Tpl - Tmn) * math.pi / 2)
        elif Tpl < Tmean <= Tpu:
            TE = 1
        elif Tpu < Tmean <= Tmx:
            TE = math.sin((Tmx - Tmean) / (Tmx - Tpu) * math.pi / 2)
        else:
            TE = 0

        return TE


    def age_stress(self, current_pdt: float):
        """生理年龄胁迫"""
        # 计算生理年龄胁迫因子
        if current_pdt <= self.PDT_SF:
            FAL = 1.0
        else:
            FAL = 1 - 0.013 * (current_pdt - self.PDT_SF)

        if current_pdt <= self.PDT_EF:
            FAP = 1.0
        else:
            FAP = 1 - 1 / 48 * (current_pdt - self.PDT_EF)

        return FAL, FAP

    @staticmethod
    def co2_stress(cx: float) -> float:
        """二氧化碳浓度胁迫
        :param Cx: 实际CO2浓度，ppm
        :return:
            FCO2, CO2影响因子
        """
        FCO2 = 1 + 0.8 * math.log(cx / 340)

        return FCO2

    @staticmethod
    def NNI_stress(TOPWT: float, NA: float) -> float:
        """N浓度胁迫
        :param TOPWT: 地上部干物质量，kg/ha
        :param NA: 地上部组织中的实际氮浓度, %
        :return:
            NNI, 氮素营养指数
        """
        DM = TOPWT / 1000  # 地上部干重单位换算成 t/ha
        if DM >= 0.88:
            NC = 4.48 * DM ** (-0.25)
        else:
            NC = 4.63

        NNI = NA / NC

        return NNI

    @staticmethod
    def WDF_stress(ETcadj: float, ETo: float) -> float:
        """
        :param ETcadj: 实际蒸散量，mm
        :param ETo: 潜在蒸散量，mm
        :return:
            WDF，水分亏缺因子"""
        WDF = ETcadj / ETo

        return WDF

class FlowerAreaIndex:
    def __init__(self, DW_onset_flower: float, FlowerArea: float = 1.99):
        """
        花面积指数(FAI)计算模型

        :param DW_onset_flower : float - 初花期地上部干物重 (kg hm⁻²)
        :param FlowerArea : float, optional - 单花面积 (cm²)，默认1.99
        """
        # 品种参数
        self.FlowerArea = FlowerArea

        # 计算开花速率
        self.Flowerrate = 41.5 * DW_onset_flower

        # 状态变量初始化
        self.FlowerNum = 0.0  # 当前花朵总数
        self.flower_queue = []  # 花朵开放队列 (剩余积温, 花朵数量)


    def _calc_senescence(self, Tmean: float) -> float:
        """计算当日凋谢花朵数"""
        senesced = 0.0
        remaining_flowers = []

        for batch in self.flower_queue:
            batch["remaining"] -= Tmean
            if batch["remaining"] <= 0:
                senesced += batch["amount"]
            else:
                remaining_flowers.append(batch)

        self.flower_queue = remaining_flowers
        return senesced


    def daily_step(self, Tmean: float) -> float:
        """
        执行单日计算

        :param Tmean : float - 当日平均温度(℃)
        :param date : str, optional - 日期字符串,月/日/年

        """
        # 当日新开花朵数
        new_flowers = self.Flowerrate * Tmean

        # 将新花朵加入追踪队列
        self.flower_queue.append({"remaining": 100.0, "amount": new_flowers})

        # 计算当日凋谢花朵数
        senesced = self._calc_senescence(Tmean)
        # 更新花朵总数
        self.FlowerNum += new_flowers - senesced

        # 计算花面积指数
        FAI = self.FlowerNum * self.FlowerArea

        return FAI


@dataclass
class GaussianCoefficients:
    """高斯积分系数容器"""
    nodes: Tuple[float, ...]
    weights: Tuple[float, ...]

class PhotosynthesisSimulator:
    # 高斯积分系数
    TIME_GAUSS = GaussianCoefficients(
        nodes=(0.112702, 0.500000, 0.887298),
        weights=(0.277778, 0.444444, 0.277778)
    )

    LAYER_GAUSS = GaussianCoefficients(
        nodes=(0.0469101, 0.2307534, 0.5000000, 0.7691465, 0.9530899),
        weights=(0.1184635, 0.2393144, 0.2844444, 0.2393144, 0.1184635)
    )
    def __init__(self, latitude: float, SH: float, PDT_BT: float, PDT_EF: float):
        self.latitude = latitude
        self.SH = SH
        self.PDT_BT = PDT_BT
        self.PDT_EF = PDT_EF

    def _cal_PARCAN(self, current_date: datetime) -> float:
        """计算冠层光合有效辐射PARCAN"""
        RAD = math.pi / 180.0
        Gsc = 1367

        # 日地距离系数计算
        dr = 1 + 0.033 * math.cos(2 * math.pi * DateUtils.cal_Doy(current_date) / 365)

        # 太阳赤纬计算
        delta = 0.409 * math.sin(2 * math.pi * DateUtils.cal_Doy(current_date) / 365 - 1.39)

        # 太阳时角计算
        ws = math.acos(-math.tan(self.latitude * RAD) * math.tan(delta))

        # 天文辐射计算
        Q = (Gsc / math.pi * dr *
             (ws * math.sin(self.latitude * RAD) * math.sin(delta) +
              math.sin(ws) * math.cos(self.latitude * RAD) * math.cos(delta)))

        # 日长计算
        DL = DateUtils.cal_DayLength(current_date, self.latitude)

        # PAR计算
        PARCAN = (0.23 + 0.48 * self.SH / DL) * Q
        return PARCAN

    def _cal_DTGA(self, current_date: datetime, GAI: float, PAR: float,
                  kx: float, px: float, PLMX: float) -> float:
        """
        计算每日光合同化物 (DTGA)

        :param current_date: 当前日期
        :param GAI: 绿色面积指数 (叶/角果面积指数)
        :param PAR: 光合有效辐射
        :param kx: 消光系数
        :param px: 反射率
        :param PLMX: 最大光合速率
        """
        RAD = math.pi / 180.0

        # 计算太阳赤纬
        delta = 0.409 * math.sin(2 * math.pi * DateUtils.cal_Doy(current_date) / 365 - 1.39)

        # 太阳位置参数
        SSIN = math.sin(self.latitude * RAD) * math.sin(delta)
        CCOS = math.cos(self.latitude * RAD) * math.cos(delta)

        # 日长计算
        DL = DateUtils.cal_DayLength(current_date, self.latitude)

        total_assimilation = 0.0

        # 时间维度高斯积分（每日时角）
        for t_node, t_weight in zip(self.TIME_GAUSS.nodes, self.TIME_GAUSS.weights):
            # 将积分点转换为实际时间
            hour = 12 + 0.5 * DL * t_node
            solar_angle = 2 * math.pi * (hour + 12) / 24
            # 计算瞬时太阳高度角正弦
            sin_beta = SSIN + CCOS * math.cos(solar_angle)
            if sin_beta < 0:
                continue  # 夜间时段

            # 冠层维度高斯积分（垂直分布）
            layer_integral = 0.0
            for l_node, l_weight in zip(self.LAYER_GAUSS.nodes, self.LAYER_GAUSS.weights):
                cumulative_gai = GAI * l_node

                # 光强衰减计算（比尔-朗伯定律）
                light_intensity = (1 - px) * PAR * 3 * sin_beta * math.exp(-kx * cumulative_gai)     # PAR偏低，加上系数

                # 光合响应模型（非直角双曲线）
                photosynthesis = PLMX * (1 - math.exp(-0.45 * light_intensity / PLMX))

                layer_integral += photosynthesis * l_weight

            # 积分结果累加
            total_assimilation += layer_integral * t_weight
        DTGA = total_assimilation * DL
        return DTGA

    def Respiratory(self, FDTGA: float, Tmean: float, Tday: float) -> tuple:
        """呼吸作用计算"""
        RmTo = 0.15
        Q10 = 2
        To = 25
        RM = RmTo * FDTGA * Q10 ** ((Tmean - To) / 10)
        RG = 0.39 * FDTGA
        RpTo = 0.33
        RP = FDTGA * RpTo * Q10 ** ((Tday - To) / 10)
        return (RM, RG, RP)

    def _cal_distribution(self, current_pdt: float, biomass: float) -> tuple:
        """干物质分配计算"""
        if current_pdt < self.PDT_BT:
            PIR = (0.18 - 0.0386) * math.sin(current_pdt / self.PDT_BT * math.pi / 2) + 0.0386
            PIFP = 0
            PIF = 0
            PIL = 0.8 if current_pdt == 0 else -0.0600 * math.log(current_pdt / self.PDT_BT) + 0.22
        else:
            PIR = -0.0011 * (current_pdt - self.PDT_BT) + 0.18
            PIFP = 0.54 * 0.1202 * (math.exp(0.0439 * (current_pdt - self.PDT_BT)) - 1)
            PIL = -0.1454 * math.log(current_pdt / self.PDT_BT) + 0.22
            PIF = 0.05 * math.exp(-(current_pdt - self.PDT_EF) ** 2 / (9.486 ** 2)) if current_pdt <= self.PDT_EF else 0

        PIS = 1 - PIR
        PIP = PIFP - PIF
        PIST = PIS - PIL - PIFP

        ROOTWT = biomass * PIR
        TOPWT = biomass * PIS
        WLV = TOPWT * PIL
        WST = TOPWT * PIST
        WAPV = TOPWT * PIP
        WFL = TOPWT * PIF

        return (ROOTWT, TOPWT, WLV, WAPV, WFL, WST)


class BiomassModel:
    def __init__(self, start: str, end: str, weather: pd.DataFrame, cultivar_params: dict,
                 management: dict, plant_params: dict, soil_params: dict):
        self.start = start
        self.end = end
        self.weather = weather
        self.cultivar_params = cultivar_params
        self.management = management
        self.plant_params = plant_params
        self.soil_params = soil_params

        self._parse_dates()
        self._init_parameters()

    def _parse_dates(self):
        """日期解析"""
        def parse_date(date_str):
            return datetime.date(*map(int, (date_str.split('/')[-1],
                                           date_str.split('/')[0],
                                           date_str.split('/')[1])))

        self.sow_date = parse_date(self.management['sowdate'])
        self.start_date = parse_date(self.start)
        self.end_date = parse_date(self.end)
        self.latitude = self.management['latitude']
        self.seed_weight = self.management['seed_weight']

    def _init_parameters(self):
        """参数初始化"""
        self.P = self.cultivar_params['P']
        self.PS = self.cultivar_params['PS']
        self.TS = self.cultivar_params['TS']
        self.PVT = self.cultivar_params['PVT']
        self.FDF = self.cultivar_params['FDF']
        self.PLMXp = self.cultivar_params['PLMXp']
        self.PLMXl = self.cultivar_params['PLMXl']
        self.SLA = self.cultivar_params['SLA']
        self.SPA = self.cultivar_params['SPA']
        self.pp = self.plant_params['pp']
        self.pl = self.plant_params['pl']
        self.kp = self.plant_params['kp']
        self.kl = self.plant_params['kl']
        self.RTGR = self.plant_params['RTGR']
        self.RFAC1 = self.plant_params['RFAC1']

        # 物候阶段阈值
        self.GDD_emergence = self.plant_params.get('GDD_emergence', 120)  # 默认120℃·d
        self.emergence_calculated = False
        self.PDT_EM = self.plant_params['1.1']
        self.PDT_BT = self.plant_params['1.5']
        self.PDT_SF = self.plant_params['2.2']
        self.PDT_EF = self.plant_params['2.4']
        self.PDT_MT = self.plant_params['3.0']

        # 温度参数
        self.Tb = self.cultivar_params['Tb']
        self.To = self.cultivar_params['To']
        self.Tm = self.cultivar_params['Tm']

        # 土壤参数
        self.SNO3_0 = self.soil_params['SNO3_0']
        self.SNH4_0 = self.soil_params['SNH4_0']
        self.BD = self.soil_params['BD']
        self.DLAYR = self.soil_params['DLAYR']
        self.NLAYR = self.soil_params['NLAYR']
        self.DS = self.soil_params['DS']
        self.DUL = self.soil_params['DUL']
        self.SAT = self.soil_params['SAT']
        self.LL = self.soil_params['LL']
        self.PH = self.soil_params['PH']
        self.SOILC = self.soil_params['SOILC']
        self.FOM = self.soil_params['FOM']
        self.FON = self.soil_params['FON']
        self.TOTN = self.soil_params['TOTN']
        self.SW0 = self.soil_params['SW0']
        self.CLAY = self.soil_params['CLAY']
        self.CN = self.soil_params['CN']

    def _get_temperature_parameters(self, current_pdt: float) -> tuple:
        """获取当前阶段的温度参数"""
        if current_pdt <= self.PDT_BT:
            return self.Tb[1], self.To[1], self.Tm[1]
        elif self.PDT_BT < current_pdt <= self.PDT_SF:
            return self.Tb[2], self.To[2], self.Tm[2]
        elif self.PDT_SF < current_pdt <= self.PDT_EF:
            return self.Tb[3], self.To[3], self.Tm[3]
        elif self.PDT_EF < current_pdt <= self.PDT_MT:
            return self.Tb[4], self.To[4], self.Tm[4]
        return self.Tb[-1], self.To[-1], self.Tm[-1]

    def _calculate_daily_effects(self, current_date: str, current_pdt: float, VD: float):
        """计算每日效应值"""
        # 光周期效应
        DPE = PhotoperiodCalculator.cal_DPE(current_date, self.latitude, self.P, self.PS)

        # 温度效应
        Tb, To, Tm = self._get_temperature_parameters(current_pdt)
        Tmin = self.weather.loc[current_date, 'Tmin']
        Tmax = self.weather.loc[current_date, 'Tmax']
        DTE = self._calculate_DTE(Tb, To, Tm, Tmin, Tmax)

        # 春化效应
        DVE = VernalizationCalculator.cal_DVE(self.PVT, VD, Tmin, Tmax)

        return DPE, DTE, DVE

    def _calculate_DTE(self, Tb: float, To: float, Tm: float, Tmin: float, Tmax: float) -> float:
        """温度效应计算"""
        T24 = TemperatureUtils.get_hourly_temps(Tmin, Tmax)
        total_TE = 0
        for Ti in T24:
            if Ti < Tb or Ti > Tm:
                TE = 0
            elif Tb <= Ti <= To:
                TE = (math.sin((Ti-Tb)/(To-Tb)*math.pi/2)) ** self.TS
            else:
                TE = (math.sin((Tm-Ti)/(Tm-To)*math.pi/2)) ** ((Tm-To)/(To-Tb)*self.TS)
            total_TE += TE
        return total_TE / 24

    def ifEmergence(self) -> str:
        """判断出苗期开始时间"""

        GDD = 0
        current_date = self.sow_date
        while GDD < self.GDD_emergence:
            date_str = current_date.strftime("%m/%d/%Y")
            try:
                Tmax = self.weather.loc[date_str, 'Tmax']
                Tmin = self.weather.loc[date_str, 'Tmin']
            except KeyError:
                current_date += datetime.timedelta(days=1)
                continue

            # 判断是否出苗
            GDD_daily = TemperatureUtils.dailyGDD(Tmax, Tmin,
                                                  Tb=self.cultivar_params['Tb'][0],Tm=self.cultivar_params['Tm'][0])
            GDD += GDD_daily

            # 记录出苗日期
            if GDD >= self.GDD_emergence:
                self.emergence_date = current_date
                self.emergence_days = (current_date - self.sow_date).days
                break

            current_date += datetime.timedelta(days=1)

            # 防止无限循环
            if (current_date - self.sow_date).days > 365:
                self.emergence_date = None
                self.emergence_days = 9999
                break

    def _cal_pdt(self, current_date: datetime.date) -> float:
        """执行模拟"""

        # 阶段1：计算出苗期
        self.ifEmergence()

        # 阶段2：出苗后发育阶段
        VD = 0  # 累积春化天数
        PDT = self.PDT_EM  # 从出苗开始发育

        # 模拟循环
        delta = (current_date - self.emergence_date).days
        for day in (self.emergence_date + datetime.timedelta(n) for n in range(delta + 1)):
            date_str = day.strftime("%m/%d/%Y")
            # 计算每日效应
            DPE, DTE, DVE = self._calculate_daily_effects(date_str, PDT, VD)
            VD += DVE

            # 计算物候发育进度
            if PDT <= self.PDT_BT:
                DPDT = DTE * (DPE * min(VD / self.PVT, 1) if self.PVT else DPE)
            else:
                DPDT = DTE * self.FDF if PDT > self.PDT_SF else DTE * DPE

            DPDT = max(0, DPDT)

            PDT += DPDT
            if PDT > self.PDT_MT:
                break

        return PDT


    def daily_DRTDEP(self, current_date: str, current_pdt: float, wdf):
        """计算每日增加的根系深度"""

        Tb, To, Tm = self._get_temperature_parameters(current_pdt)
        Tmin = self.weather.loc[current_date, 'Tmin']
        Tmax = self.weather.loc[current_date, 'Tmax']

        GDD = TemperatureUtils.dailyGDD(Tmax, Tmin, Tb, Tm)
        DRTDEP = self.RTGR * GDD * wdf  # DRTDEP(cm)

        return DRTDEP

    def RLV(self, RTDEP: float, rootwt: float):
        """
        计算单位体积土壤中根系的总长度，单位为 cm/cm3

        :param RTDEP: 根深，单位为 cm
        :param rootwt: 地下部生物量(干重)，单位为 kg 干物质/ha/d
        :return:
            RLV: Root Length Volume, 指单位体积土壤中根系的总长度，单位为 cm/cm3
        """
        CUMDEP = 0
        RLV_list = []  # cm/cm3
        for L in range(self.NLAYR):
            DEP = min(RTDEP - CUMDEP, self.DLAYR[L])
            RLINIT = rootwt / 10000 / 10000 * 1000 * self.RFAC1 * DEP / (RTDEP)
            CUMDEP += DEP
            RLV = round(RLINIT / self.DLAYR[L], 3)
            RLV_list.append(RLV)

            if CUMDEP > RTDEP:
                break

        return RLV_list

    def biomass(self, day: datetime.date, Biomass0: float, TE: int = 0, FA: int = 0, FCO2: int = 0, NNI: int = 0,
                WDF: int = 0) -> tuple:
        """
        执行模拟

        :param TE: 温度胁迫， 0 - 无胁迫（默认），1 - 有胁迫
        :param FA: 生长年龄胁迫， 0 - 无胁迫（默认），1 - 有胁迫
        :param FCO2: CO2胁迫， 0 - 无胁迫（默认），1 - 有胁迫
        :param NNI: 氮素胁迫， 0 - 无胁迫（默认），1 - 有胁迫
        :param WDF: 水分胁迫， 0 - 无胁迫（默认），1 - 有胁迫
        """

        # 计算出苗期，从出苗后开始发育
        self.ifEmergence()

        date_str = day.strftime("%m/%d/%Y")
        try:
            Tmax = self.weather.loc[date_str, 'Tmax']
            Tmin = self.weather.loc[date_str, 'Tmin']
            Tmean = (Tmax + Tmin) / 2
            SH = self.weather.loc[date_str, 'SH']
        except KeyError:
            print(f"气象数据中缺少 {date_str} 的记录，跳过该天。")

        PDT = self._cal_pdt(day)

        # 开始计算

        photo = PhotosynthesisSimulator(self.latitude, SH, self.PDT_BT, self.PDT_EF)

        # 地上部、地下部、叶片、角果、花、茎生物量
        ROOTWT, TOPWT, WLV, WAPV, WFL, WST = photo._cal_distribution(PDT, Biomass0)

        # 胁迫判断
        stess = StressCalculator(self.plant_params)
        if FA == 0:
            FAL, FAP = 1.0, 1.0
        else:
            FAL, FAP = stess.age_stress(PDT)

        if TE == 0:
            tem_stress = 1.0
        else:
            tem_stress = StressCalculator.temperature_stress(Tmin, Tmax)

        if FCO2 == 0:
            co2 = 1.0
        else:
            co2 = StressCalculator.co2_stress(cx=340)

        if WDF == 0:
            wdf = 1.0
        # else:
        #     wdf = StressCalculator.WDF_stress(ETcadj, ETo)

        if NNI == 0:
            nni = 1.0
        # else:
        #     nni = StressCalculator.NNI_stress(TOPWT, NA)

        PLMXL = self.PLMXl * co2 * FAL * tem_stress * min(wdf, nni)
        PLMXP = self.PLMXp * co2 * FAP * tem_stress * min(wdf, nni)

        # 花面积指数
        startflowering = False
        if not startflowering and PDT >= self.PDT_SF:
            OnesetFlower = WFL
            startflowering = True

        if self.PDT_SF <= PDT <= self.PDT_EF:
            flowermodel = FlowerAreaIndex(DW_onset_flower=OnesetFlower)
            FAI = flowermodel.daily_step(Tmean)
        else:
            FAI = 0

        # 叶面积指数
        if PDT <= self.PDT_BT:
            SLAv = self.SLA
        else:
            SLAv = self.SLA * (-0.0009 * (PDT - self.PDT_BT) ** 2 + 0.0683 * (PDT - self.PDT_BT) + 1.177)

        LAI = SLAv * WLV

        # 角果面积指数
        if PDT >= self.PDT_SF:
            SPAv = self.SPA * (0.00017 * (PDT - self.PDT_SF) ** 2 - 0.05554 * (PDT - self.PDT_SF) + 2.91887)
            PAI = SPAv * WAPV
        else:
            PAI = 0

        # 计算光合同化物
        if PDT < self.PDT_SF:
            Tp = photo._cal_PARCAN(date_str)  # 无角果遮挡
        else:
            Tp = (1 - self.pp) * photo._cal_PARCAN(date_str) * math.exp(-self.kp * PAI)  # 有角果遮挡

        Tl = (1 - self.pl) * Tp * math.exp(-self.kl * LAI)  # radiation passing through the leaves
        Rl = 0.05 * Tp  # leaf reflection
        Rs = 0.2 * Tl  # soil reflection

        if PDT >= self.PDT_SF:
            PARP = (1 - self.pp) * photo._cal_PARCAN(date_str) + Rl
        else:
            PARP = 0
        PARL = (1 - self.pl) * Tp + Rs

        DTGAp = photo._cal_DTGA(date_str, PAI, PARP, self.kp, self.pp, PLMXP)
        DTGAl = photo._cal_DTGA(date_str, LAI, PARL, self.kl, self.pl, PLMXL)
        # print('DTGAp, DTGAl: ', DTGAp, DTGAl)

        FDTGA = (DTGAl + DTGAp) * min(wdf, nni)
        # print('FDTGA: ', FDTGA)

        # RmTo = 0.03*WLV+0.015*WST+0.015*ROOTWT+0.01*WAPV
        # print('RmTo: ', RmTo)
        Tday = TemperatureUtils.get_day_temps(Tmin, Tmax)
        RM, RG, RP = photo.Respiratory(FDTGA, Tmean, Tday)
        # print('RM, RG, RP: ', RM, RG, RP)

        PND = FDTGA - RM - RG
        # print('PND: ', PND)

        TDRW = 30 / 44 * 0.95 * PND / (1 - 0.05)

        TBiomass = TDRW + Biomass0

        return TBiomass, LAI, TOPWT, ROOTWT


class SoilWaterModel:
    def __init__(self, start: str,
                 end: str,
                 weather: pd.DataFrame,
                 cultivar_params: dict,
                 plant_params: dict,
                 soil_params: dict,
                 management: dict,
                 n_params: dict):
        self.start = start
        self.end = end
        self.weather = weather
        self.cultivar_params = cultivar_params
        self.plant_params = plant_params
        self.soil_params = soil_params
        self.management = management
        self.n_params = n_params

        self._parse_dates()
        self._init_parameters()

    def _parse_dates(self):
        """日期解析"""
        def parse_date(date_str):
            return datetime.date(*map(int, (date_str.split('/')[-1],
                                           date_str.split('/')[0],
                                           date_str.split('/')[1])))

        self.sow_date = parse_date(self.management['sowdate'])
        self.start_date = parse_date(self.start)
        self.end_date = parse_date(self.end)

    def _init_parameters(self):
        """参数初始化"""

        # 管理措施参数
        self.latitude = self.management['latitude']
        self.irrigation = self.management['irrigation']
        self.elevation = self.management['elevation']
        self.seed_weight = self.management['seed_weight']

        # 作物参数
        self.GDD_emergence = self.plant_params.get('GDD_emergence', 120)  # 默认120℃·d
        self.emergence_calculated = False
        self.PDT_EM = self.plant_params['1.1']
        self.PDT_BT = self.plant_params['1.5']
        self.PDT_SF = self.plant_params['2.2']
        self.PDT_EF = self.plant_params['2.4']
        self.PDT_MT = self.plant_params['3.0']
        self.RTGR = self.plant_params['RTGR']
        self.Kcbini = self.plant_params['Kcbini']
        self.Kcbmid = self.plant_params['Kcbmid']
        self.Kcbend = self.plant_params['Kcbend']
        self.h = self.plant_params['h']
        self.kl = self.plant_params['kl']

        # 品种参数
        self.Po = self.cultivar_params['P']
        self.PS = self.cultivar_params['PS']
        self.TS = self.cultivar_params['TS']
        self.PVT = self.cultivar_params['PVT']
        self.FDF = self.cultivar_params['FDF']
        self.Tb = self.cultivar_params['Tb']
        self.To = self.cultivar_params['To']
        self.Tm = self.cultivar_params['Tm']

        # 土壤参数
        self.SNO3_0 = self.soil_params['SNO3_0']
        self.SNH4_0 = self.soil_params['SNH4_0']
        self.BD = self.soil_params['BD']
        self.DLAYR = self.soil_params['DLAYR']
        self.NLAYR = self.soil_params['NLAYR']
        self.DS = self.soil_params['DS']
        self.DUL = self.soil_params['DUL']
        self.SAT = self.soil_params['SAT']
        self.LL = self.soil_params['LL']
        self.PH = self.soil_params['PH']
        self.SOILC = self.soil_params['SOILC']
        self.FOM = self.soil_params['FOM']
        self.FON = self.soil_params['FON']
        self.TOTN = self.soil_params['TOTN']
        self.SW0 = self.soil_params['SW0']
        self.CLAY = self.soil_params['CLAY']
        self.CN = self.soil_params['CN']

        # 氮相关参数
        self.Nd = self.n_params['Nd']

    def interception(self, P: float, LAI: float, Fw: float, SW0: list):
        """Interception, 作物叶片截留, 单位为cm3/cm3"""

        eta = math.exp(- self.kl * LAI)
        f = 0.2  # f means 0.2kg water/kg fresh weight,f为单位鲜重生物量的截留能力（0.2kg水/kg 鲜重）
        intcap = (1-eta) * f * Fw/10000/10000*1000    # Fw means fresh weight of the crop,its unit is kg/ha. turns to g/cm2 ???? Fw = 地上部鲜重？ ,intercap(cm)
        t = 1   # 模拟步长
        Interception = min(P, intcap/t)
        if Interception <= 0:  ###补充代码，避免出现截留量<0的情况
            # print('interception < 0, please check')
            Interception = 0

        return Interception, SW0

    def runoff(self, P: float, I: float, Interception: float, SW0: list, Flood: bool=False):
        """
        计算Runoff, 地表径流量，cm

        :param Flood: 是否淹水，布尔类型，默认False
        """
        SW1 = [0 for i in range(self.NLAYR)]
        if Flood == True:
            FldH = 4   # flooding depth, cm  水田淹水深度
            Pinf = P + I - Interception
            FldH = FldH + Pinf
            BundH = 30   # height of ridge, cm  田埂高
            if FldH > BundH:
                Runoff = FldH - BundH   #### SW会改变吗？？？？？？？？！！！！！！！！
                depth = BundH
            if FldH <= BundH:
                Runoff = 0
                depth  = FldH
            # SW = []
            for L in range(self.NLAYR):
                SW1[L] = self.DUL[L]

        if Flood == False:

            # 土壤储水量 Soil storage
            # SMX = 254.0 * (100.0 / CN - 1.0) # (mm)
            SMX = 254.0 * (100.0 / self.CN - 1.0) / 10  #(cm)

            # 初始吸收比率 Initial abstraction ratio
            # 地表径流与土壤地表两层含水量的平均值有关。Runoff is related to the average soil water content of the top two layers of soil
            SWABI = 0.15 * ((self.SAT[0] - SW0[0]) / (self.SAT[0] - self.LL[0] * 0.5) + (
                        self.SAT[1] - SW0[1]) / (self.SAT[1] - self.LL[1] * 0.5))
            SWABI = max(0.0, SWABI)

            # 暂不考虑土壤表层的覆盖物
            IABS = SWABI

            # WATAVL, Water available for infiltration or runoff (rainfall plus irrigation) (cm/d)
            # 可用于入渗或径流的水量（= 降雨量+灌溉量）单位为cm/d，与降雨和灌溉的单位一致
            WATAVL = P + I

            # PB, Determines threshold amount of rainfall that will occur before runoff starts (mm/d)
            # 确定在径流开始之前，将发生的降雨或灌溉阈值（毫米/天）
            PB = WATAVL - IABS * SMX
            if WATAVL > 0.001:
                if PB > 0:
                    Runoff = PB ** 2 / (WATAVL + (1.0 - IABS) * SMX)
                else:
                    Runoff = 0.0
            else:
                Runoff = 0.0

            for L in range(self.NLAYR):
                # print(L)
                if L <= 1:
                    SW1[L] = self.DUL[L]  # 土壤前两层为饱和含水量
                else:
                    SW1[L] = self.SW0[L]

        return Runoff, SW1

    def pinf(self, P: float, I: float, Interception: float, SW1: list):
        """降水或灌溉的入渗，cm"""

        SWSAT =[]
        # SWFC = DUL  # SWFC means field water capacity,the same as DUL.
        Ksw = []
        Pinf = [0 for i in range(self.NLAYR)]
        SW2 = [0 for i in range(self.NLAYR)]

        # Interception = Soil_Interception.Interception(WEATHER, t, PDT, HI, TBiomass, SLA)
        sum_Pinf = [0 for i in range(self.NLAYR)]
        sum_Pinf[0] = P + I - Interception
        # print('sum_Pinf = ',sum_Pinf)

        for L in range(self.NLAYR):
            SWSAT.append(1 - self.BD[L] / 2.65)  #SWSAT == SAT
            # SWSAT = SAT
            Ksw.append(75 * ((SWSAT[L] - self.DUL[L]) / self.DUL[L]) ** 2)  #Its unit is cm/h, KSAT means Soil Saturated Hydraulic Conductivity. the same as Ksw.饱和导水率
            SW2[L] = (sum_Pinf[L]/self.DLAYR[L] + SW1[L])
            if SW2[L] < self.DUL[L]:
                Pinf[L] = sum_Pinf[L]
                sum_Pinf[L] = 0

                continue
            if SW2[L] == self.DUL[L]:
                if L == 0:
                    Pinf[L] = sum_Pinf[L]
                if L > 0:
                    sum_Pinf[L] = sum_Pinf[L - 1]
                continue
            if SW2[L] > self.DUL[L]:
                Pinf[L] = (self.DUL[L]-SW1[L])*self.DLAYR[L]
                SW2[L] = self.DUL[L]
                sum_Pinf[L] = sum_Pinf[L] - Pinf[L] #sum_Pinf 表示到下一层的入渗量

            if SW2[L] <= self.LL[L]:   ###补充代码，农田的含水量始终高于萎蔫含水量，避免出现含水量为0的情况
                SW2[L] == self.LL[L]
                print('SW2 < 0, please check')

        return SW2

    def capillary_rise(self, SW2: list, water_table_presence: bool=1):
        """
        计算地下水的毛细管上升量

        :param water_table_presence: 是否存在地下水（0，1）
        :return: 地下水位对毛细管的补给量
        """

        Soil_fshape_cr = 16.0  # 毛细管上升形状系数，用于控制毛细管上升水量的比例
        FluxOut = [0.04, 0.03, 0.02, 0.01]      # 从每层流出量
        z_gw = 25
        Init_water_Adj = [0, 0, 0, 0]

        zMid_1 = np.cumsum(self.DLAYR)
        zMid_2 = [item / 2 for item in self.DLAYR]
        zMid = np.subtract(zMid_1, zMid_2)

        # 计算每一层的饱和导水率 (mm/day) 饱和导水率
        Ksat = []
        for ii in range(len(self.SAT)):
            lmbda = 1 / ((np.log(1500) - np.log(33)) / (np.log(self.DUL[ii]) - np.log(self.LL[ii])))
            K1 = (1930 * (self.SAT[ii] - self.DUL[ii]) ** (3 - lmbda)) * 24
            Ksat.append(K1)
        # print(Ksat)

        # 模块输入系数总体上bCR、aCR,需要的系数
        bCR = np.zeros(self.NLAYR)
        aCR = np.zeros(self.NLAYR)
        for jj in range(0, self.NLAYR):
            if (
                    (self.LL[jj] >= 0.04)
                    and (self.LL[jj] <= 0.15)
                    and (self.DUL[jj] >= 0.09)
                    and (self.DUL[jj] <= 0.28)
                    and (self.SAT[jj] >= 0.32)
                    and (self.SAT[jj] <= 0.51)
            ):
                # Sandy soil class
                if (Ksat[jj] >= 200) and (Ksat[jj] <= 2000):
                    aCR[jj] = -0.3112 - (Ksat[jj] * (1e-5))
                    bCR[jj] = -1.4936 + (0.2416 * np.log(Ksat[jj]))
                elif Ksat[jj] < 200:
                    aCR[jj] = -0.3112 - (200 * (1e-5))
                    bCR[jj] = -1.4936 + (0.2416 * np.log(200))
                elif Ksat[jj] > 2000:
                    aCR[jj] = -0.3112 - (2000 * (1e-5))
                    bCR[jj] = -1.4936 + (0.2416 * np.log(2000))

            elif (
                    (self.LL[jj] >= 0.06)
                    and (self.LL[jj] <= 0.20)
                    and (self.DUL[jj] >= 0.23)
                    and (self.DUL[jj] <= 0.42)
                    and (self.SAT[jj] >= 0.42)
                    and (self.SAT[jj] <= 0.55)
            ):
                # Loamy soil class
                if (Ksat[jj] >= 100) and (Ksat[jj] <= 750):
                    aCR[jj] = -0.4986 + (9 * (1e-5) * Ksat[jj])
                    bCR[jj] = -2.132 + (0.4778 * np.log(Ksat[jj]))
                elif Ksat[jj] < 100:
                    aCR[jj] = -0.4986 + (9 * (1e-5) * 100)
                    bCR[jj] = -2.132 + (0.4778 * np.log(100))
                elif Ksat[jj] > 750:
                    aCR[jj] = -0.4986 + (9 * (1e-5) * 750)
                    bCR[jj] = -2.132 + (0.4778 * np.log(750))

            elif (
                    (self.LL[jj] >= 0.16)
                    and (self.LL[jj] <= 0.34)
                    and (self.DUL[jj] >= 0.25)
                    and (self.DUL[jj] <= 0.45)
                    and (self.SAT[jj] >= 0.40)
                    and (self.SAT[jj] <= 0.53)
            ):
                # Sandy clayey soil class
                if (Ksat[jj] >= 5) and (Ksat[jj] <= 150):
                    aCR[jj] = -0.5677 - (4 * (1e-5) * Ksat[jj])
                    bCR[jj] = -3.7189 + (0.5922 * np.log(Ksat[jj]))
                elif Ksat[jj] < 5:
                    aCR[jj] = -0.5677 - (4 * (1e-5) * 5)
                    bCR[jj] = -3.7189 + (0.5922 * np.log(5))
                elif Ksat[jj] > 150:
                    aCR[jj] = -0.5677 - (4 * (1e-5) * 150)
                    bCR[jj] = -3.7189 + (0.5922 * np.log(150))

            elif (
                    (self.LL[jj] >= 0.20)
                    and (self.LL[jj] <= 0.42)
                    and (self.DUL[jj] >= 0.40)
                    and (self.DUL[jj] <= 0.58)
                    and (self.SAT[jj] >= 0.49)
                    and (self.SAT[jj] <= 0.58)
            ):

                # Silty clayey soil class
                if (Ksat[jj] >= 1) and (Ksat[jj] <= 150):
                    # print(aCR)
                    aCR[jj] = -0.6366 + (8 * (1e-4) * Ksat[jj])
                    bCR[jj] = -1.9165 + (0.7063 * np.log(Ksat[jj]))
                elif Ksat[jj] < 1:
                    aCR[jj] = -0.6366 + (8 * (1e-4) * 1)
                    bCR[jj] = -1.9165 + (0.7063 * np.log(1))
                elif Ksat[jj] > 150:
                    aCR[jj] = -0.6366 + (8 * (1e-4) * 150)
                    bCR[jj] = -1.9165 + (0.7063 * np.log(150))

            for i in range(len(aCR)):
                aCR[i] = max(aCR[i], 0)
                bCR[i] = max(bCR[i], 0)

        # print(aCR, bCR)

        ## Get groundwater table elevation on current day ##
        z_gw = z_gw

        ## Calculate capillary rise ##
        if water_table_presence == 0:  # No water table present
            # Capillary rise is zero
            CrTot = 0
        elif water_table_presence == 1:  # Water table present
            # Get maximum capillary rise for bottom compartment
            zBot = zMid_1[-1]
            zBotMid = zMid[-1]
            if (Ksat[-1] > 0) and (z_gw > 0) and ((z_gw - zBotMid) < 4):
                if zBotMid >= z_gw:
                    MaxCR = 99
                else:
                    MaxCR = np.exp((np.log(z_gw - zBotMid) - bCR[-1]) / aCR[-1])
                    if MaxCR > 99:
                        MaxCR = 99

            else:
                MaxCR = 0


            compi = self.NLAYR - 1  # 从土层最底部开始
            WCr = 0  # 毛细管上升计数器

            Init_water = SW2

            while (round(MaxCR * 1000) > 0) and (compi > -1) and (round(FluxOut[compi] * 1000) == 0):
                #  向上运动，达到土壤表面，或遇到当天已经发生向下排水渗透的隔室，找到当前土层，计算驱动力
                if (Init_water[compi] >= self.LL[compi]) and (Soil_fshape_cr > 0):
                    Df = 1 - (
                        (
                            (Init_water[compi] - self.LL[compi])
                            / (Init_water_Adj[compi] - self.LL[compi])
                        )
                        ** Soil_fshape_cr
                    )
                    if Df > 1:
                        Df = 1
                    elif Df < 0:
                        Df = 0

                else:
                    Df = 1

                #  计算相对水力传导系数
                thThr = (self.LL[compi] + self.DUL[compi]) / 2
                if Init_water[compi] < thThr:
                    if (Init_water[compi] <= self.LL[compi]) or (thThr <= self.LL[compi]):
                        Krel = 0
                    else:
                        Krel = (Init_water[compi] - self.LL[compi]) / (thThr - self.LL[compi])

                else:
                    Krel = 1

                #   检查是否有空间储存毛细管上升的水
                dth = round(Init_water_Adj[compi] - Init_water[compi],4)

                # Store water if room is available  如果空间有效的话，存储水量
                if (dth > 0) and ((zBot - self.DLAYR[compi] / 2) < z_gw):
                    dthMax = Krel * Df * MaxCR / (10 * self.DLAYR[compi])
                    if dth >= dthMax:
                        Init_water[compi] = Init_water[compi] + dthMax
                        CRcomp = dthMax * 10 * self.DLAYR[compi]
                        MaxCR = 0
                    else:
                        Init_water[compi] = Init_water_Adj[compi]
                        CRcomp = dth * 10 * self.DLAYR[compi]
                        MaxCR = (Krel * MaxCR) - CRcomp

                    WCr = WCr + CRcomp

                #  更新各层底部高度
                zBot = zBot - self.DLAYR[compi]
                #  更新隔间和层计数器
                compi = compi - 1
                #   更新毛细管最大上升限制
                if compi > -1:

                    zBotMid = zBot - (self.DLAYR[compi] / 2)
                    if (Ksat[compi] > 0) and (z_gw > 0) and ((z_gw - zBotMid) < 4):
                        if zBotMid >= z_gw:
                            LimCR = 99
                        else:
                            LimCR = np.exp((np.log(z_gw - zBotMid) - bCR[compi]) / aCR[compi])
                            if LimCR > 99:
                                LimCR = 99

                    else:
                        LimCR = 0

                    if MaxCR > LimCR:
                        MaxCR = LimCR

            # Store total depth of capillary rise  毛细管上升总深度
            CrTot = WCr

        return CrTot, SW2

    def e0(self, T: float):
        """
        :param T: 气温，℃
        :return:
            e: 气温为T时的水汽压，kPa
        """
        e = 0.6108*(math.exp((17.27 * T / (T + 237.3))))
        return e

    def cal_Rn(self, doy: int, SH: float, Tmax: float, Tmin: float, RH: float) -> float:
        """作物表面上的净辐射，MJ/(m2·day)"""

        a = 0.23       # 冠层反射率，绿色植被覆盖的地表的反射率约为0.20-0.25。以绿草为参考作物时，a=0.23。
        Lat = self.latitude * math.pi / 180.0
        Z = self.elevation
        delta = 0.409 * math.sin(2 * math.pi * doy / 365 - 1.39)        # 太阳赤纬
        ws = math.acos(-math.tan(Lat*math.pi / 180.0)*math.tan(delta))        # 太阳时角
        dr = 1 + 0.033 * math.cos(2 * math.pi * doy / 365)      # 日地距离系数
        Ra = ((24 * 60 * 0.082) / math.pi * dr *
              (ws * math.sin(Lat*math.pi / 180.0) * math.sin(delta) + math.sin(ws) * math.cos(Lat*math.pi / 180.0) * math.cos(delta)))      # 天文辐射
        Rso = (0.75 + 2 * 10**(-5) * Z) * Ra        # 晴天辐射
        N = 24 / math.pi * ws
        Rs = (0.25 + 0.5 * SH / N) * Ra         # 太阳短波辐射
        Rns = (1 - a) * Rs  # 净短波辐射

        es = (self.e0(Tmax) + self.e0(Tmin)) / 2
        ea = es * RH / 100
        Rnl = (4.903 * 10**(-9) * (((Tmax + 273.16)**4 - (Tmin + 273.16)**4) / 4) *
               (0.34 - 0.14 * ea**0.5) * (1.35 * Rs / Rso - 0.35))      # 净长波辐射
        Rn = Rns - Rnl
        return Rn

    def PM_ET0(self, Tmax, Tmin, RH, P, u2, Rn):
        """作物潜在蒸散量，mm"""

        Tmean = (Tmax + Tmin) / 2
        # 土壤热通量
        G = 0
        # 饱和水汽压
        es = (self.e0(Tmax) + self.e0(Tmin)) / 2
        # 实际水汽压
        ea = es * RH / 100
        # 饱和水汽压曲线斜率
        delta = 4098 * self.e0(Tmean) / (Tmean + 237.3) ** 2
        # 湿度计常数
        gama = 0.665 * 0.001 * P
        ET0 = (0.408*delta*(Rn-G) + gama*900*u2*(es-ea)/(Tmean+273)) / (delta + gama*(1 + 0.34*u2))

        return ET0

    def cal_Kcb(self, current_pdt: float, RHmin: float, u2: float) -> float:
        """Kcb，作物蒸腾系数"""

        PDT_1 = self.PDT_EM
        PDT_2 = self.PDT_BT
        PDT_3 = self.PDT_EF
        PDT_4 = self.PDT_MT

        if RHmin != 45 or u2 != 2:
            if self.Kcbmid >= 0.45:
                kcbmid_new = self.Kcbmid + (0.04*(u2-2)-0.004*(RHmin-45))*(self.h/3)**0.3
            else:
                kcbmid_new = self.Kcbmid

            if self.Kcbend >= 0.45:
                kcbend_new = self.Kcbend + (0.04*(u2-2)-0.004*(RHmin-45))*(self.h/3)**0.3
            else:
                kcbend_new = self.Kcbend
        else:
            kcbmid_new = self.Kcbmid
            kcbend_new = self.Kcbend

        if current_pdt <= PDT_1:
            Kcb = self.Kcbini
        elif PDT_1 < current_pdt <= PDT_2:
            Kcb = self.Kcbini + (kcbmid_new - self.Kcbini)/(PDT_2 - PDT_1) * (current_pdt - PDT_1)
        elif PDT_2 < current_pdt <= PDT_3:
            Kcb = kcbmid_new
        elif PDT_3 < current_pdt < PDT_4:
            Kcb = self.Kcbmid + (kcbend_new - kcbmid_new)/(PDT_4 - PDT_3) * (current_pdt - PDT_3)
        else:
            Kcb = kcbend_new
        return Kcb

    def ETc(self, current_date: str, current_pdt: float, rootdepth: float, Runoff: float, I: float, CrTot: float, De0: float, SW0: list):
        """ETcadj: 胁迫条件下的作物蒸腾量，mm"""

        Rain = self.weather.loc[current_date, 'RAIN']
        Tmin = self.weather.loc[current_date, 'Tmin']
        Tmax = self.weather.loc[current_date, 'Tmax']
        SH = self.weather.loc[current_date, 'SH']
        u2 = self.weather.loc[current_date, 'u2']
        RH = self.weather.loc[current_date, 'RHave']
        RHmin = self.weather.loc[current_date, 'RHmin']
        P = self.weather.loc[current_date, 'Atm']

        # 计算ETo
        doy = DateUtils.cal_Doy(current_date)
        Rn = self.cal_Rn(doy, SH, Tmax, Tmin, RH)
        ETo = self.PM_ET0(Tmax, Tmin, RH, P, u2, Rn)         # 参考作物蒸散量，mm/d

        kcb = self.cal_Kcb(current_pdt, RHmin, u2)      # 作物蒸腾系数
        # print(kcb)
        kcmin = 0.15        # 无地表覆盖的干燥土壤最小Kc值，约为 0.15－0.20
        kcmax = max((1.2 + (0.04 * (u2 - 2) - 0.004 * (RHmin - 45)) * (self.h / 3) ** 0.3), (kcb + 0.5))
        if kcb >= kcmin:
            fc = ((kcb - kcmin) / (kcmax - kcmin)) ** (1 + 0.5 * self.h)  # 地表遮盖比例
        else:
            fc = 0
        fw = 1.0  # 地表湿润比例，降雨、喷灌、漫灌、畦灌条件下为1.0
        few = min((1 - fc), fw)

        # 计算蒸发减少系数kr
        Ze = 0.10  # 由于蒸发而变干的表土层厚度，可取0.10~0.15m
        TEW = 1000 * (self.DUL[0] - 0.5 * self.LL[0]) * Ze  # 总蒸发水量
        REW = 8  # 易蒸发水量
        if De0 <= REW:
            kr = 1
        else:
            kr = (TEW - De0) / (TEW - REW)

        # 计算地表蒸发系数ke和地表蒸发量E
        ke = min(kr*(kcmax-kcb), few*kcb)
        E = kr * ke * ETo       # 地表蒸发量，mm

        Tew = 0        # 通常可忽略

        # 当有强降水或灌溉时，De可认为是0
        if Rain >= 3 or I > 0:
            De = 0
            DPe = (Rain - Runoff) + I / fw - De0
            SW = SW0
            SW[0] = self.DUL[0]
        else:
            DPe = 0
            De = De0 - (Rain - Runoff) - I/fw + E/few + Tew + DPe
            if De < 0:
                De = 0
            elif De > TEW:
                De = TEW
            SW = SW0
            SW[0] = SW0[0] - De / (Ze * 1000)

        ETcadj = 0
        if rootdepth > 0:

            # ======================== 计算 ET =========================
            Zr = rootdepth / 10      # 每日有效根深，m
            # Zr = 0.5

            Depth = 0
            if Zr > 0.4:
                i = 3
                Zr = 0.4
            else:
                i = 0
                while i < self.NLAYR:
                    Depth += self.DLAYR[i] / 100
                    if Zr > Depth:
                        i += 1
                    else:
                        break

            # print(i)

            TAW = 1000 * (self.DUL[i] - self.LL[i]) * Zr  # 总有效水量，mm
            pct = 0.5  # 不遭受水分胁迫时，作物从根系层中吸收的有效水量与总有效水量 TAW 之比，通常取0.5

            Dr0 = 1000 * (self.DUL[i] - SW0[i]) * Zr
            if Dr0 > TAW:
                Dr0 = TAW
            if Dr0 <= pct*TAW:
                ks = 1.0
            else:
                ks = (TAW - Dr0) / ((1 - pct) * TAW)
            ETc = (kcb + ke) * ETo
            # ETcadj = (ks*kcb + ke) * ETo
            ETcadj = (kcb + ke) * ETo
            # 当有强降水或灌溉时，De可认为是0
            if Rain >= 3 or I > 0:
                DP = (Rain - Runoff) + I - ETcadj - Dr0
            else:
                DP = 0

            Dr = Dr0 - (Rain - Runoff) - I - CrTot + ETcadj + DP
            if Dr > TAW:
                Dr = TAW
            SW[i] = SW0[i] - Dr / (Zr * 1000)

        return E, ETcadj, De, SW

    def simulate(self):
        """根系吸水量，单位为 cm3/cm3"""

        biomassmodel = BiomassModel(self.start, self.end, self.weather, self.cultivar_params, self.management,
                                    self.plant_params, self.soil_params)
        biomassmodel.ifEmergence()

        results = {}

        delta = (self.end_date - self.start_date).days
        SW = self.SW0
        De0 = 0
        Biomass0 = self.seed_weight / 2
        rootdepth = 0
        for day in (self.start_date + datetime.timedelta(n) for n in range(delta + 1)):
            date_str = day.strftime("%m/%d/%Y")

            Biomass0, LAI, Fw, ROOTWT = BiomassModel(self.start, date_str, self.weather, self.cultivar_params, self.management,
                                           self.plant_params, self.soil_params).biomass(day, Biomass0)

            irri_date = self.irrigation['date']
            if date_str == irri_date:
                I = self.irrigation['amount']
            else:
                I = 0.0
            try:
                Rain = self.weather.loc[date_str, 'RAIN']
                Tmin = self.weather.loc[date_str, 'Tmin']
                Tmax = self.weather.loc[date_str, 'Tmax']
                SH = self.weather.loc[date_str, 'SH']
                u2 = self.weather.loc[date_str, 'u2']
                RH = self.weather.loc[date_str, 'RHave']
                RHmin = self.weather.loc[date_str, 'RHmin']
                # P = self.weather.loc[date_str, 'Atm']
            except KeyError:
                print(f"气象数据中缺少 {date_str} 的记录，跳过该天。")
                continue
            if day < biomassmodel.emergence_date:
                print("作物还未出苗。")
                PDT = 0
                Interception = 0
                Runoff, SW1 = self.runoff(Rain, I, Interception, SW)
                SW2 = self.pinf(Rain, I, Interception, SW1)
                CrTot, SW2 = self.capillary_rise(SW2)
                E, ETcadj, De, SW3 = self.ETc(date_str, PDT, rootdepth, Runoff, I, CrTot, De0, SW2)
                Si = [0 for i in range(self.NLAYR)]
                SW4 = [0 for i in range(self.NLAYR)]
                for L in range(self.NLAYR):
                    SW4[L] = min(max(self.LL[L], SW3[L] - Si[L]), self.SAT[L])  ### 土壤水分含量最低为LL
            else:
                Interception, SW0 = self.interception(Rain, LAI, Fw, SW)
                PDT = biomassmodel._cal_pdt(day)

                daily_rootdepth = biomassmodel.daily_DRTDEP(date_str, PDT, wdf=1)
                rootdepth += daily_rootdepth
                # 达到最大有效根深后根系深度不再增加
                if rootdepth == 150:
                    rootdepth = 150
                RLV = biomassmodel.RLV(rootdepth, ROOTWT)

                Runoff, SW1 = self.runoff(Rain, I, Interception, SW0)
                SW2 = self.pinf(Rain, I, Interception, SW1)
                CrTot, SW2 = self.capillary_rise(SW2)
                E, ETcadj, De, SW3 = self.ETc(date_str, PDT, rootdepth, Runoff, I, CrTot, De0, SW2)

                Tp = ETcadj - E  # Tp means potential transpiration of plants

                Si = []
                SW4 = [0 for i in range(self.NLAYR)]
                sum_RLV = sum(RLV)
                for L in range(self.NLAYR):
                    ccc = RLV[L] / sum_RLV
                    Smi = Tp * ccc
                    WFi = 0.4  # WFi means the effct on Soil water transpiration, and its range is 0-1.
                    # WFi = 1
                    # print(min(SW3[L] - self.LL[L], WFi * Smi))
                    Si.append(max(min(SW3[L] - self.LL[L], WFi * Smi), 0))  ### 作物仅能吸收LL以上的土壤水分
                    # sum_Ta = sum_Ta + Si
                    SW4[L] = min(max(self.LL[L], SW3[L] - Si[L]), self.SAT[L])  ### 土壤水分含量最低为LL

            results[date_str] = Si
            SW = SW4
            De0 = De
        return results


# Pydantic models for API
class SimulationRequest(BaseModel):
    start_date: str
    end_date: str
    cultivar_params: Dict[str, Any]
    plant_params: Dict[str, Any]
    soil_params: Dict[str, Any]
    management: Dict[str, Any]
    n_params: Dict[str, Any]
    weather_data: List[Dict[str, Any]]  # Weather data as list of dicts

class HealthCheckResponse(BaseModel):
    status: str
    message: str

# FastAPI app
app = FastAPI(
    title="Soil Water Balance Model API",
    description="API for calculating soil water balance and root water uptake",
    version="1.0.0"
)

@app.get("/health", response_model=HealthCheckResponse)
async def health_check():
    """Health check endpoint"""
    return HealthCheckResponse(status="healthy", message="Soil water balance model API is running")

@app.post("/simulate", response_model=Dict[str, Any])
async def simulate_soil_water_balance(request: SimulationRequest):
    """
    Run soil water balance simulation

    Parameters:
    - start_date: Start date in MM/DD/YYYY format
    - end_date: End date in MM/DD/YYYY format
    - cultivar_params: Cultivar parameters
    - plant_params: Plant parameters
    - soil_params: Soil parameters
    - management: Management parameters
    - n_params: Nitrogen parameters
    - weather_data: Weather data as list of dictionaries

    Returns:
    - Dictionary containing simulation results with dates as keys and root water uptake (Si) as values
    """
    try:
        # Convert weather data list to pandas DataFrame
        weather_df = pd.DataFrame(request.weather_data)
        if 'date' in weather_df.columns:
            weather_df.set_index('date', inplace=True)
            weather_df.index = pd.to_datetime(weather_df.index)

        # Initialize the model
        model = SoilWaterModel(
            start=request.start_date,
            end=request.end_date,
            weather=weather_df,
            cultivar_params=request.cultivar_params,
            plant_params=request.plant_params,
            soil_params=request.soil_params,
            management=request.management,
            n_params=request.n_params
        )

        # Run simulation
        results = model.simulate()

        return {
            "status": "success",
            "results": results,
            "start_date": request.start_date,
            "end_date": request.end_date
        }

    except Exception as e:
        raise HTTPException(status_code=500, detail=f"Simulation failed: {str(e)}")

@app.post("/simulate_from_files")
async def simulate_from_files(
    start_date: str,
    end_date: str,
    weather_file: UploadFile = File(...),
    management_file: UploadFile = File(...),
    soil_params_file: UploadFile = File(...),
    n_params_file: UploadFile = File(...),
    plant_params_file: UploadFile = File(...),
    cultivar_params_file: UploadFile = File(...)
):
    """
    Run simulation using uploaded parameter files

    Parameters:
    - start_date: Start date in MM/DD/YYYY format
    - end_date: End date in MM/DD/YYYY format
    - weather_file: CSV file with weather data
    - management_file: JSON file with management parameters
    - soil_params_file: JSON file with soil parameters
    - n_params_file: JSON file with nitrogen parameters
    - plant_params_file: JSON file with plant parameters
    - cultivar_params_file: JSON file with cultivar parameters

    Returns:
    - Dictionary containing simulation results
    """
    try:
        # Read uploaded files
        weather_data = pd.read_csv(weather_file.file, index_col='date', parse_dates=True)

        management = json.load(management_file.file)
        soil_params = json.load(soil_params_file.file)
        n_params = json.load(n_params_file.file)
        plant_params = json.load(plant_params_file.file)
        cultivar_params = json.load(cultivar_params_file.file)

        # Initialize the model
        model = SoilWaterModel(
            start=start_date,
            end=end_date,
            weather=weather_data,
            cultivar_params=cultivar_params,
            plant_params=plant_params,
            soil_params=soil_params,
            management=management,
            n_params=n_params
        )

        # Run simulation
        results = model.simulate()

        return {
            "status": "success",
            "results": results,
            "start_date": start_date,
            "end_date": end_date
        }

    except Exception as e:
        raise HTTPException(status_code=500, detail=f"Simulation failed: {str(e)}")

if __name__ == '__main__':
    # Run with uvicorn
    uvicorn.run(app, host="0.0.0.0", port=8000)