Merge pull request #30 from ChristopherRabotin/issue-29

Fixed the SOI

src/cmd/continuous/conditions.go

@@ -19,15 +19,15 @@ func InitialEarthOrbit() *dynamics.Orbit {
 
 // FromEarthWaypoints returns the waypoints.
 func FromEarthWaypoints() []dynamics.Waypoint {
-	ref2Sun := &dynamics.WaypointAction{Type: dynamics.REFSUN, Cargo: nil}
+	//ref2Sun := &dynamics.WaypointAction{Type: dynamics.REFSUN, Cargo: nil}
 	ref2Mars := &dynamics.WaypointAction{Type: dynamics.REFMARS, Cargo: nil}
 	return []dynamics.Waypoint{
 		// Loiter for 12 hours (eg. IOT)
 		dynamics.NewLoiter(time.Duration(12)*time.Hour, nil),
 		// Change the inclination by 12 degrees
 		//dynamics.NewRelativeOrbitTarget(nil, []dynamics.RelativeOE{dynamics.RelativeOE{Law: dynamics.OptiΔiCL, Value: 12.0}}),
 		// Leave Earth
-		dynamics.NewOutwardSpiral(dynamics.Earth, ref2Sun),
+		dynamics.NewOutwardSpiral(dynamics.Earth, nil),
 		dynamics.NewLoiter(time.Duration(12)*time.Hour, nil),
 		// Go straight to Mars destination
 		dynamics.NewPlanetTarget(dynamics.Mars, time.Date(2016, 3+7, 14, 9, 31, 0, 0, time.UTC), ref2Mars),

src/dynamics/mission.go

@@ -7,6 +7,7 @@ import (
 	"time"
 
 	"github.com/ChristopherRabotin/ode"
+	"github.com/gonum/floats"
 )
 
 const (
@@ -172,14 +173,6 @@ func (a *Mission) SetState(t float64, s []float64) {
 	a.Orbit.Ω = s[3]
 	a.Orbit.ω = s[4]
 	a.Orbit.ν = s[5]
-	// Let's execute any function which is in the queue of this time step.
-	for _, f := range a.Vehicle.FuncQ {
-		if f == nil {
-			continue
-		}
-		f()
-	}
-	a.Vehicle.FuncQ = make([]func(), 5) // Clear the queue.
 
 	// Orbit sanity checks and warnings.
 	if !a.collided && a.Orbit.GetRNorm() < a.Orbit.Origin.Radius {
@@ -189,15 +182,26 @@ func (a *Mission) SetState(t float64, s []float64) {
 		// Now further from the 1% dead zone
 		a.collided = false
 		a.Vehicle.logger.Log("level", "critical", "subsys", "astro", "revived", a.Orbit.Origin.Name, "dt", a.CurrentDT)
-	} else if a.Orbit.GetRNorm() > a.Orbit.Origin.SOI {
-		a.Vehicle.ToXCentric(Sun, a.CurrentDT, a.Orbit)
+	} else if (a.Orbit.GetRNorm() > a.Orbit.Origin.SOI || floats.EqualWithinAbs(a.Orbit.e, 1, eccentricityε)) && !a.Orbit.Origin.Equals(Sun) {
+		fmt.Printf("%s => %v\n", a.Orbit.Origin, !a.Orbit.Origin.Equals(Sun))
+		a.Vehicle.FuncQ = append(a.Vehicle.FuncQ, a.Vehicle.ToXCentric(Sun, a.CurrentDT, a.Orbit))
 	}
 
 	// Propulsion sanity check
 	if a.Vehicle.FuelMass > 0 && s[6] <= 0 {
 		a.Vehicle.logger.Log("level", "critical", "subsys", "prop", "fuel(kg)", s[6])
 	}
 	a.Vehicle.FuelMass = s[6]
+
+	// Let's execute any function which is in the queue of this time step.
+	for _, f := range a.Vehicle.FuncQ {
+		if f == nil {
+			continue
+		}
+		f()
+	}
+	a.Vehicle.FuncQ = make([]func(), 5) // Clear the queue.
+
 }
 
 // Func is the integration function using Gaussian VOP as per Ruggiero et al. 2011.
@@ -241,9 +245,7 @@ func (a *Mission) Func(t float64, f []float64) (fDot []float64) {
 	}
 	for i := 0; i < 7; i++ {
 		if math.IsNaN(fDot[i]) {
-			Rcur, Vcur := a.Orbit.GetRV()
-			Rtmp, Vtmp := tmpOrbit.GetRV()
-			panic(fmt.Errorf("fDot[%d]=NaN @ dt=%s\np=%f\th=%f\tdv=%+v\ntmp:%s\ncur:%s\nR_cur=%+v\tV_cur=%+v\nR_tmp=%+v\tV_tmp=%+v", i, a.CurrentDT, p, h, Δv, tmpOrbit, a.Orbit, Rcur, Vcur, Rtmp, Vtmp))
+			panic(fmt.Errorf("fDot[%d]=NaN @ dt=%s\ntmp:%s\ncur:%s\n%.20f\n", i, a.CurrentDT, tmpOrbit, a.Orbit, a.Orbit.e-1))
 		}
 	}
 	return

src/dynamics/spacecraft.go

@@ -77,11 +77,6 @@ func (sc *Spacecraft) Accelerate(dt time.Time, o *Orbit) (Δv []float64, fuel fl
 			sc.logger.Log("level", "info", "subsys", "astro", "date", dt, "thrust", clType, "reason", ctrl.Reason(), "v(km/s)", norm(o.GetV()))
 			sc.prevCL = &clType
 		}
-		// Check if we're in a parabolic orbit and if so, we're activating the action NOW.
-		if p := o.GetSemiParameter(); p < 0 { // TODO: check if this ever happens anymore.
-			sc.logger.Log("level", "critical", "subsys", "astro", "date", dt, "p", p, "action", wp.Action())
-			reached = true
-		}
 		if reached {
 			sc.logger.Log("level", "notice", "subsys", "astro", "waypoint", wp, "status", "completed", "r(km)", norm(o.GetR()), "v (km/s)", norm(o.GetV()))
 			// Handle waypoint action

src/dynamics/waypoints.go

@@ -38,19 +38,6 @@ type Waypoint interface {
 
 // NewOutwardSpiral defines a new outward spiral from a celestial object.
 func NewOutwardSpiral(body CelestialObject, action *WaypointAction) *ReachDistance {
-	/*
-		Add the following possibilities, where:
-		the delta V vector for me is f_r, f_theta, f_h.
-		+Pinkham's spiral:
-		f_r = (q*sqrt(mu*p)*v_r)/(2*r**2)
-		f_theta = (q*sqrt(mu*p)*v_theta)/(2*r**2)
-		p = p_s * exp*(q0)
-		r = (p*(1+q**1))/(1+exp(q0)*(1+q**2)*k*cos(theta-omega))
-
-		+ Lawden's spiral
-		r = (r_s*(sin alpha)^6)/(1-3*(sin(alpha))^2)
-		theta = theta_0 - 4*alpha-3*cotan(alpha) -- cotan == cot
-	*/
 	return &ReachDistance{body.SOI, action, false}
 }