Improve tide corrections

README.md

@@ -19,3 +19,4 @@ Astronomy algorithms are a combination of custom made algorithms and implementat
 
 - [osgb](https://github.com/kylecorry31/osgb/blob/master/LICENSE)
 - [Geo-Coordinate-Conversion-Java](https://github.com/kylecorry31/Geo-Coordinate-Conversion-Java/blob/master/GDAL_License.TXT)
+- Tide corrections are computed using [uptide](https://github.com/stephankramer/uptide/blob/master/LICENSE)

src/main/kotlin/com/kylecorry/sol/science/oceanography/TideConstituentCorrections.kt

@@ -0,0 +1,110 @@
+package com.kylecorry.sol.science.oceanography
+
+import com.kylecorry.sol.math.SolMath.cube
+import com.kylecorry.sol.math.SolMath.power
+import com.kylecorry.sol.math.SolMath.square
+
+class ConstituentCorrection(val f: Float, val uv: Float)
+
+class TideConstituentCorrections(
+    private val Q1: ConstituentCorrection,
+    private val O1: ConstituentCorrection,
+    private val P1: ConstituentCorrection,
+    private val K1: ConstituentCorrection,
+    private val N2: ConstituentCorrection,
+    private val M2: ConstituentCorrection,
+    private val S2: ConstituentCorrection,
+    private val K2: ConstituentCorrection,
+    private val L2: ConstituentCorrection,
+    private val M1: ConstituentCorrection,
+    private val J1: ConstituentCorrection,
+    private val MM: ConstituentCorrection,
+    private val MF: ConstituentCorrection,
+    private val M3: ConstituentCorrection
+) {
+
+    fun getUV(constituent: TideConstituent): Float {
+        return when (constituent) {
+            TideConstituent.M2 -> M2.uv
+            TideConstituent.S2 -> S2.uv
+            TideConstituent.N2 -> N2.uv
+            TideConstituent.K1 -> K1.uv
+            TideConstituent.M4 -> M2.uv * 2
+            TideConstituent.O1 -> O1.uv
+            TideConstituent.M6 -> M2.uv * 3
+            TideConstituent.MK3 -> M2.uv + K1.uv
+            TideConstituent.S4 -> S2.uv * 2
+            TideConstituent.MN4 -> M2.uv + N2.uv
+            TideConstituent.NU2 -> M2.uv
+            TideConstituent.S6 -> S2.uv * 3
+            TideConstituent.MU2 -> M2.uv * 2 - S2.uv
+            TideConstituent._2N2 -> M2.uv
+            TideConstituent.LAM2 -> M2.uv
+            TideConstituent.S1 -> 0f
+            TideConstituent.M1 -> M1.uv
+            TideConstituent.J1 -> J1.uv
+            TideConstituent.MM -> MM.uv
+            TideConstituent.SSA -> 0f
+            TideConstituent.SA -> 0f
+            TideConstituent.MSF -> S2.uv - M2.uv
+            TideConstituent.MF -> MF.uv
+            TideConstituent.RHO -> M2.uv - K1.uv // NU2 - K1
+            TideConstituent.Q1 -> Q1.uv
+            TideConstituent.T2 -> 0f
+            TideConstituent.R2 -> 0f
+            TideConstituent._2Q1 -> N2.uv - J1.uv
+            TideConstituent.P1 -> P1.uv
+            TideConstituent._2SM2 -> 2 * S2.uv - M2.uv
+            TideConstituent.M3 -> M3.uv
+            TideConstituent.L2 -> L2.uv
+            TideConstituent._2MK3 -> M2.uv + O1.uv
+            TideConstituent.K2 -> K2.uv
+            TideConstituent.M8 -> M2.uv * 4
+            TideConstituent.MS4 -> M2.uv + S2.uv
+            TideConstituent.Z0 -> 0f
+        }
+    }
+
+    fun getF(constituent: TideConstituent): Float {
+        return when (constituent) {
+            TideConstituent.M2 -> M2.f
+            TideConstituent.S2 -> S2.f
+            TideConstituent.N2 -> N2.f
+            TideConstituent.K1 -> K1.f
+            TideConstituent.M4 -> square(M2.f)
+            TideConstituent.O1 -> O1.f
+            TideConstituent.M6 -> cube(M2.f)
+            TideConstituent.MK3 -> M2.f * K1.f
+            TideConstituent.S4 -> square(S2.f)
+            TideConstituent.MN4 -> M2.f * N2.f
+            TideConstituent.NU2 -> M2.f
+            TideConstituent.S6 -> cube(S2.f)
+            TideConstituent.MU2 -> square(M2.f) * S2.f
+            TideConstituent._2N2 -> M2.f
+            TideConstituent.LAM2 -> M2.f
+            TideConstituent.S1 -> 1f
+            TideConstituent.M1 -> M1.f
+            TideConstituent.J1 -> J1.f
+            TideConstituent.MM -> MM.f
+            TideConstituent.SSA -> 1f
+            TideConstituent.SA -> 1f
+            TideConstituent.MSF -> S2.f * M2.f
+            TideConstituent.MF -> MF.f
+            TideConstituent.RHO -> M2.f * K1.f // NU2 * K1
+            TideConstituent.Q1 -> Q1.f
+            TideConstituent.T2 -> 1f
+            TideConstituent.R2 -> 1f
+            TideConstituent._2Q1 -> N2.f * J1.f
+            TideConstituent.P1 -> P1.f
+            TideConstituent._2SM2 -> square(S2.f) * M2.f
+            TideConstituent.M3 -> M3.f
+            TideConstituent.L2 -> L2.f
+            TideConstituent._2MK3 -> M2.f * O1.f
+            TideConstituent.K2 -> K2.f
+            TideConstituent.M8 -> power(M2.f, 4)
+            TideConstituent.MS4 -> M2.f * S2.f
+            TideConstituent.Z0 -> 1f
+        }
+    }
+
+}

src/main/kotlin/com/kylecorry/sol/science/oceanography/waterlevel/HarmonicWaterLevelCalculator.kt

@@ -1,34 +1,87 @@
 package com.kylecorry.sol.science.oceanography.waterlevel
 
 import com.kylecorry.sol.math.SolMath
-import com.kylecorry.sol.science.oceanography.AstronomicalArgumentCalculator
-import com.kylecorry.sol.science.oceanography.NodeFactorCalculator
+import com.kylecorry.sol.science.oceanography.ConstituentCorrection
 import com.kylecorry.sol.science.oceanography.TidalHarmonic
-import com.kylecorry.sol.science.oceanography.TideConstituent
-import com.kylecorry.sol.time.Time
+import com.kylecorry.sol.science.oceanography.TideConstituentCorrections
 import java.time.Duration
 import java.time.LocalDateTime
 import java.time.ZoneId
 import java.time.ZonedDateTime
+import kotlin.math.abs
 
 class HarmonicWaterLevelCalculator(private val harmonics: List<TidalHarmonic>) :
     IWaterLevelCalculator {
+
     override fun calculate(time: ZonedDateTime): Float {
-        // TODO: This is a hack to get the correct year for the astronomical arguments
-        val year = 2022
-        val start = ZonedDateTime.of(
-            LocalDateTime.of(year, 1, 1, 0, 0),
-            ZoneId.of("UTC")
-        )
+        val start = getStartDate(time.year)
+        val corrections = getCorrections(time.year)
         val t = Duration.between(start, time).seconds / 3600f
         val heights = harmonics.map {
-            val constituentPhase = AstronomicalArgumentCalculator.get(it.constituent, year)
-            val nodeFactor = NodeFactorCalculator.get(
-                it.constituent,
-                year
+            corrections.getF(it.constituent) * it.amplitude * SolMath.cosDegrees(
+                it.constituent.speed * t + corrections.getUV(
+                    it.constituent
+                ) - it.phase
             )
-            nodeFactor * it.amplitude * SolMath.cosDegrees(it.constituent.speed * t + constituentPhase - it.phase)
         }
         return heights.sum()
     }
+
+    companion object {
+
+        private fun getStartDate(year: Int): ZonedDateTime {
+            return ZonedDateTime.of(
+                LocalDateTime.of(getClosestYearWithData(year), 1, 1, 0, 0),
+                ZoneId.of("UTC")
+            )
+        }
+
+        private fun getCorrections(year: Int): TideConstituentCorrections {
+            return when (getClosestYearWithData(year)) {
+                2024 -> corrections2024
+                2025 -> corrections2025
+                else -> corrections2025
+            }
+        }
+
+        private fun getClosestYearWithData(year: Int): Int {
+            val yearsWithData = listOf(2024, 2025)
+            return yearsWithData.minByOrNull { abs(it - year) } ?: yearsWithData.first()
+        }
+
+        private val corrections2024 = TideConstituentCorrections(
+            ConstituentCorrection(1.1837230626740127f, 65.7634306026157f),
+            ConstituentCorrection(1.1837230626740127f, 241.93791654123925f),
+            ConstituentCorrection(1.0f, 349.84105553831614f),
+            ConstituentCorrection(1.1134500118048742f, 6.987321819855424f),
+            ConstituentCorrection(0.9654291266366927f, 71.3253040141426f),
+            ConstituentCorrection(0.9654291266366927f, 247.49978995299898f),
+            ConstituentCorrection(1.0f, 0.0f),
+            ConstituentCorrection(1.2912235076190786f, 194.01027966105903f),
+            ConstituentCorrection(0.9654291266366927f, 243.67427589156432f),
+            ConstituentCorrection(0.9827145633183463f, 123.74989497649949f),
+            ConstituentCorrection(1.1539047999567282f, 181.72210479492787f),
+            ConstituentCorrection(0.8785347692640552f, 176.17448593885638f),
+            ConstituentCorrection(1.4298200424975473f, 303.6239568705205f),
+            ConstituentCorrection(0.948143689955039f, 191.24968492938206f),
+        )
+
+        private val corrections2025 = TideConstituentCorrections(
+            ConstituentCorrection(1.195936272630187f, 36.03341654129326f),
+            ConstituentCorrection(1.195936272630187f, 313.99553812877275f),
+            ConstituentCorrection(1.0f, 349.0941103290388f),
+            ConstituentCorrection(1.1209608093715056f, 10.673557053567492f),
+            ConstituentCorrection(0.963012609158733f, 46.602554438402876f),
+            ConstituentCorrection(0.963012609158733f, 324.56467602588236f),
+            ConstituentCorrection(1.0f, 0.0f),
+            ConstituentCorrection(1.3099025346108744f, 201.34972457219556f),
+            ConstituentCorrection(0.963012609158733f, 62.52679761359468f),
+            ConstituentCorrection(0.9815063045793665f, 342.2823380129412f),
+            ConstituentCorrection(1.1649424068155168f, 288.53021528816316f),
+            ConstituentCorrection(0.8700443024496025f, 277.9621215874795f),
+            ConstituentCorrection(1.4568589137374195f, 236.57359976903535f),
+            ConstituentCorrection(0.9445189137380996f, 126.84701403905638f),
+        )
+    }
+
 }