orbitalkiller: another attempt to implement algorithm to choose dt fo…

…r real trajectory

orbitalkiller/src/main/scala/com/github/dunnololda/scageprojects/orbitalkiller/OrbitDataUpdater.scala

@@ -130,6 +130,7 @@ object OrbitDataUpdater {
             drawFilledCircle(p, 3 / globalScale, orbit_color)
           })*/
           drawSlidingLines(RealTrajectory2.realTrajectory, ORANGE)
+          drawSlidingLines(RealTrajectory3.realTrajectory, ORANGE_RED)
         }
       }
       /*drawLine(new_o.f*scale, new_o.center*scale, GRAY)
@@ -253,6 +254,7 @@ object OrbitDataUpdater {
             drawFilledCircle(p, 3 / globalScale, orbit_color)
           })*/
           drawSlidingLines(RealTrajectory2.realTrajectory, ORANGE)
+          drawSlidingLines(RealTrajectory3.realTrajectory, ORANGE_RED)
         }
       }
       if(InterfaceHolder.namesSwitcher.showNames) {

orbitalkiller/src/main/scala/com/github/dunnololda/scageprojects/orbitalkiller/OrbitalKiller.scala

@@ -165,12 +165,14 @@ object OrbitalKiller extends ScageScreenAppDMT("Orbital Killer", property("scree
       _update_orbits = true
       RealTrajectory.init()
       RealTrajectory2.init()
+      //RealTrajectory3.init()
     }
   }
 
   actionDynamicPeriodIgnorePause(500 / timeMultiplier) {
     RealTrajectory.continue()
     RealTrajectory2.continue()
+    //RealTrajectory3.continue()
   }
 
   val sun = new Star(

orbitalkiller/src/main/scala/com/github/dunnololda/scageprojects/orbitalkiller/RealTrajectory.scala

@@ -6,18 +6,16 @@ import scala.collection.immutable
 import scala.collection.mutable.ArrayBuffer
 import OrbitalKiller._
 
-object RealTrajectory extends RealTrajectoryC
-object RealTrajectory2 extends RealTrajectoryC {
-  override protected def chooseDt:Double = base_dt
-}
+object RealTrajectory extends RealTrajectoryC(None)
+object RealTrajectory2 extends RealTrajectoryC(Some(1))
+object RealTrajectory3 extends RealTrajectoryC(Some(100))
 
-class RealTrajectoryC {
+class RealTrajectoryC(max_multiplier:Option[Double]) {
   private var real_trajectory:ArrayBuffer[DVec] = ArrayBuffer[DVec]()
   var curPoints:Long = 0
   private var dropped = 0
   private var system_evolution_copy:SystemEvolution = _
   private var celestials:Seq[(CelestialBody, MutableBodyState)] = _
-
   private val angle_diff = 1
 
   def init(): Unit = {
@@ -33,21 +31,48 @@ class RealTrajectoryC {
     }).values.toSeq
   }
 
+  private var min_m:Double = Double.MaxValue
+  private var max_m:Double = 0
+
   protected def chooseDt:Double = {
-    // выберем dt по критерию расстояния корабля до ближайшей поверхности: в 100 км от поверхности объект, движущийся со скоростью 30 км/сек
-    // за dt должен пролететь не более 500 метров (это справедливо для base_dt). Чем дальше, тем dt может быть больше.
-    // но пока работают двигатели, dt должен быть равен base_dt, иначе неверно работают формулы.
-    if(player_ship.engines.exists(_.stopMomentTacts > system_evolution_copy.tacts)) {
-      base_dt
+    if(player_ship.engines.exists(_.stopMomentTacts >= system_evolution_copy.tacts)) {
+      base_dt // пока работают двигатели, dt должен быть равен base_dt, иначе неверно работают формулы.
     } else {
-      system_evolution_copy.bodyState(player_ship.thisOrActualProxyShipIndex).map(bs => {
+      if(max_multiplier.exists(_ == 1)) base_dt
+      else {
+        (for {
+          ss <- system_evolution_copy.bodyState(sun.index)
+          es <- system_evolution_copy.bodyState(earth.index)
+          ms <- system_evolution_copy.bodyState(moon.index)
+          ps <- system_evolution_copy.bodyState(player_ship.index)
+        } yield {
+          // http://arxiv.org/pdf/1105.1082.pdf
+          // N-body simulations of gravitational dynamics, Walter Dehnen and Justin I. Read,
+          // p. 7, 2.2.3 The choise of time-step, formula (20)
+          // Считаем гравитационный потенциал в точке, от модуля извлекаем квадртаный корень и делим на модуль ускорения корабля.
+          // И умножаем на коэффициент. Эмпирически подобрал коэффициент 1/25565. При этом минимальный dt будет base_dt (1/63 секунды),
+          // максимальный в 6000 больше.
+          val x = math.sqrt(G * (/*sun.mass/ss.coord.dist(ps.coord) +
+                             */ earth.mass / es.coord.dist(ps.coord) +
+            moon.mass / ms.coord.dist(ps.coord))) / ps.acc.norma
+          //val m = x/25565/base_dt
+          val m = x / 5600 / base_dt
+          if (m < min_m) min_m = m
+          if (m > max_m) max_m = m
+          println(f"x = $m*base_dt, min_m = $min_m, max_m = $max_m")
+          if(max_multiplier.isEmpty) m * base_dt
+          else math.min(m, max_multiplier.get)*base_dt
+        }).getOrElse(base_dt)
+      }
+      /*val multiplier = system_evolution_copy.bodyState(player_ship.thisOrActualProxyShipIndex).map(bs => {
         val (min_dist_to_any_surface, planet_mass) = celestials.map(x => (bs.coord.dist(x._2.coord) - x._1.radius, x._2.mass)).minBy(_._1)
         if(min_dist_to_any_surface >= 10000000) {
-          math.max(base_dt, 0.005 * min_dist_to_any_surface / bs.vel.norma)
+          math.max(1.0, math.min(30000.0/100000 * min_dist_to_any_surface / bs.vel.norma, max_multiplier.getOrElse(Double.MaxValue)))
         } else {
-          math.max(base_dt, 0.005 * min_dist_to_any_surface / bs.vel.norma * moon.mass / planet_mass)
+          math.max(1.0, math.min(30000.0/100000 * min_dist_to_any_surface / bs.vel.norma * moon.mass / planet_mass, max_multiplier.getOrElse(Double.MaxValue)))
         }
-      }).getOrElse(base_dt)
+      }).getOrElse(1.0)
+      multiplier*base_dt*/
     }
   }
 

orbitalkiller/src/main/scala/com/github/dunnololda/scageprojects/orbitalkiller/package.scala

@@ -1863,8 +1863,8 @@ package object orbitalkiller {
    */
   def calculateOrbit(planet_mass: Double, planet_coord: DVec, body_mass: Double, body_relative_coord: DVec, body_relative_velocity: DVec, G: Double): KeplerOrbit = {
     //https://ru.wikipedia.org/wiki/Гравитационный_параметр
-    val mu = (planet_mass + body_mass) * G // гравитационный параметр
-    //val mu = planet_mass*G // гравитационный параметр
+    //val mu = (planet_mass + body_mass) * G // гравитационный параметр
+    val mu = planet_mass*G // гравитационный параметр
 
     //http://ru.wikipedia.org/wiki/Кеплеровы_элементы_орбиты
     //val a = body_relative_coord.norma*k/(2*k*k - body_relative_coord.norma*body_velocity.norma2)