diff --git a/src/main/java/com/gunshi/project/xyt/service/ForecastResultsService.java b/src/main/java/com/gunshi/project/xyt/service/ForecastResultsService.java index fade810..a1ac830 100644 --- a/src/main/java/com/gunshi/project/xyt/service/ForecastResultsService.java +++ b/src/main/java/com/gunshi/project/xyt/service/ForecastResultsService.java @@ -26,6 +26,7 @@ import com.gunshi.project.xyt.model.StRsvrR; import com.gunshi.project.xyt.model.StStbprpB; import com.gunshi.project.xyt.model.StZqrlB; import com.gunshi.project.xyt.model.StZvarlB; +import com.gunshi.project.xyt.util.BigdecimalUtil; import com.gunshi.project.xyt.util.DataHandleUtil; import lombok.extern.slf4j.Slf4j; import org.springframework.beans.factory.annotation.Autowired; @@ -760,4 +761,106 @@ public class ForecastResultsService extends ServiceImpl + * @auther: cxw + * @date: 2024-11-06, 周三, 10:16:26 + */ + private List calcU(BigDecimal areaF, BigDecimal lengthL, BigDecimal j, double dt, BigDecimal tl){ + // 计算m1、n + Map m1AndNMap = calcM1AndN(areaF, lengthL, j); + BigDecimal m1 = m1AndNMap.get("m1"); + BigDecimal n = m1AndNMap.get("n"); + // k = m1 / n。如果m1有非线性改正,则为改正后的值。算例保留1位小数 + BigDecimal k = m1.divide(n, 1, BigDecimal.ROUND_HALF_UP); + /* + * 瞬时单位线u(0,t)=(1/(k*Γ(n)))*((t/k)^(n-1))*e^(-t/k) + */ + // 计算伽马函数 + BigDecimal nGamma = BigDecimal.ONE.divide(k.multiply(BigdecimalUtil.bigDecimalGamma(n)), 4, BigDecimal.ROUND_HALF_UP); + List uList = new ArrayList<>(); + for (int i = 1; i <= 11; i++) { + BigDecimal t = BigDecimal.valueOf(dt * i); + if (t.compareTo(tl) > 0) { + break; + } + // 计算t、k、n + BigDecimal tkn = BigdecimalUtil.bigDecimalExponentiation(t.divide(k, 4, BigDecimal.ROUND_HALF_UP), n.subtract(BigDecimal.ONE).doubleValue()).setScale(2, BigDecimal.ROUND_HALF_UP); + // 计算 e 的指数 + BigDecimal eToThePower = BigdecimalUtil.bigDecimalExponentiation(BigDecimal.valueOf(Math.E), t.divide(k, 4, BigDecimal.ROUND_HALF_UP).multiply(BigDecimal.ONE.negate()).setScale(2, BigDecimal.ROUND_HALF_UP).doubleValue()); + // 计算u(0,t) + BigDecimal u0t = nGamma.multiply(tkn).multiply(eToThePower).setScale(2, BigDecimal.ROUND_HALF_UP); + uList.add(u0t); + } + System.out.println(uList); + // TODO S(t)=0∫t*u(0,t)dt为参数n、k的函数,有S(t)查算表? + + // TODO 时段单位线公式u(∆t,t)=[S(t)-S(t-∆t)] + + // TODO h毫米净雨时段单位线q(∆t,t)=Fh/(3.6*∆t)*[u(∆t,t)] + + + return new ArrayList<>(); + } + + + /** + * @description: 计算m1和n参数 + * @param areaF + * @param lengthL + * @param j + * @return: void + * @auther: cxw + * @date: 2024-11-06, 周三, 13:16:24 + */ + private Map calcM1AndN(BigDecimal areaF, BigDecimal lengthL, BigDecimal j) { + BigDecimal m1 = BigDecimal.ZERO; + BigDecimal n = BigDecimal.ZERO; + // 所属片区(全省山丘区瞬时单位线分三个片进行参数的地区综合,第一片包括水文分区 1、2、4区即京广线两侧及鄂东黄冈、咸宁地区一带,(江汉平原湖区在外);第二片包括水文分区6、8、9、11区即鄂北,鄂西北及宜昌地区长江以北一带;第三片包括7、10水文分区即清江流域,恩施地区) + String slice = "I"; + if (slice.toUpperCase().equals("I")) {// I片(1、2、4区) + if (areaF.compareTo(new BigDecimal("30")) > 0) {// areaF > 30 + // m1 = 0.82*F^0.29*L^0.23*j^-0.20 + m1 = new BigDecimal("0.82").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.29)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.23)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.20)); + } else if (areaF.compareTo(new BigDecimal("30")) <= 0) {// areaF <= 30 + // m1 = 1.38*F^0.27*L^0.216*j^-0.185 + m1 = new BigDecimal("1.38").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.27)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.216)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.185)); + } + if (j.compareTo(new BigDecimal("5")) > 0) {// j > 5 + // n = 0.34*F^0.35*j^0.1 + n = new BigDecimal("0.34").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.35)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.1)); + } else if (j.compareTo(new BigDecimal("5")) <= 0) {// j <= 5 + // n = 1.04*F^0.3*j^0.1 + n = new BigDecimal("1.04").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.3)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.1)); + } + } else if (slice.toUpperCase().equals("II")) {// II片(6、8、9、11区) + // m1 = 1.64*F^0.231*L^0.131*j^-0.08 + m1 = new BigDecimal("1.64").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.231)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.131)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.08)); + // n = 0.529*F^0.25*j^0.20 + n = new BigDecimal("0.529").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.25)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.20)); + } else if (slice.toUpperCase().equals("III")) {// III片(7、10区) + // m1 = 0.8*F^0.3*L^0.1*j^-0.06 + m1 = new BigDecimal("0.8").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.3)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.1)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.06)); + // n = 0.69*F^0.224*j^0.092 + n = new BigDecimal("0.69").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.224)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.092)); + } + // TODO 特殊流域适用的公式 + + // TODO 当设计雨强超过10mm/h时,除熔岩地区外一般应作参数m1外延的非线性改正 + + BigDecimal finalM1 = m1.setScale(1, BigDecimal.ROUND_HALF_UP); + BigDecimal finalN = n.setScale(1, BigDecimal.ROUND_HALF_UP); + return new HashMap<>() {{ + put("m1", finalM1); + put("n", finalN); + }}; + } } diff --git a/src/main/java/com/gunshi/project/xyt/util/BigdecimalUtil.java b/src/main/java/com/gunshi/project/xyt/util/BigdecimalUtil.java new file mode 100644 index 0000000..cf1284f --- /dev/null +++ b/src/main/java/com/gunshi/project/xyt/util/BigdecimalUtil.java @@ -0,0 +1,51 @@ +package com.gunshi.project.xyt.util; + +import java.math.BigDecimal; +import java.math.MathContext; +import java.math.RoundingMode; + +/** + * @author cxw + * @description + * @classname BigdecimalUtil.java + * @create 2024-11-06, 星期三, 10:33:41 + */ +public class BigdecimalUtil { + + + /** + * @description: 小数指数计算 + * @param base 底数 + * @param exponent 指数 + * @return: java.math.BigDecimal + * @auther: cxw + * @date: 2024-11-06, 周三, 10:34:47 + */ + public static BigDecimal bigDecimalExponentiation(BigDecimal base, double exponent) { + BigDecimal result = BigDecimal.valueOf(Math.pow(base.doubleValue(), exponent)).round(new MathContext(10, RoundingMode.HALF_UP)); + return result; + } + + + /** + * @description: 计算x的伽马函数 + * @param x + * @return: java.math.BigDecimal + * @auther: cxw + * @date: 2024-11-06, 周三, 11:31:24 + */ + public static BigDecimal bigDecimalGamma(BigDecimal x) { + BigDecimal result = BigDecimal.valueOf(gamma(x.doubleValue())); + return result.setScale(x.scale(), RoundingMode.HALF_UP); + } + public static double gamma(double x) { + return Math.exp(logGamma(x)); + } + public static double logGamma(double x) { + double tmp = (x - 0.5) * Math.log(x + 4.5) - (x + 4.5); + double ser = 1.0 + 76.18009173 / (x + 0) - 86.50532033 / (x + 1) + + 24.01409822 / (x + 2) - 1.231739516 / (x + 3) + + 0.00120858003 / (x + 4) - 0.00000536382 / (x + 5); + return tmp + Math.log(ser * Math.sqrt(2 * Math.PI)); + } +}