orbitalkiller: fixed future velocity direction calculation, updated k…

…epler equation solving: number of iterations is not fixed
delorum · Jun 7, 2016 · 80d198b · 80d198b
1 parent 688b168
commit 80d198b
Show file tree

Hide file tree

Showing 2 changed files with 30 additions and 24 deletions.
diff --git a/orbitalkiller/src/main/scala/com/github/dunnololda/scageprojects/orbitalkiller/package.scala b/orbitalkiller/src/main/scala/com/github/dunnololda/scageprojects/orbitalkiller/package.scala
@@ -1477,22 +1477,27 @@ package object orbitalkiller {
     }
 
     private def _iteration(Ex: Double, M: Double): Double = {
-      e * math.sin(Ex) + M
+      M + e * math.sin(Ex)
     }
 
     private val default_num_iterations = 30
 
     // Балк М.Б. Элементы динамики космического полета. Гл. III, параграф 3 "Решение уравнения Кеплера", стр. 111
     // http://pskgu.ru/ebooks/astro3/astro3_03_03.pdf
+    // https://en.wikipedia.org/wiki/Kepler%27s_equation
     def orbitalPointAfterTimeCCW(point1: DVec, time_sec: Long, num_iterations: Int = default_num_iterations): DVec = {
       val t1 = tetaDeg360InPoint(point1)
       val time_from_r_p = travelTimeOnOrbitMsecCCW(0, t1 /*, print_variant = true*/)
       val all_time = time_from_r_p / 1000 + time_sec
       val M = 1 / inv_n * all_time
-      val E7 = (1 to num_iterations).foldLeft(M) {
-        case (res, i) => _iteration(res, M)
+      val E1 = M + e * math.sin(M)
+      def solver(prev_E:Double = E1, i:Int = 0):(Double, Int) = {
+        val diff = (prev_E - e*math.sin(prev_E) - M).abs
+        if(diff < 1E-15) (prev_E, i)
+        else solver(M + e*math.sin(prev_E), i+1)
       }
-      val tg_half_teta_res_rad = math.sqrt((1 + e) / (1 - e)) * math.tan(E7 / 2)
+      val (res_E, i) = solver()
+      val tg_half_teta_res_rad = math.sqrt((1 + e) / (1 - e)) * math.tan(res_E / 2)
       val teta_res_rad = math.atan(tg_half_teta_res_rad) * 2
       val teta_res_deg = teta_res_rad / math.Pi * 180
       orbitalPointByTrueAnomalyDeg(teta_res_deg)
@@ -1503,10 +1508,14 @@ package object orbitalkiller {
       val time_from_r_p_msec = travelTimeOnOrbitMsecCW(0, t1 /*, print_variant = true*/)
       val all_time = t.toLong - (time_from_r_p_msec / 1000 + time_sec)
       val M = 1 / inv_n * all_time
-      val E7 = (1 to num_iterations).foldLeft(M) {
-        case (res, i) => _iteration(res, M)
+      val E1 = M + e * math.sin(M)
+      def solver(prev_E:Double = E1, i:Int = 0):(Double, Int) = {
+        val diff = (prev_E - e*math.sin(prev_E) - M).abs
+        if(diff < 1E-15) (prev_E, i)
+        else solver(M + e*math.sin(prev_E), i+1)
       }
-      val tg_half_teta_res_rad = math.sqrt((1 + e) / (1 - e)) * math.tan(E7 / 2)
+      val (resE, iterations) = solver()
+      val tg_half_teta_res_rad = math.sqrt((1 + e) / (1 - e)) * math.tan(resE / 2)
       val teta_res_rad = math.atan(tg_half_teta_res_rad) * 2
       val teta_res_deg = teta_res_rad / math.Pi * 180
       orbitalPointByTrueAnomalyDeg(teta_res_deg)

diff --git a/...killer/src/main/scala/com/github/dunnololda/scageprojects/orbitalkiller/ships/Ship4.scala b/...killer/src/main/scala/com/github/dunnololda/scageprojects/orbitalkiller/ships/Ship4.scala
@@ -722,6 +722,20 @@ class Ship4(index: Int,
         // current velocity
         drawArrow(DVec.zero.actualPosBeforeRotation, DVec.zero.actualPosBeforeRotation + linearVelocity.n * radius, colorIfPlayerAliveOrRed(BLUE))
         drawArrow(DVec.zero.actualPosBeforeRotation, DVec.zero.actualPosBeforeRotation + relativeLinearVelocity.n * radius, colorIfPlayerAliveOrRed(InterfaceHolder.linearVelocityInfo.color))
+        // velocity direction at stop moment
+        if(_stop_after_number_of_tacts > 0 && InterfaceHolder.orbParams.calculationOn) {
+          thisOrActualProxyShipOrbitData.filter(!_.is_landed).foreach(or => {
+            or.ellipseOrbit.foreach(e => {
+              val time_to_stop_sec = (_stop_after_number_of_tacts * base_dt).toLong
+              val position_at_stop_moment = e.orbitalPointAfterTime(coord, time_to_stop_sec, or.ccw)
+              val (vt, vr) = e.orbitalVelocityInPoint(position_at_stop_moment)
+              val r = if (or.ccw) (position_at_stop_moment - or.planet.coord).n else -(position_at_stop_moment - or.planet.coord).n
+              val t = r.perpendicular
+              val v = (vt*t + vr*r).n
+              drawDashedArrow(DVec.zero.actualPosBeforeRotation, DVec.zero.actualPosBeforeRotation + v*radius, 1, colorIfPlayerAliveOrRed(CYAN))
+            })
+          })
+        }
       }
       if (!InterfaceHolder.sunRelativeInfo.isMinimized) {
         // direction to earth
@@ -740,23 +754,6 @@ class Ship4(index: Int,
           drawArrow(DVec.zero.actualPosBeforeRotation, DVec.zero.actualPosBeforeRotation + (si.monitoring_ship.coord - coord).n * radius, colorIfPlayerAliveOrRed(si.color))
         }
       })
-      if(_stop_after_number_of_tacts > 0 && InterfaceHolder.orbParams.calculationOn) {
-        thisOrActualProxyShipOrbitData match {
-          case Some(or) =>
-            if (!or.is_landed) {
-              or.ellipseOrbit match {
-                case Some(e) =>
-                  val time_to_stop_sec = (_stop_after_number_of_tacts * base_dt).toLong
-                  val position_when_stop_moment = e.orbitalPointAfterTime(coord, time_to_stop_sec, or.ccw)
-                  val r = if (or.ccw) (position_when_stop_moment - or.planet.coord).n else -(position_when_stop_moment - or.planet.coord).n
-                  val t = r.perpendicular
-                  drawDashedArrow(DVec.zero.actualPosBeforeRotation, DVec.zero.actualPosBeforeRotation + t*radius, 1, colorIfPlayerAliveOrRed(CYAN))
-                case None =>
-              }
-            }
-          case None =>
-        }
-      }
     }
 
     /*val pa = (earth.coord - coord).n*(coord.dist(earth.coord) - earth.radius) + (earth.coord - coord).p*70000