|
| 1 | +package com.lallimaven.eclipsetimer.wear |
| 2 | + |
| 3 | +import kotlin.math.PI |
| 4 | +import kotlin.math.asin |
| 5 | +import kotlin.math.atan |
| 6 | +import kotlin.math.atan2 |
| 7 | +import kotlin.math.cos |
| 8 | +import kotlin.math.max |
| 9 | +import kotlin.math.sin |
| 10 | +import kotlin.math.tan |
| 11 | + |
| 12 | +object LocalSunMoonCalculator { |
| 13 | + data class LiveGeometry( |
| 14 | + val showMoon: Boolean, |
| 15 | + val moonRadiusNorm: Float, |
| 16 | + val moonCenterXNorm: Float, |
| 17 | + val moonCenterYNorm: Float, |
| 18 | + ) |
| 19 | + |
| 20 | + private data class EquatorialPosition( |
| 21 | + val rightAscensionRad: Double, |
| 22 | + val declinationRad: Double, |
| 23 | + val distanceKm: Double, |
| 24 | + val angularRadiusRad: Double, |
| 25 | + ) |
| 26 | + |
| 27 | + private data class HorizontalPosition( |
| 28 | + val altitudeRad: Double, |
| 29 | + val azimuthRad: Double, |
| 30 | + ) |
| 31 | + |
| 32 | + private data class Vec3( |
| 33 | + val x: Double, |
| 34 | + val y: Double, |
| 35 | + val z: Double, |
| 36 | + ) |
| 37 | + |
| 38 | + fun calculateLiveGeometry( |
| 39 | + latitudeDeg: Double, |
| 40 | + longitudeDeg: Double, |
| 41 | + epochMillis: Long, |
| 42 | + ): LiveGeometry { |
| 43 | + val julianDay = epochMillis / MILLIS_PER_DAY + JULIAN_UNIX_EPOCH |
| 44 | + val daysSinceJ2000 = julianDay - JULIAN_J2000 |
| 45 | + val centuriesSinceJ2000 = daysSinceJ2000 / 36525.0 |
| 46 | + val obliquityRad = toRadians(23.439291 - 0.0130042 * centuriesSinceJ2000) |
| 47 | + |
| 48 | + val sunEquatorial = computeSunEquatorial(daysSinceJ2000, obliquityRad) |
| 49 | + val moonEquatorial = computeMoonEquatorial(daysSinceJ2000, obliquityRad) |
| 50 | + val moonTopocentric = applyMoonTopocentricCorrection( |
| 51 | + moonEquatorial = moonEquatorial, |
| 52 | + latitudeDeg = latitudeDeg, |
| 53 | + longitudeDeg = longitudeDeg, |
| 54 | + julianDay = julianDay, |
| 55 | + ) |
| 56 | + |
| 57 | + val localSiderealRad = computeLocalSiderealRad(julianDay, longitudeDeg) |
| 58 | + val latitudeRad = toRadians(latitudeDeg) |
| 59 | + val sunHorizontal = convertToHorizontal(sunEquatorial, localSiderealRad, latitudeRad) |
| 60 | + val moonHorizontal = convertToHorizontal(moonTopocentric, localSiderealRad, latitudeRad) |
| 61 | + val sunVector = horizontalToEnuVector(sunHorizontal) |
| 62 | + val moonVector = horizontalToEnuVector(moonHorizontal) |
| 63 | + val centerSeparationRad = angularSeparationRad(sunVector, moonVector) |
| 64 | + |
| 65 | + val sunAboveHorizon = sunHorizontal.altitudeRad > MIN_BODY_ALTITUDE_RAD |
| 66 | + val moonAboveHorizon = moonHorizontal.altitudeRad > MIN_BODY_ALTITUDE_RAD |
| 67 | + val overlapThresholdRad = |
| 68 | + sunEquatorial.angularRadiusRad + moonTopocentric.angularRadiusRad + ECLIPSE_MARGIN_RAD |
| 69 | + val shouldShowMoon = sunAboveHorizon && moonAboveHorizon && centerSeparationRad <= overlapThresholdRad |
| 70 | + if (!shouldShowMoon || sunEquatorial.angularRadiusRad <= 0.0) { |
| 71 | + return LiveGeometry( |
| 72 | + showMoon = false, |
| 73 | + moonRadiusNorm = 0f, |
| 74 | + moonCenterXNorm = DEFAULT_CENTER_NORM.toFloat(), |
| 75 | + moonCenterYNorm = DEFAULT_CENTER_NORM.toFloat(), |
| 76 | + ) |
| 77 | + } |
| 78 | + |
| 79 | + val tangentBasisEast = Vec3( |
| 80 | + x = cos(sunHorizontal.azimuthRad), |
| 81 | + y = -sin(sunHorizontal.azimuthRad), |
| 82 | + z = 0.0, |
| 83 | + ) |
| 84 | + val tangentBasisNorth = Vec3( |
| 85 | + x = -sin(sunHorizontal.altitudeRad) * sin(sunHorizontal.azimuthRad), |
| 86 | + y = -sin(sunHorizontal.altitudeRad) * cos(sunHorizontal.azimuthRad), |
| 87 | + z = cos(sunHorizontal.altitudeRad), |
| 88 | + ) |
| 89 | + |
| 90 | + val forward = max(dot(moonVector, sunVector), MIN_FORWARD_DOT) |
| 91 | + val offsetEastRad = dot(moonVector, tangentBasisEast) / forward |
| 92 | + val offsetNorthRad = dot(moonVector, tangentBasisNorth) / forward |
| 93 | + val normPerRad = DEFAULT_SUN_RADIUS_NORM / sunEquatorial.angularRadiusRad |
| 94 | + val moonRadiusNorm = |
| 95 | + (DEFAULT_SUN_RADIUS_NORM * (moonTopocentric.angularRadiusRad / sunEquatorial.angularRadiusRad)) |
| 96 | + .coerceIn(MIN_MOON_RADIUS_NORM, MAX_MOON_RADIUS_NORM) |
| 97 | + val moonCenterXNorm = (DEFAULT_CENTER_NORM + offsetEastRad * normPerRad).coerceIn(0.0, 1.0) |
| 98 | + val moonCenterYNorm = (DEFAULT_CENTER_NORM - offsetNorthRad * normPerRad).coerceIn(0.0, 1.0) |
| 99 | + |
| 100 | + return LiveGeometry( |
| 101 | + showMoon = true, |
| 102 | + moonRadiusNorm = moonRadiusNorm.toFloat(), |
| 103 | + moonCenterXNorm = moonCenterXNorm.toFloat(), |
| 104 | + moonCenterYNorm = moonCenterYNorm.toFloat(), |
| 105 | + ) |
| 106 | + } |
| 107 | + |
| 108 | + private fun computeSunEquatorial( |
| 109 | + daysSinceJ2000: Double, |
| 110 | + obliquityRad: Double, |
| 111 | + ): EquatorialPosition { |
| 112 | + val meanAnomalyRad = toRadians(normalizeDegrees(357.5291 + 0.98560028 * daysSinceJ2000)) |
| 113 | + val equationOfCenterRad = |
| 114 | + toRadians(1.9148) * sin(meanAnomalyRad) + |
| 115 | + toRadians(0.0200) * sin(2.0 * meanAnomalyRad) + |
| 116 | + toRadians(0.0003) * sin(3.0 * meanAnomalyRad) |
| 117 | + val perihelionRad = toRadians(102.9372) |
| 118 | + val eclipticLongitudeRad = meanAnomalyRad + equationOfCenterRad + perihelionRad + PI |
| 119 | + |
| 120 | + val rightAscensionRad = atan2( |
| 121 | + sin(eclipticLongitudeRad) * cos(obliquityRad), |
| 122 | + cos(eclipticLongitudeRad), |
| 123 | + ) |
| 124 | + val declinationRad = asin(sin(eclipticLongitudeRad) * sin(obliquityRad)) |
| 125 | + val distanceAu = |
| 126 | + 1.00014 - |
| 127 | + 0.01671 * cos(meanAnomalyRad) - |
| 128 | + 0.00014 * cos(2.0 * meanAnomalyRad) |
| 129 | + val angularRadiusRad = toRadians(0.2666 / distanceAu) |
| 130 | + |
| 131 | + return EquatorialPosition( |
| 132 | + rightAscensionRad = rightAscensionRad, |
| 133 | + declinationRad = declinationRad, |
| 134 | + distanceKm = distanceAu * ASTRONOMICAL_UNIT_KM, |
| 135 | + angularRadiusRad = angularRadiusRad, |
| 136 | + ) |
| 137 | + } |
| 138 | + |
| 139 | + private fun computeMoonEquatorial( |
| 140 | + daysSinceJ2000: Double, |
| 141 | + obliquityRad: Double, |
| 142 | + ): EquatorialPosition { |
| 143 | + val meanLongitudeRad = toRadians(normalizeDegrees(218.316 + 13.176396 * daysSinceJ2000)) |
| 144 | + val meanAnomalyRad = toRadians(normalizeDegrees(134.963 + 13.064993 * daysSinceJ2000)) |
| 145 | + val argumentOfLatitudeRad = toRadians(normalizeDegrees(93.272 + 13.229350 * daysSinceJ2000)) |
| 146 | + val eclipticLongitudeRad = meanLongitudeRad + toRadians(6.289) * sin(meanAnomalyRad) |
| 147 | + val eclipticLatitudeRad = toRadians(5.128) * sin(argumentOfLatitudeRad) |
| 148 | + val distanceKm = 385_001.0 - 20_905.0 * cos(meanAnomalyRad) |
| 149 | + |
| 150 | + val rightAscensionRad = atan2( |
| 151 | + sin(eclipticLongitudeRad) * cos(obliquityRad) - |
| 152 | + tan(eclipticLatitudeRad) * sin(obliquityRad), |
| 153 | + cos(eclipticLongitudeRad), |
| 154 | + ) |
| 155 | + val declinationRad = asin( |
| 156 | + sin(eclipticLatitudeRad) * cos(obliquityRad) + |
| 157 | + cos(eclipticLatitudeRad) * sin(obliquityRad) * sin(eclipticLongitudeRad), |
| 158 | + ) |
| 159 | + val angularRadiusRad = asin((MOON_RADIUS_KM / distanceKm).coerceIn(-1.0, 1.0)) |
| 160 | + |
| 161 | + return EquatorialPosition( |
| 162 | + rightAscensionRad = rightAscensionRad, |
| 163 | + declinationRad = declinationRad, |
| 164 | + distanceKm = distanceKm, |
| 165 | + angularRadiusRad = angularRadiusRad, |
| 166 | + ) |
| 167 | + } |
| 168 | + |
| 169 | + private fun applyMoonTopocentricCorrection( |
| 170 | + moonEquatorial: EquatorialPosition, |
| 171 | + latitudeDeg: Double, |
| 172 | + longitudeDeg: Double, |
| 173 | + julianDay: Double, |
| 174 | + ): EquatorialPosition { |
| 175 | + val latitudeRad = toRadians(latitudeDeg) |
| 176 | + val localSiderealRad = computeLocalSiderealRad(julianDay, longitudeDeg) |
| 177 | + val hourAngleRad = normalizeSignedRadians(localSiderealRad - moonEquatorial.rightAscensionRad) |
| 178 | + val u = atan(0.99664719 * tan(latitudeRad)) |
| 179 | + val observerX = cos(u) |
| 180 | + val observerY = 0.99664719 * sin(u) |
| 181 | + val sinParallax = EARTH_RADIUS_KM / moonEquatorial.distanceKm |
| 182 | + val rightAscensionCorrectionRad = atan2( |
| 183 | + -observerX * sinParallax * sin(hourAngleRad), |
| 184 | + cos(moonEquatorial.declinationRad) - observerX * sinParallax * cos(hourAngleRad), |
| 185 | + ) |
| 186 | + val correctedRightAscensionRad = moonEquatorial.rightAscensionRad + rightAscensionCorrectionRad |
| 187 | + val correctedDeclinationRad = atan2( |
| 188 | + (sin(moonEquatorial.declinationRad) - observerY * sinParallax) * cos(rightAscensionCorrectionRad), |
| 189 | + cos(moonEquatorial.declinationRad) - observerX * sinParallax * cos(hourAngleRad), |
| 190 | + ) |
| 191 | + |
| 192 | + return moonEquatorial.copy( |
| 193 | + rightAscensionRad = correctedRightAscensionRad, |
| 194 | + declinationRad = correctedDeclinationRad, |
| 195 | + ) |
| 196 | + } |
| 197 | + |
| 198 | + private fun computeLocalSiderealRad(julianDay: Double, longitudeDeg: Double): Double { |
| 199 | + val centuriesSinceJ2000 = (julianDay - JULIAN_J2000) / 36525.0 |
| 200 | + val gmstDeg = |
| 201 | + 280.46061837 + |
| 202 | + 360.98564736629 * (julianDay - JULIAN_J2000) + |
| 203 | + 0.000387933 * centuriesSinceJ2000 * centuriesSinceJ2000 - |
| 204 | + (centuriesSinceJ2000 * centuriesSinceJ2000 * centuriesSinceJ2000) / 38710000.0 |
| 205 | + return normalizeRadians(toRadians(gmstDeg + longitudeDeg)) |
| 206 | + } |
| 207 | + |
| 208 | + private fun convertToHorizontal( |
| 209 | + equatorial: EquatorialPosition, |
| 210 | + localSiderealRad: Double, |
| 211 | + latitudeRad: Double, |
| 212 | + ): HorizontalPosition { |
| 213 | + val hourAngleRad = normalizeSignedRadians(localSiderealRad - equatorial.rightAscensionRad) |
| 214 | + val sinAltitude = |
| 215 | + sin(latitudeRad) * sin(equatorial.declinationRad) + |
| 216 | + cos(latitudeRad) * cos(equatorial.declinationRad) * cos(hourAngleRad) |
| 217 | + val altitudeRad = asin(sinAltitude.coerceIn(-1.0, 1.0)) |
| 218 | + val azimuthRad = atan2( |
| 219 | + -sin(hourAngleRad), |
| 220 | + tan(equatorial.declinationRad) * cos(latitudeRad) - sin(latitudeRad) * cos(hourAngleRad), |
| 221 | + ) |
| 222 | + return HorizontalPosition( |
| 223 | + altitudeRad = altitudeRad, |
| 224 | + azimuthRad = normalizeRadians(azimuthRad), |
| 225 | + ) |
| 226 | + } |
| 227 | + |
| 228 | + private fun horizontalToEnuVector(horizontal: HorizontalPosition): Vec3 { |
| 229 | + val cosAltitude = cos(horizontal.altitudeRad) |
| 230 | + return Vec3( |
| 231 | + x = cosAltitude * sin(horizontal.azimuthRad), |
| 232 | + y = cosAltitude * cos(horizontal.azimuthRad), |
| 233 | + z = sin(horizontal.altitudeRad), |
| 234 | + ) |
| 235 | + } |
| 236 | + |
| 237 | + private fun angularSeparationRad(first: Vec3, second: Vec3): Double { |
| 238 | + val cosine = dot(first, second).coerceIn(-1.0, 1.0) |
| 239 | + return kotlin.math.acos(cosine) |
| 240 | + } |
| 241 | + |
| 242 | + private fun dot(first: Vec3, second: Vec3): Double = |
| 243 | + first.x * second.x + first.y * second.y + first.z * second.z |
| 244 | + |
| 245 | + private fun normalizeRadians(value: Double): Double { |
| 246 | + val period = 2.0 * PI |
| 247 | + var normalized = value % period |
| 248 | + if (normalized < 0.0) { |
| 249 | + normalized += period |
| 250 | + } |
| 251 | + return normalized |
| 252 | + } |
| 253 | + |
| 254 | + private fun normalizeSignedRadians(value: Double): Double { |
| 255 | + var normalized = normalizeRadians(value) |
| 256 | + if (normalized > PI) { |
| 257 | + normalized -= 2.0 * PI |
| 258 | + } |
| 259 | + return normalized |
| 260 | + } |
| 261 | + |
| 262 | + private fun normalizeDegrees(value: Double): Double { |
| 263 | + var normalized = value % 360.0 |
| 264 | + if (normalized < 0.0) { |
| 265 | + normalized += 360.0 |
| 266 | + } |
| 267 | + return normalized |
| 268 | + } |
| 269 | + |
| 270 | + private fun toRadians(valueDeg: Double): Double = valueDeg * PI / 180.0 |
| 271 | + |
| 272 | + private const val MILLIS_PER_DAY = 86_400_000.0 |
| 273 | + private const val JULIAN_UNIX_EPOCH = 2_440_587.5 |
| 274 | + private const val JULIAN_J2000 = 2_451_545.0 |
| 275 | + private const val ASTRONOMICAL_UNIT_KM = 149_597_870.7 |
| 276 | + private const val EARTH_RADIUS_KM = 6_378.14 |
| 277 | + private const val MOON_RADIUS_KM = 1_737.4 |
| 278 | + private const val DEFAULT_SUN_RADIUS_NORM = 0.24 |
| 279 | + private const val DEFAULT_CENTER_NORM = 0.5 |
| 280 | + private const val MIN_MOON_RADIUS_NORM = 0.05 |
| 281 | + private const val MAX_MOON_RADIUS_NORM = 0.45 |
| 282 | + private const val MIN_FORWARD_DOT = 0.000001 |
| 283 | + private const val MIN_BODY_ALTITUDE_RAD = -1.5 * PI / 180.0 |
| 284 | + private const val ECLIPSE_MARGIN_RAD = 0.08 * PI / 180.0 |
| 285 | +} |
0 commit comments