From 6a78cf6a412412235f65f561dfce254ab354eb07 Mon Sep 17 00:00:00 2001 From: minjaesong Date: Sat, 2 Sep 2023 23:43:50 +0900 Subject: [PATCH] reverting akima intp: derivatives are not being preserved at the sample points --- .../torvald/terrarum/weather/Weatherbox.kt | 167 +++++++++++------- 1 file changed, 106 insertions(+), 61 deletions(-) diff --git a/src/net/torvald/terrarum/weather/Weatherbox.kt b/src/net/torvald/terrarum/weather/Weatherbox.kt index be831207a..1761b30ba 100644 --- a/src/net/torvald/terrarum/weather/Weatherbox.kt +++ b/src/net/torvald/terrarum/weather/Weatherbox.kt @@ -1,15 +1,13 @@ package net.torvald.terrarum.weather import com.jme3.math.FastMath -import net.torvald.terrarum.App.printdbg import net.torvald.terrarum.floorToInt import net.torvald.terrarum.gameworld.GameWorld import net.torvald.terrarum.gameworld.fmod -import org.apache.commons.math3.analysis.interpolation.AkimaSplineInterpolator -import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction +import org.apache.commons.math3.analysis.UnivariateFunction +import org.apache.commons.math3.analysis.interpolation.NevilleInterpolator +import org.apache.commons.math3.analysis.interpolation.SplineInterpolator import java.util.* -import kotlin.math.max -import kotlin.math.pow data class WeatherSchedule(val weather: BaseModularWeather = WeatherMixer.DEFAULT_WEATHER, val duration: Long = 3600) @@ -93,22 +91,37 @@ open class WeatherStateBox( // pM1 and p4 only exists for the sake of better weather forecasting // - removing p4 and beyond: for faster response to the changing weather schedule and make the forecasting less accurate like irl ) { - - protected lateinit var polynomial: PolynomialSplineFunction - protected val interpolator = AkimaSplineInterpolator() +// protected lateinit var polynomial: UnivariateFunction +// protected val interpolator = SplineInterpolator() open val value: Float; get() = valueAt(x) - open fun valueAt(x: Float) = polynomial.value(x + 1.0).toFloat() + open fun valueAt(x: Float) = when (x.floorToInt()) { + -2 -> interpolate(x + 2, pM2,pM1, p0, p1) + -1 -> interpolate(x + 1, pM1, p0, p1, p2) + 0 -> interpolate(x - 0, p0, p1, p2, p3) + 1 -> interpolate(x - 1, p1, p2, p3, p3) + 2 -> interpolate(x - 2, p2, p3, p3, p3) + 3 -> interpolate(x - 3, p3, p3, p3, p3) + else -> throw IllegalArgumentException() + } + + open protected fun extrapolate(): Float { + val d1 = p2 - p1 + val d2 = p3 - p2 + + // if two slopes are monotonic + return if (d1 * d2 >= 0f) + (d2 + d1) / 2f + p3 + else + (d2 + d1) / 2f + p3 + } open protected fun updatePolynomials() { - polynomial = interpolator.interpolate( - doubleArrayOf(-2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0), - doubleArrayOf(pM2.toDouble(), pM1.toDouble(), p0.toDouble(), p1.toDouble(), p2.toDouble(), p3.toDouble(), p3.toDouble()) - ) + } open fun update(xdelta: Float, next: () -> Float) { - if (!::polynomial.isInitialized) updatePolynomials() +// if (!::polynomial.isInitialized) updatePolynomials() synchronized(WeatherMixer.RNG) { x += xdelta @@ -121,9 +134,21 @@ open class WeatherStateBox( p2 = p3 p3 = next() } - updatePolynomials() +// updatePolynomials() } } + + protected fun interpolate(u: Float, p0: Float, p1: Float, p2: Float, p3: Float): Float { + val T = FastMath.interpolateLinear(u, p1, p2).div(maxOf(p0, p3).coerceAtLeast(1f)).toDouble().coerceIn(0.0, 0.5) +// if (u == x) printdbg(this, "u=$u, p1=$p1, p2=$p2; T=$T") + + val c1 = p1.toDouble() + val c2 = -1.0 * T * p0 + T * p2 + val c3 = 2 * T * p0 + (T - 3) * p1 + (3 - 2 * T) * p2 + -T * p3 + val c4 = -T * p0 + (2 - T) * p1 + (T - 2) * p2 + T * p3 + + return (((c4 * u + c3) * u + c2) * u + c1).toFloat() + } } /** @@ -138,9 +163,8 @@ class WeatherDirBox( p2: Float = 0f, p3: Float = 0f, ) : WeatherStateBox(x, pM2, pM1, p0, p1, p2, p3) { - override fun valueAt(x: Float) = polynomial.value(x + 1.0).plus(2.0).fmod(4.0).minus(2.0).toFloat() - override fun updatePolynomials() { + override fun valueAt(x: Float): Float { var pM2 = pM2 var pM1 = pM1 var p0 = p0 @@ -148,50 +172,55 @@ class WeatherDirBox( var p2 = p2 var p3 = p3 - - if (pM1 - pM2 > 2f) { - pM2 -= 4f - pM1 -= 4f - p0 -= 4f - p1 -= 4f - p2 -= 4f - p3 -= 4f - } - else if (pM1 - pM2 < -2f) { - pM2 += 4f - pM1 += 4f - p0 += 4f - p1 += 4f - p2 += 4f - p3 += 4f + if (x < -2f) { + if (pM1 - pM2 > 2f) { + pM2 -= 4f + pM1 -= 4f + p0 -= 4f + p1 -= 4f + p2 -= 4f + p3 -= 4f + } + else if (pM1 - pM2 < -2f) { + pM2 += 4f + pM1 += 4f + p0 += 4f + p1 += 4f + p2 += 4f + p3 += 4f + } } - if (pM1 - pM2 > 2f) { - pM1 -= 4f - p0 -= 4f - p1 -= 4f - p2 -= 4f - p3 -= 4f - } - else if (pM1 - pM2 < -2f) { - pM1 += 4f - p0 += 4f - p1 += 4f - p2 += 4f - p3 += 4f + if (x < -1f) { + if (pM1 - pM2 > 2f) { + pM1 -= 4f + p0 -= 4f + p1 -= 4f + p2 -= 4f + p3 -= 4f + } + else if (pM1 - pM2 < -2f) { + pM1 += 4f + p0 += 4f + p1 += 4f + p2 += 4f + p3 += 4f + } } - if (p0 - pM1 > 2f) { - p0 -= 4f - p1 -= 4f - p2 -= 4f - p3 -= 4f - } - else if (p0 - pM1 < -2f) { - p0 += 4f - p1 += 4f - p2 += 4f - p3 += 4f + if (x < 0f) { + if (p0 - pM1 > 2f) { + p0 -= 4f + p1 -= 4f + p2 -= 4f + p3 -= 4f + } + else if (p0 - pM1 < -2f) { + p0 += 4f + p1 += 4f + p2 += 4f + p3 += 4f + } } if (p1 - p0 > 2f) { @@ -221,9 +250,25 @@ class WeatherDirBox( p3 += 4f } - polynomial = interpolator.interpolate( - doubleArrayOf(-2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0), - doubleArrayOf(pM2.toDouble(), pM1.toDouble(), p0.toDouble(), p1.toDouble(), p2.toDouble(), p3.toDouble(), p3.toDouble()) - ) + return when (x.floorToInt()) { + -2 -> interpolate2(x + 2, pM2,pM1, p0, p1) + -1 -> interpolate2(x + 1, pM1, p0, p1, p2) + 0 -> interpolate2(x - 0, p0, p1, p2, p3) + 1 -> interpolate2(x - 1, p1, p2, p3, p3) + 2 -> interpolate2(x - 2, p2, p3, p3, p3) + 3 -> interpolate2(x - 3, p3, p3, p3, p3) + else -> throw IllegalArgumentException() + }.plus(2f).fmod(4f).minus(2f) + } + + private fun interpolate2(u: Float, p0: Float, p1: Float, p2: Float, p3: Float): Float { + val T = 0.5 + + val c1 = p1.toDouble() + val c2 = -1.0 * T * p0 + T * p2 + val c3 = 2 * T * p0 + (T - 3) * p1 + (3 - 2 * T) * p2 + -T * p3 + val c4 = -T * p0 + (2 - T) * p1 + (T - 2) * p2 + T * p3 + + return (((c4 * u + c3) * u + c2) * u + c1).toFloat() } } \ No newline at end of file