reverting akima intp: derivatives are not being preserved at the sample points

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()
     }
 }