minor additions

This commit is contained in:
Sergey Chernov 2026-04-01 02:50:14 +03:00
parent c140567e0c
commit 446c8d9a6e
2 changed files with 68 additions and 71 deletions

View File

@ -1,28 +1,13 @@
/*
Рассчитывает глубину провала по времени падения камня и прихода звука.
@param T измеренное полное время (с)
@param m масса камня (кг)
@param d диаметр камня (м) (предполагается сферическая форма)
@param rho плотность воздуха (кг/м³), по умолчанию 1.2
@param c скорость звука (м/с), по умолчанию 340.0
@param g ускорение свободного падения (м/с²), по умолчанию 9.81
@param Cd коэффициент лобового сопротивления, по умолчанию 0.5
@param epsilon относительная точность (м), по умолчанию 1e-3
@param maxIter максимальное число итераций, по умолчанию 20
@return глубина h (м), или null если расчёт не сошёлся
*/
fun calculateDepth( fun calculateDepth(
T: Real, T: Real,
m: Real, m: Real,
d: Real, d: Real,
rho: Real = 1.2, rho: Real = 1.2,
c: Real = 340.0, c: Real = 340.0,
g: Real = 9.81, g: Real = 9.81,
Cd: Real = 0.5, Cd: Real = 0.5,
epsilon: Real = 1e-3, eps: Real = 1e-3,
maxIter: Int = 20 maxIter: Int = 100
): Real? { ): Real? {
// Площадь миделя // Площадь миделя
val r = d / 2.0 val r = d / 2.0
@ -36,69 +21,56 @@ fun calculateDepth(
// Функция времени падения с высоты h // Функция времени падения с высоты h
fun tFall(h: Real): Real { fun tFall(h: Real): Real {
// Для численной стабильности при больших h используем логарифмическую форму
val arg = exp(g * h / (vTerm * vTerm)) val arg = exp(g * h / (vTerm * vTerm))
// arcosh(x) = ln(x + sqrt(x^2 - 1)) // arcosh(x) = ln(x + sqrt(x^2 - 1))
val arcosh = ln(arg + sqrt(arg * arg - 1.0)) val arcosh = ln(arg + sqrt(arg * arg - 1.0))
return vTerm / g * arcosh return vTerm / g * arcosh
} }
// Производная времени падения по h // Полное расчётное время
fun dtFall_dh(h: Real): Real {
val expArg = exp(2.0 * g * h / (vTerm * vTerm))
return 1.0 / (vTerm * sqrt(expArg - 1.0))
}
// Полное расчётное время T_calc(h) = tFall(h) + h/c
fun Tcalc(h: Real): Real = tFall(h) + h / c fun Tcalc(h: Real): Real = tFall(h) + h / c
// Производная T_calc по h // Находим интервал, содержащий корень
fun dTcalc_dh(h: Real): Real = dtFall_dh(h) + 1.0 / c // Нижняя граница: глубина не может быть отрицательной
var lo = 0.0
// Начальное приближение (без сопротивления) // Верхняя граница: сначала попробуем оценку по свободному падению (без звука)
val term = 1.0 + g * T / c var hi = 0.5 * g * T * T // максимальная глубина, если бы не было сопротивления и звука
val sqrtTerm = sqrt(1.0 + 2.0 * g * T / c) // Уточним hi, чтобы Tcalc(hi) было заведомо больше T
var h = (c * c / g) * (term - sqrtTerm) while (Tcalc(hi) < T && hi < 1e4) {
hi *= 2.0
// Проверка на валидность начального приближения
if (h.isNaN() || h <= 0.0) {
// Если формула дала некорректный результат, используем оценку по свободному падению
h = 0.5 * g * T * T // грубая оценка, всё равно будет уточняться
if (h.isNaN() || h <= 0.0) return null
} }
// Проверка, что hi достаточно велико
if (Tcalc(hi) < T) return null // слишком большая глубина, не укладываемся в разумное
// Итерации Ньютона // Бисекция
var iter = 0 var iter = 0
while (iter < maxIter) { var h = (lo + hi) / 2.0
while (iter < maxIter && (hi - lo) > eps) {
val f = Tcalc(h) - T val f = Tcalc(h) - T
val df = dTcalc_dh(h) if (abs(f) < eps) break
if (f > 0) {
// Если производная близка к нулю, выходим hi = h
if (abs(df) < 1e-12) return null } else {
lo = h
val hNew = h - f / df
// Проверка сходимости
if (abs(hNew - h) < epsilon) {
return hNew
} }
h = (lo + hi) / 2.0
h = hNew
iter++ iter++
println("iter: $iter: $h")
} }
// Не сошлось за maxIter return h
return null
} }
// Пример использования // Пример: T=12 секунд
val T = 6.0 // секунды val T = 26.0
val m = 1.0 // кг val m = 1.0 // кг
val d = 0.1 // м (10 см) val d = 0.1 // м
val depth = calculateDepth(T, m, d) val depth = calculateDepth(T, m, d)
if (depth != null) { if (depth != null) {
println("Глубина: %.2f м"(depth)) println("Глубина: %.2f м"(depth))
} else { // Для проверки выведем теоретическое время при найденной глубине
println("Расчёт не сошёлся") // (можно добавить функцию для самопроверки)
} } else {
println("Расчёт не сошёлся")
}

View File

@ -284,6 +284,31 @@ class DecimalModuleTest {
) )
} }
@Test
fun testDecimalTruncateToTwoFractionDigitsViaGlobalRound() = runTest {
val scope = Script.newScope()
scope.eval(
"""
import lyng.decimal
fun trunc2(x: Decimal): Decimal {
val scaled = x * 100.d
val whole = if (scaled >= 0.d) {
floor(scaled) as Decimal
} else {
ceil(scaled) as Decimal
}
whole / 100.d
}
assertEquals("12.34", trunc2("12.349".d).toStringExpanded())
assertEquals("12.34", trunc2("12.340".d).toStringExpanded())
assertEquals("-12.34", trunc2("-12.349".d).toStringExpanded())
assertEquals("-12.34", trunc2("-12.340".d).toStringExpanded())
""".trimIndent()
)
}
@Test @Test
fun testDecimalMathHelpersFallbackThroughRealTemporarily() = runTest { fun testDecimalMathHelpersFallbackThroughRealTemporarily() = runTest {
val scope = Script.newScope() val scope = Script.newScope()