From e5204472cedec6a21f18fdac1ee78d3be72eb69f Mon Sep 17 00:00:00 2001
From: mclrc <moritzamando@proton.me>
Date: Mon, 16 Oct 2023 15:59:28 +0300
Subject: [PATCH] Add Euler's method

---
 src/cli.rs         |  2 +-
 src/propagate.rs   |  2 ++
 src/run.rs         |  1 +
 src/solvers.rs     | 51 +++++++++++++++++++++++++++++++++++++++++++++-
 src/spice_utils.rs |  2 +-
 src/tests.rs       | 19 ++++++++++++++++-
 6 files changed, 73 insertions(+), 4 deletions(-)

diff --git a/src/cli.rs b/src/cli.rs
index 7818704..9a5f105 100644
--- a/src/cli.rs
+++ b/src/cli.rs
@@ -77,6 +77,6 @@ pub struct Args {
 	)]
 	pub cb_id: Option<i32>,
 
-	#[clap(long, value_name = "rk4|dopri45", help = "Integration method")]
+	#[clap(long, value_name = "rk4|dopri45|euler", help = "Integration method")]
 	pub method: Option<String>,
 }
diff --git a/src/propagate.rs b/src/propagate.rs
index 6bbc283..fe50c56 100644
--- a/src/propagate.rs
+++ b/src/propagate.rs
@@ -5,6 +5,7 @@ use ndarray::Array1;
 
 pub enum SolverConfig {
 	Rk4 { h: f64 },
+	Euler { h: f64 },
 	Dopri45 { h: f64, atol: f64, rtol: f64 },
 }
 
@@ -61,6 +62,7 @@ pub fn propagate(
 	// Create solver object based on config on the heap (since exact type is unknown)
 	let mut solver: Box<dyn solvers::Solver> = match solver {
 		SolverConfig::Rk4 { h } => Box::new(solvers::Rk4::new(f, h, et0, &y0, etfinal)),
+		SolverConfig::Euler { h } => Box::new(solvers::Euler::new(f, h, et0, &y0, etfinal)),
 		SolverConfig::Dopri45 { h, atol, rtol } => {
 			Box::new(solvers::Dopri45::new(f, h, et0, &y0, etfinal, atol, rtol))
 		}
diff --git a/src/run.rs b/src/run.rs
index 290b613..1226ea3 100644
--- a/src/run.rs
+++ b/src/run.rs
@@ -41,6 +41,7 @@ pub fn run(
 	// Create solver config based on CLI args
 	let solver = match method.as_deref() {
 		Some("rk4") | None => propagate::SolverConfig::Rk4 { h },
+		Some("euler") => propagate::SolverConfig::Euler { h },
 		Some("dopri45") => propagate::SolverConfig::Dopri45 {
 			h,
 			atol: atol.unwrap_or(50000f64),
diff --git a/src/solvers.rs b/src/solvers.rs
index 1278443..d5977ee 100644
--- a/src/solvers.rs
+++ b/src/solvers.rs
@@ -8,7 +8,7 @@ pub mod step_fns {
 		x: f64,
 		y: &Array1<f64>,
 		h: f64,
-		) -> Result<(f64, Array1<f64>), String> {
+	) -> Result<(f64, Array1<f64>), String> {
 		let k1 = f(x, y)?;
 		let k2 = f(x + 0.5 * h, &(y + h * &k1*0.5))?;
 		let k3 = f(x + 0.5 * h, &(y + h * &k2*0.5))?;
@@ -18,6 +18,16 @@ pub mod step_fns {
 
 		Ok((x + h, next_y))
 	}
+	pub fn euler(
+		f: impl Fn(f64, &Array1<f64>) -> Result<Array1<f64>, String>,
+		x: f64,
+		y: &Array1<f64>,
+		h: f64,
+	) -> Result<(f64, Array1<f64>), String> {
+		let dy = f(x, y)?;
+		let next_y = y + dy * h;
+		Ok((x + h, next_y))
+	}
 	#[allow(clippy::too_many_arguments)]
 	pub fn dopri(
 		f: impl Fn(f64, &Array1<f64>) -> Result<Array1<f64>, String>,
@@ -95,6 +105,45 @@ where
 	}
 }
 
+pub struct Euler<F> {
+	f: F,
+	x: f64,
+	y: Array1<f64>,
+	h: f64,
+	xmax: f64,
+}
+
+impl<F> Euler<F> {
+	pub fn new(f: F, h: f64, x0: f64, y0: &Array1<f64>, xmax: f64) -> Self {
+		Self {
+			f,
+			x: x0,
+			y: y0.clone(),
+			h,
+			xmax,
+		}
+	}
+}
+
+impl<F> Solver for Euler<F>
+where
+	F: Fn(f64, &Array1<f64>) -> Result<Array1<f64>, String>,
+{
+	fn next_state(&mut self) -> Result<Option<(f64, Array1<f64>)>, String> {
+		if self.x >= self.xmax {
+			return Ok(None);
+		} else if self.x + self.h > self.xmax {
+			self.h = self.xmax - self.x;
+		}
+
+		let p = step_fns::euler(|x, y| (self.f)(x, y), self.x, &self.y, self.h)?;
+		self.x = p.0;
+		self.y = p.1.clone();
+
+		Ok(Some(p.clone()))
+	}
+}
+
 pub struct Dopri45<F> {
 	f: F,
 	h: f64,
diff --git a/src/spice_utils.rs b/src/spice_utils.rs
index 6fd2f39..4df540f 100644
--- a/src/spice_utils.rs
+++ b/src/spice_utils.rs
@@ -8,7 +8,7 @@ pub fn set_error_handling(action: &str, len: &str, dev: &str) {
 	unsafe {
 		spice::c::errprt_c(spice::cstr!("set"), 0, spice::cstr!(len));
 		spice::c::erract_c(spice::cstr!("set"), 20, spice::cstr!(action));
-		spice::c::errdev_c(spice::cstr!("SET"), 0, spice::cstr!(dev));
+		spice::c::errdev_c(spice::cstr!("set"), 0, spice::cstr!(dev));
 	}
 }
 
diff --git a/src/tests.rs b/src/tests.rs
index b3245f3..b858670 100644
--- a/src/tests.rs
+++ b/src/tests.rs
@@ -49,6 +49,23 @@ fn run_test_scenario(
 	println!();
 }
 
+#[test]
+#[serial]
+fn maven_cruise_euler() {
+	println!("starting euler test");
+	run_test_scenario(
+		"spice/maven_cruise.bsp",
+		"2013-NOV-20",
+		Some(&["Sun", "Earth", "Jupiter Barycenter", "Mars"]),
+		Some(&["Maven"]),
+		None,
+		"2014-SEP-21",
+		100f64,
+		"euler",
+		None,
+	)
+}
+
 #[test]
 #[serial]
 fn maven_cruise_rk4() {
@@ -83,7 +100,7 @@ fn voyager2_flyby_dopri45() {
 
 #[test]
 #[serial]
-fn spkattractors() {
+fn spk_attractors() {
 	run_test_scenario(
 		"spice/voyager2_flyby.bsp",
 		"1978-JAN-23",