package com.dy.pipIrrModel.modelCalculate;

import net.objecthunter.exp4j.Expression;
import net.objecthunter.exp4j.ExpressionBuilder;

/**
 * @Author: liurunyu
 * @Date: 2025/8/18 14:51
 * @Description 哈格里夫斯（Hargreaves）公式 计算作物实际蒸散量（mm/day）
 */
public class Hargreaves {

    /**
     * 计算弧度
     * 弧度=角度× π/180
     * @param lat 地理纬度
     * @return
     */
    public static Double rad(Double lat) {
        Expression expression = new ExpressionBuilder("x * pi / 180")
                .variables("x")
                .build();
        // 设置变量值
        expression.setVariable("x", lat);
        Double result = expression.evaluate();
        return result ;
    }

    /**
     * 太阳磁偏角
     * @param dayIndex 为年内某天的日序数（比如1月1日为1，12月31日为365或366）
     * @return
     */
    public static Double sunMagnetismAngular(Integer dayIndex){
        Expression expression = new ExpressionBuilder("0.409 * sin(((2 * pi) / 365) * x - 1.39)")
                .variables("x")
                .build();
        // 设置变量值
        expression.setVariable("x", dayIndex);
        Double result = expression.evaluate();
        return result ;
    }

    /**
     * 日地间相对距离的倒数
     * @param dayIndex 为年内某天的日序数（比如1月1日为1，12月31日为365或366）
     * @return
     */
    public static Double sunEarthDistance(Integer dayIndex){
        Expression expression = new ExpressionBuilder("1 + (0.033 * cos(((2 * pi) / 365) * x))")
                .variables("x")
                .build();
        // 设置变量值
        expression.setVariable("x", dayIndex);
        Double result = expression.evaluate();
        return result ;
    }


    /**
     * 太阳时角
     * 0.409 * sin(2 * pi / 365 * x - 1.39)
     * @param fai 地理弧度
     * @param sunMagnetismAngular 太阳磁偏角
     * @return
     */
    public static Double sunTimeAngular(Double fai, Double sunMagnetismAngular){
        Expression expression = new ExpressionBuilder("acos(-tan(x) * tan(y))")
                .variables("x", "y")
                .build();
        // 设置变量值
        expression.setVariable("x", fai);
        expression.setVariable("y", sunMagnetismAngular);
        Double result = expression.evaluate();
        return result ;
    }

    /**
     * 计算天顶辐射
     * @param sunEarthDistance 日地间相对距离的倒数
     * @param sunTimeAngular 太阳时偏角
     * @param fai 地理弧度
     * @param sunMagnetismAngular 太阳磁偏角
     * @return
     */
    public static Double zenithRadiation(Double sunEarthDistance,
                                         Double sunTimeAngular,
                                         Double fai,
                                         Double sunMagnetismAngular){
        Expression expression = new ExpressionBuilder("((24 * 60) / pi) * 0.082 * a * ((b * sin(c) * sin(d)) + (cos(c) * cos(d) * sin(b)))")
                .variables("a", "b", "c", "d")
                .build();
        // 设置变量值
        expression.setVariable("a", sunEarthDistance);
        expression.setVariable("b", sunTimeAngular);
        expression.setVariable("c", fai);
        expression.setVariable("d", sunMagnetismAngular);
        Double result = expression.evaluate();
        return result ;
    }

    /**
     * 格里夫斯（Hargreaves）公式计算作物蒸散量（mm/day）
     * @param kc 作物系数
     * @param maxT 一日内最高温度
     * @param maxT 一日内最高温度
     * @param zenithRadiation 天顶辐射
     * @return
     */
    public static Double ET0(Double kc, Double maxT, Double minT, Double zenithRadiation){
        Expression expression = new ExpressionBuilder("x * (0.0023 * (((a + b) / 2) + 17.8) * ((a - b)^(1/2)) * c * 0.408)")
                .variables("x", "a", "b", "c")
                .build();
        // 设置变量值
        expression.setVariable("x", kc);
        expression.setVariable("a", maxT);
        expression.setVariable("b", minT);
        expression.setVariable("c", zenithRadiation);
        Double result = expression.evaluate();
        return result ;
    }


    public static void main(String[] args) {
        Double lat = 38.56064095143279 ;//纬度
        Integer dayIndex = 239 ;//年内日序数（比如1月1日为1，12月31日为365或366）
        Double kc = 0.41 ;//作物系数
        Double maxT = 40.1 ;//一日内最高温度
        Double minT = 28.1 ;//一日内最低温度

        Double fai = rad(lat);
        System.out.println("弧度:" + fai);

        Double sunMagnetismAngular = sunMagnetismAngular(dayIndex);
        System.out.println("太阳磁偏角:" + sunMagnetismAngular);

        Double sunTimeAngular = sunTimeAngular(fai, sunMagnetismAngular);
        System.out.println("太阳时角:" + sunTimeAngular);

        Double sunEarthDistance = sunEarthDistance(dayIndex);
        System.out.println("日地间相对距离的倒数:" + sunEarthDistance);

        Double zenithRadiation = zenithRadiation(sunEarthDistance, sunTimeAngular, fai, sunMagnetismAngular);
        System.out.println("天顶辐射:" + zenithRadiation);

        Double et0 = ET0(kc, maxT, minT, zenithRadiation);
        System.out.println("蒸散量:" + et0);
    }

}