Add different epsilon based on distance from Sun

examples/thesis-inbound/main.go

@@ -13,7 +13,7 @@ func main() {
 	baseDepart := time.Date(2018, 11, 8, 0, 0, 0, 0, time.UTC)
 	// Est arrival is best according to lambert solver
 	estArrival := time.Date(2019, 03, 23, 0, 0, 0, 0, time.UTC)
-	maxPropDT := estArrival.Add(time.Duration(6*31*24) * time.Hour)
+	maxPropDT := estArrival.Add(time.Duration(24*31*24) * time.Hour)
 
 	// Find target hyperbola
 	hypDepart := estArrival.Add(time.Duration(-6*31*24) * time.Hour)
@@ -22,7 +22,7 @@ func main() {
 	hyp.Propagate()
 
 	// Now let's grab the final hyperbolic orbit as the target.
-	inb := smd.NewMission(InboundSpacecraft("inbSC", *hyp.Orbit), InitialMarsOrbit(), baseDepart, maxPropDT, smd.Cartesian, smd.Perturbations{}, smd.ExportConfig{AsCSV: false, Cosmo: true, Filename: "inb"})
+	inb := smd.NewPreciseMission(InboundSpacecraft("inbSC", *hyp.Orbit), InitialMarsOrbit(), baseDepart, maxPropDT, smd.Cartesian, smd.Perturbations{}, time.Minute, smd.ExportConfig{AsCSV: false, Cosmo: true, Filename: "inb"})
 	inb.Propagate()
 }
 

mission.go

@@ -46,13 +46,19 @@ type Mission struct {
 	CurrentDT      time.Time
 	propMethod     Propagator
 	perts          Perturbations
+	step           time.Duration // time step
 	stopChan       chan (bool)
 	histChan       chan<- (MissionState)
 	done, collided bool
 }
 
-// NewMission returns a new Mission instance from the position and velocity vectors.
+// NewMission is the same as NewPreciseMission with the default step size.
 func NewMission(s *Spacecraft, o *Orbit, start, end time.Time, method Propagator, perts Perturbations, conf ExportConfig) *Mission {
+	return NewPreciseMission(s, o, start, end, method, perts, StepSize, conf)
+}
+
+// NewPreciseMission returns a new Mission instance with custom provided time step.
+func NewPreciseMission(s *Spacecraft, o *Orbit, start, end time.Time, method Propagator, perts Perturbations, step time.Duration, conf ExportConfig) *Mission {
 	// If no filepath is provided, then no output will be written.
 	var histChan chan (MissionState)
 	if !conf.IsUseless() {
@@ -73,7 +79,7 @@ func NewMission(s *Spacecraft, o *Orbit, start, end time.Time, method Propagator
 		end = end.UTC()
 	}
 
-	a := &Mission{s, o, start, end, start, method, perts, make(chan (bool), 1), histChan, false, false}
+	a := &Mission{s, o, start, end, start, method, perts, step, make(chan (bool), 1), histChan, false, false}
 	// Write the first data point.
 	if histChan != nil {
 		histChan <- MissionState{a.CurrentDT, *s, *o}
@@ -105,7 +111,7 @@ func (a *Mission) Propagate() {
 		}
 	}()
 	vInit := norm(a.Orbit.V())
-	ode.NewRK4(0, StepSize.Seconds(), a).Solve() // Blocking.
+	ode.NewRK4(0, a.step.Seconds(), a).Solve() // Blocking.
 	vFinal := norm(a.Orbit.V())
 	a.done = true
 	duration := a.CurrentDT.Sub(a.StartDT)
@@ -135,7 +141,7 @@ func (a *Mission) Stop(t float64) bool {
 		}
 		return true // Stop because there is a request to stop.
 	default:
-		a.CurrentDT = a.CurrentDT.Add(StepSize)
+		a.CurrentDT = a.CurrentDT.Add(a.step)
 		if a.EndDT.Before(a.StartDT) {
 			// Check if any waypoint still needs to be reached.
 			for _, wp := range a.Vehicle.WayPoints {

orbit.go

@@ -10,10 +10,15 @@ import (
 )
 
 const (
+	// Precise ε
 	eccentricityε = 5e-5                         // 0.00005
 	angleε        = (5e-3 / 360) * (2 * math.Pi) // 0.005 degrees
 	distanceε     = 3e1                          // 20 km
-	velocityε     = 1e-6                         // in km/s
+	// Coarse ε (for interplanetary flight)
+	eccentricityLgε = 1e-2                         // 0.01
+	angleLgε        = (5e-1 / 360) * (2 * math.Pi) // 0.5 degrees
+	distanceLgε     = 5e2                          // 500 km
+	velocityε       = 1e-6                         // in km/s
 )
 
 // Orbit defines an orbit via its orbital elements.
@@ -212,6 +217,14 @@ func (o Orbit) String() string {
 	return fmt.Sprintf("a=%.1f e=%.4f i=%.3f Ω=%.3f ω=%.3f ν=%.3f λ=%.3f u=%.3f", o.a, o.e, Rad2deg(o.i), Rad2deg(o.Ω), Rad2deg(o.ω), Rad2deg(o.ν), Rad2deg(o.TrueLongλ()), Rad2deg(o.ArgLatitudeU()))
 }
 
+// epsilons returns the epsilons used to determine equality.
+func (o Orbit) epsilons() (float64, float64, float64) {
+	if o.RNorm() > 0.5*AU {
+		return distanceLgε, eccentricityLgε, angleLgε
+	}
+	return distanceε, eccentricityε, angleε
+}
+
 // Equals returns whether two orbits are identical with free true anomaly.
 // Use StrictlyEquals to also check true anomaly.
 func (o Orbit) Equals(o1 Orbit) (bool, error) {

prop.go

@@ -379,6 +379,7 @@ func (cl *OptimalΔOrbit) Control(o Orbit) []float64 {
 		// Note that, as described in Hatten MSc. thesis, the summing method only
 		// works one way (because of the δO^2) per OE. So I added the sign function
 		// *every here and there* as needed that to fix it.
+		dε, eε, aε := o.epsilons()
 		for _, ctrl := range cl.controls {
 			var weight, δO float64
 			p := o.SemiParameter()
@@ -387,19 +388,19 @@ func (cl *OptimalΔOrbit) Control(o Orbit) []float64 {
 			switch ctrl.Type() {
 			case OptiΔaCL:
 				δO = o.a - cl.oTgt.a
-				if math.Abs(δO) < distanceε {
+				if math.Abs(δO) < dε {
 					δO = 0
 				}
 				weight = sign(-δO) * math.Pow(h, 2) / (4 * math.Pow(o.a, 4) * math.Pow(1+o.e, 2))
 			case OptiΔeCL:
 				δO = o.e - cl.oTgt.e
-				if math.Abs(δO) < eccentricityε {
+				if math.Abs(δO) < eε {
 					δO = 0
 				}
 				weight = sign(-δO) * math.Pow(h, 2) / (4 * math.Pow(p, 2))
 			case OptiΔiCL:
 				δO = o.i - cl.oTgt.i
-				if math.Abs(δO) < angleε {
+				if math.Abs(δO) < aε {
 					δO = 0
 				}
 				weight = sign(-δO) * math.Pow((h+o.e*h*math.Cos(o.ω+math.Asin(o.e*sinω)))/(p*(math.Pow(o.e*sinω, 2)-1)), 2)
@@ -409,7 +410,7 @@ func (cl *OptimalΔOrbit) Control(o Orbit) []float64 {
 					// Enforce short path to correct angle.
 					δO *= -1
 				}
-				if math.Abs(δO) < angleε {
+				if math.Abs(δO) < aε {
 					δO = 0
 				}
 				weight = sign(-δO) * math.Pow((h*math.Sin(o.i)*(o.e*math.Sin(o.ω+math.Asin(o.e*cosω))-1))/(p*(1-math.Pow(o.e*cosω, 2))), 2)
@@ -419,7 +420,7 @@ func (cl *OptimalΔOrbit) Control(o Orbit) []float64 {
 					// Enforce short path to correct angle.
 					δO *= -1
 				}
-				if math.Abs(δO) < angleε {
+				if math.Abs(δO) < aε {
 					δO = 0
 				}
 				weight = sign(-δO) * (math.Pow(o.e*h, 2) / (4 * math.Pow(p, 2))) * (1 - math.Pow(o.e, 2)/4)