Export thrusting angles

export.go

@@ -184,8 +184,8 @@ func createInterpolatedFile(filename string, stamped bool, stateDT time.Time) *o
 	return f
 }
 
-// createAsCSVCSVFile returns a file which requires a defer close statement!
-func createAsCSVCSVFile(filename string, conf ExportConfig, stateDT time.Time) *os.File {
+// createAsCSVFile returns a file which requires a defer close statement!
+func createAsCSVFile(filename string, conf ExportConfig, stateDT time.Time) *os.File {
 	if conf.Timestamp {
 		t := time.Now()
 		filename = fmt.Sprintf("%s/orbital-elements-%s-%d-%02d-%02dT%02d.%02d.%02d.csv", smdConfig().outputDir, filename, t.Year(), t.Month(), t.Day(), t.Hour(), t.Minute(), t.Second())
@@ -258,7 +258,7 @@ func StreamStates(conf ExportConfig, stateChan <-chan (State)) {
 					}
 				}
 				if conf.AsCSV {
-					fAsCSV = createAsCSVCSVFile(fmt.Sprintf("%s-%d", conf.Filename, fileNo), conf, state.DT)
+					fAsCSV = createAsCSVFile(fmt.Sprintf("%s-%d", conf.Filename, fileNo), conf, state.DT)
 				}
 				fileNo++
 			} else {
@@ -283,7 +283,7 @@ func StreamStates(conf ExportConfig, stateChan <-chan (State)) {
 					}
 					if conf.AsCSV {
 						fAsCSV.WriteString(fmt.Sprintf("\n# Simulation time end (UTC): %s\n", state.DT.UTC()))
-						fAsCSV = createAsCSVCSVFile(fmt.Sprintf("%s-%d", conf.Filename, fileNo), conf, state.DT)
+						fAsCSV = createAsCSVFile(fmt.Sprintf("%s-%d", conf.Filename, fileNo), conf, state.DT)
 					}
 					fileNo++
 					// Force writing this data point now instead of creating N new files.
@@ -380,3 +380,55 @@ type ExportConfig struct {
 func (c ExportConfig) IsUseless() bool {
 	return !c.Cosmo && !c.AsCSV
 }
+
+// ThurstAngleExport configures the exporting of the simulation. Exports in CSV
+// in the format of "datetime,alpha,beta"
+type ThurstAngleExport struct {
+	enabled   bool
+	timestamp bool
+	latestDT  time.Time
+	fhdl      *os.File
+}
+
+func (t *ThurstAngleExport) createFile(filename string) {
+	if !t.enabled {
+		return
+	}
+	osfn := smdConfig().outputDir
+	if osfn == "" {
+		fmt.Println("\nWARNING: SMD_CONFIG outputDir undefined. Forced to `./`")
+		osfn = "./"
+	}
+
+	if t.timestamp {
+		dt := time.Now()
+		osfn += fmt.Sprintf("/thrusting-angles-%s-%d-%02d-%02dT%02d.%02d.%02d.csv", filename, dt.Year(), dt.Month(), dt.Day(), dt.Hour(), dt.Minute(), dt.Second())
+	} else {
+		osfn += fmt.Sprintf("/thrusting-angles-%s.csv", filename)
+	}
+	f, err := os.Create(osfn)
+	if err != nil {
+		panic(err)
+	}
+	// Header
+	f.WriteString(fmt.Sprintf(`# Creation date (UTC): %s
+# Records are datetime, angle, beta. All angles are in degrees.
+datetime,alpha,beta,`, time.Now()))
+	f.WriteString("\n")
+	t.fhdl = f
+}
+
+// Store appends a new tuple of values for the given time. Note that if the previous time is the same as the new one, the new entry is *ignored*.
+func (t *ThurstAngleExport) Store(dt time.Time, alpha, beta float64) {
+	if t.latestDT == dt {
+		return
+	}
+	t.fhdl.WriteString(fmt.Sprintf("\"%s\",%3.6f,%3.6f,\n", dt, alpha, beta))
+	t.latestDT = dt
+}
+
+func NewThurstAngleExport(filename string, timestamped bool) *ThurstAngleExport {
+	txp := ThurstAngleExport{true, timestamped, time.Date(1900, 1, 1, 0, 0, 0, 0, time.UTC), nil}
+	txp.createFile(filename)
+	return &txp
+}

mission_test.go

@@ -201,9 +201,10 @@ func TestCorrectOEa(t *testing.T) {
 		EPThrusters := []EPThruster{new(PPS1350)}
 		dryMass := 300.0
 		fuelMass := 67.0
-		otgt := NewOrbitTarget(*oTarget, nil, meth, OptiΔaCL)
-		otgt.SetEpsilons(1e-1, 1e-5, Deg2rad(0.001))
-		sc := NewSpacecraft("COE", dryMass, fuelMass, eps, EPThrusters, false, []*Cargo{}, []Waypoint{otgt})
+		wpOrbitTgt := NewOrbitTarget(*oTarget, nil, meth, OptiΔaCL)
+		wpOrbitTgt.SetEpsilons(1e-1, 1e-5, Deg2rad(0.001))
+		wpOrbitTgt.SetExport(NewThurstAngleExport("correctOE-sma", true))
+		sc := NewSpacecraft("COE", dryMass, fuelMass, eps, EPThrusters, false, []*Cargo{}, []Waypoint{wpOrbitTgt})
 		start := time.Date(2017, 1, 1, 0, 0, 0, 0, time.UTC)
 		end := start.Add(time.Duration(45*24) * time.Hour)
 		astro := NewMission(sc, oInit, start, end, Perturbations{}, false, ExportConfig{Filename: fmt.Sprintf("ruggOEa-%s", meth), Cosmo: smdConfig().testExport, AsCSV: smdConfig().testExport})
@@ -258,7 +259,9 @@ func TestCorrectOEi(t *testing.T) {
 		EPThrusters := []EPThruster{new(PPS1350)}
 		dryMass := 300.0
 		fuelMass := 67.0
-		sc := NewSpacecraft("COE", dryMass, fuelMass, eps, EPThrusters, false, []*Cargo{}, []Waypoint{NewOrbitTarget(*oTarget, nil, meth, OptiΔiCL)})
+		wpOrbitTgt := NewOrbitTarget(*oTarget, nil, meth, OptiΔiCL)
+		wpOrbitTgt.SetExport(NewThurstAngleExport("correctOE-inc", false))
+		sc := NewSpacecraft("COE", dryMass, fuelMass, eps, EPThrusters, false, []*Cargo{}, []Waypoint{wpOrbitTgt})
 		start := time.Date(2017, 1, 1, 0, 0, 0, 0, time.UTC)
 		end := start.Add(time.Duration(55*24) * time.Hour)
 		astro := NewMission(sc, oInit, start, end, Perturbations{}, false, ExportConfig{})

prop.go

@@ -178,6 +178,13 @@ func unitΔvFromAngles(α, β float64) []float64 {
 	return []float64{sinα * cosβ, cosα * cosβ, sinβ}
 }
 
+func anglesFromUnitΔv(Δv []float64) (α, β float64) {
+	β = math.Asin(Δv[2])
+	cosβ := math.Cos(β)
+	α = math.Asin(Δv[0] / cosβ)
+	return
+}
+
 // OptimalThrust is an optimal thrust.
 type OptimalThrust struct {
 	ctrl func(o Orbit) []float64

waypoints.go

@@ -143,6 +143,7 @@ type OrbitTarget struct {
 	target  Orbit
 	ctrl    *OptimalΔOrbit
 	action  *WaypointAction
+	xprt    *ThurstAngleExport
 	cleared bool
 }
 
@@ -177,15 +178,24 @@ func (wp *OrbitTarget) ThrustDirection(o Orbit, dt time.Time) (ThrustControl, bo
 		fmt.Printf("[WARNING] OrbitTarget reached @%s *but* %s: %s\n", dt, err, o.String())
 		wp.cleared = true
 	}
+	if wp.xprt != nil {
+		α, β := anglesFromUnitΔv(Unit(wp.ctrl.Control(o)))
+		wp.xprt.Store(dt, Rad2deg(α), Rad2deg(β))
+	}
 	return wp.ctrl, wp.cleared
 }
 
+// ThrustDirection implements the optimal orbit target.
+func (wp *OrbitTarget) SetExport(t *ThurstAngleExport) {
+	wp.xprt = t
+}
+
 // NewOrbitTarget defines a new orbit target.
 func NewOrbitTarget(target Orbit, action *WaypointAction, meth ControlLawType, laws ...ControlLaw) *OrbitTarget {
 	if target.Periapsis() < target.Origin.Radius || target.Apoapsis() < target.Origin.Radius {
 		fmt.Printf("[WARNING] Target orbit on collision course with %s\n", target.Origin)
 	}
-	return &OrbitTarget{target, NewOptimalΔOrbit(target, meth, laws), action, false}
+	return &OrbitTarget{target, NewOptimalΔOrbit(target, meth, laws), action, nil, false}
 }
 
 // HohmannTransfer allows to perform an Hohmann transfer.