From 33e75242617e8d1f214461fe6304019533123dc4 Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Thu, 30 Nov 2023 18:30:02 +0100
Subject: [PATCH 01/16] sv relativistic effects compensation

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/candidate.rs | 18 +++---------------
 src/cfg.rs       | 10 +++++-----
 src/solver.rs    | 44 ++++++++++++++++++++++++++++++++------------
 3 files changed, 40 insertions(+), 32 deletions(-)

diff --git a/src/candidate.rs b/src/candidate.rs
index 745be6d..a8cacfd 100644
--- a/src/candidate.rs
+++ b/src/candidate.rs
@@ -3,6 +3,7 @@
 use gnss::prelude::SV;
 use hifitime::Unit;
 use log::debug;
+use map_3d::deg2rad;
 use nyx::cosmic::SPEED_OF_LIGHT;
 use nyx::linalg::{DVector, MatrixXx4};
 
@@ -13,6 +14,7 @@ use crate::{
     bias::{IonosphericBias, TropoModel, TroposphericBias},
     Error,
 };
+use nyx::md::prelude::Frame;
 
 /// Signal observation to attach to each candidate
 #[derive(Debug, Default, Clone)]
@@ -153,26 +155,12 @@ impl Candidate {
         })
     }
      */
-    /*
-     * Returns IONOD possibly impacting
-    pub(crate) fn ionod_model(&self, frequency: f64) -> Option<Duration> {
-        self.modeled_ionod
-            .iter()
-            .filter_map(|ionod| {
-                if ionod.frequency == frequency {
-                    Some(ionod.time_delay)
-                } else {
-                    None
-                }
-            })
-            .reduce(|k, _| k)
-    }
-     */
     /*
      * Computes signal transmission Epoch
      * returns (t_tx, dt_ttx)
      * "t_tx": Epoch in given timescale
      * "dt_ttx": elapsed duration in seconds in given timescale
+     * "frame": Solid body reference Frame
      */
     pub(crate) fn transmission_time(&self, cfg: &Config) -> Result<(Epoch, f64), Error> {
         let (t, ts) = (self.t, self.t.time_scale);
diff --git a/src/cfg.rs b/src/cfg.rs
index 8530650..688d559 100644
--- a/src/cfg.rs
+++ b/src/cfg.rs
@@ -75,8 +75,8 @@ fn default_earth_rot() -> bool {
     true
 }
 
-fn default_rel_clock_corr() -> bool {
-    false
+fn default_relativistic_effect() -> bool {
+    true
 }
 
 fn default_lsq_weight() -> Option<LSQWeight> {
@@ -106,6 +106,8 @@ pub struct Modeling {
     #[cfg_attr(feature = "serde", serde(default))]
     pub sv_clock_bias: bool,
     #[cfg_attr(feature = "serde", serde(default))]
+    pub sv_relativistic_effect: bool,
+    #[cfg_attr(feature = "serde", serde(default))]
     pub tropo_delay: bool,
     #[cfg_attr(feature = "serde", serde(default))]
     pub iono_delay: bool,
@@ -113,8 +115,6 @@ pub struct Modeling {
     pub sv_total_group_delay: bool,
     #[cfg_attr(feature = "serde", serde(default))]
     pub earth_rotation: bool,
-    #[cfg_attr(feature = "serde", serde(default))]
-    pub relativistic_clock_corr: bool,
 }
 
 impl Default for Modeling {
@@ -125,7 +125,7 @@ impl Default for Modeling {
             tropo_delay: default_tropo(),
             sv_total_group_delay: default_sv_tgd(),
             earth_rotation: default_earth_rot(),
-            relativistic_clock_corr: default_rel_clock_corr(),
+            sv_relativistic_effect: default_relativistic_effect(),
         }
     }
 }
diff --git a/src/solver.rs b/src/solver.rs
index 19dea88..5005527 100644
--- a/src/solver.rs
+++ b/src/solver.rs
@@ -1,12 +1,14 @@
 //! PVT solver
 use std::collections::HashMap;
 
-use hifitime::Epoch;
+use hifitime::{Epoch, Unit};
 use log::{debug, error, warn};
+use map_3d::deg2rad;
 use thiserror::Error;
 
 use nyx::cosmic::eclipse::{eclipse_state, EclipseState};
 use nyx::cosmic::Orbit;
+use nyx::cosmic::SPEED_OF_LIGHT;
 use nyx::md::prelude::{Arc, Cosm};
 use nyx::md::prelude::{Bodies, Frame, LightTimeCalc};
 
@@ -116,7 +118,16 @@ impl InterpolationResult {
     pub(crate) fn orbit(&self, dt: Epoch, frame: Frame) -> Orbit {
         let (x, y, z) = self.position();
         let (v_x, v_y, v_z) = self.velocity().unwrap_or((0.0_f64, 0.0_f64, 0.0_f64));
-        Orbit::cartesian(x, y, z, v_x / 1000.0, v_y / 1000.0, v_z / 1000.0, dt, frame)
+        Orbit::cartesian(
+            x / 1000.0,
+            y / 1000.0,
+            z / 1000.0,
+            v_x / 1000.0,
+            v_y / 1000.0,
+            v_z / 1000.0,
+            dt,
+            frame,
+        )
     }
 }
 
@@ -159,14 +170,9 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         let cosmic = Cosm::de438();
         let sun_frame = cosmic.frame("Sun J2000");
         let earth_frame = cosmic.frame("EME2000");
-
         /*
          * print some infos on latched config
          */
-        if cfg.modeling.relativistic_clock_corr {
-            warn!("relativistic clock corr. is not feasible at the moment");
-        }
-
         if mode.is_spp() && cfg.min_sv_sunlight_rate.is_some() {
             warn!("eclipse filter is not meaningful when using spp strategy");
         }
@@ -221,23 +227,37 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         /* interpolate positions */
         let mut pool: Vec<Candidate> = pool
             .iter()
-            .filter_map(|c| {
+            .filter_map(|mut c| {
                 let (t_tx, dt_ttx) = c.transmission_time(&self.cfg).ok()?;
 
                 if let Some(mut interpolated) = (self.interpolator)(t_tx, c.sv, interp_order) {
                     let mut c = c.clone();
-
+                    if modeling.sv_relativistic_effect {
+                        if interpolated.velocity.is_some() { // otherwise, following calaculations would diverge
+                            let orbit = interpolated.orbit(t_tx, self.earth_frame);
+                            const EARTH_SEMI_MAJOR_AXIS_WGS84: f64 = 6378137.0_f64;
+                            const EARTH_GRAVITATIONAL_CONST: f64 = 3986004.418 * 10.0E8;
+                            let orbit = interpolated.orbit(t_tx, self.earth_frame);
+                            let ea_rad = deg2rad(orbit.ea_deg());
+                            let gm = (EARTH_SEMI_MAJOR_AXIS_WGS84 * EARTH_GRAVITATIONAL_CONST).sqrt();
+                            let bias = -2.0_f64 * orbit.ecc() * ea_rad.sin() * gm / SPEED_OF_LIGHT / SPEED_OF_LIGHT * Unit::Second;
+                            debug!(
+                                "{:?} ({}) : relativistic effect: {}",
+                                t_tx, c.sv, bias
+                            );
+                            c.clock_corr += bias;
+                        } else {
+                            warn!("{:?} ({}) : relativistic effect can't be compensated without velocity interpolation", t_tx, c.sv);
+                        }
+                    }
                     if modeling.earth_rotation {
                         const EARTH_OMEGA_E_WGS84: f64 = 7.2921151467E-5;
-
                         let we = EARTH_OMEGA_E_WGS84 * dt_ttx;
                         let (we_cos, we_sin) = (we.cos(), we.sin());
-
                         let r = Matrix3::<f64>::new(
                             we_cos, we_sin, 0.0_f64, -we_sin, we_cos, 0.0_f64, 0.0_f64, 0.0_f64,
                             1.0_f64,
                         );
-
                         interpolated.position = r * interpolated.position;
                     }
 

From bb2184bdd9222b529146cad0ad96fbc817e3cdc5 Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Thu, 30 Nov 2023 22:00:53 +0100
Subject: [PATCH 02/16] relativistic effect compensation

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/candidate.rs |  5 +++-
 src/cfg.rs       | 25 +++++++++--------
 src/solver.rs    | 73 ++++++++++++++++++++++++++++++++++--------------
 3 files changed, 70 insertions(+), 33 deletions(-)

diff --git a/src/candidate.rs b/src/candidate.rs
index a8cacfd..e22c506 100644
--- a/src/candidate.rs
+++ b/src/candidate.rs
@@ -229,7 +229,10 @@ impl Candidate {
         g[(row_index, 2)] = (z0 - sv_z) / rho;
         g[(row_index, 3)] = 1.0_f64;
 
-        let mut models = -clock_corr * SPEED_OF_LIGHT;
+        let mut models = 0.0_f64;
+        if cfg.modeling.sv_clock_bias {
+            models -= clock_corr * SPEED_OF_LIGHT;
+        }
 
         /*
          * Possible delay compensations
diff --git a/src/cfg.rs b/src/cfg.rs
index 688d559..b94a3ba 100644
--- a/src/cfg.rs
+++ b/src/cfg.rs
@@ -75,10 +75,14 @@ fn default_earth_rot() -> bool {
     true
 }
 
-fn default_relativistic_effect() -> bool {
+fn default_relativistic_clock_bias() -> bool {
     true
 }
 
+fn default_relativistic_path_range() -> bool {
+    false
+}
+
 fn default_lsq_weight() -> Option<LSQWeight> {
     None
 }
@@ -106,14 +110,16 @@ pub struct Modeling {
     #[cfg_attr(feature = "serde", serde(default))]
     pub sv_clock_bias: bool,
     #[cfg_attr(feature = "serde", serde(default))]
-    pub sv_relativistic_effect: bool,
+    pub sv_total_group_delay: bool,
+    #[cfg_attr(feature = "serde", serde(default))]
+    pub relativistic_clock_bias: bool,
+    #[cfg_attr(feature = "serde", serde(default))]
+    pub relativistic_path_range: bool,
     #[cfg_attr(feature = "serde", serde(default))]
     pub tropo_delay: bool,
     #[cfg_attr(feature = "serde", serde(default))]
     pub iono_delay: bool,
     #[cfg_attr(feature = "serde", serde(default))]
-    pub sv_total_group_delay: bool,
-    #[cfg_attr(feature = "serde", serde(default))]
     pub earth_rotation: bool,
 }
 
@@ -125,7 +131,8 @@ impl Default for Modeling {
             tropo_delay: default_tropo(),
             sv_total_group_delay: default_sv_tgd(),
             earth_rotation: default_earth_rot(),
-            sv_relativistic_effect: default_relativistic_effect(),
+            relativistic_clock_bias: default_relativistic_clock_bias(),
+            relativistic_path_range: default_relativistic_path_range(),
         }
     }
 }
@@ -135,11 +142,7 @@ impl From<Mode> for Modeling {
         let s = Self::default();
         match mode {
             Mode::PPP => {
-                // TODO: unlock this please
-                // s.earth_rotation = true;
-                // TODO : unlock this please
-                // s.relativistic_clock_corr = true;
-                // TODO:
+                // TODO ?
             },
             _ => {},
         }
@@ -208,7 +211,7 @@ impl Config {
                 interp_order: default_interp(),
                 code_smoothing: default_smoothing(),
                 min_sv_sunlight_rate: None,
-                min_sv_elev: Some(10.0),
+                min_sv_elev: Some(15.0),
                 min_snr: Some(30.0),
                 modeling: Modeling::default(),
                 max_sv: default_max_sv(),
diff --git a/src/solver.rs b/src/solver.rs
index 5005527..8ca441e 100644
--- a/src/solver.rs
+++ b/src/solver.rs
@@ -154,6 +154,8 @@ where
     earth_frame: Frame,
     /* Sun frame */
     sun_frame: Frame,
+    /* prev. state vector for speed determination */
+    prev_sv_state: HashMap<SV, (Epoch, Vector3<f64>)>,
     /* prev. solution for internal logic */
     prev_pvt: Option<(Epoch, PVTSolution)>,
     /* prev. estimate for recursive impl. */
@@ -170,12 +172,16 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         let cosmic = Cosm::de438();
         let sun_frame = cosmic.frame("Sun J2000");
         let earth_frame = cosmic.frame("EME2000");
+
         /*
          * print some infos on latched config
          */
         if mode.is_spp() && cfg.min_sv_sunlight_rate.is_some() {
             warn!("eclipse filter is not meaningful when using spp strategy");
         }
+        if cfg.modeling.relativistic_path_range {
+            warn!("relativistic path range cannot be modeled at the moment");
+        }
 
         Ok(Self {
             mode,
@@ -185,8 +191,9 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
             apriori,
             interpolator,
             cfg: cfg.clone(),
-            prev_pvt: Option::<(Epoch, PVTSolution)>::None,
+            prev_sv_state: HashMap::new(),
             prev_state: Option::<Estimate>::None,
+            prev_pvt: Option::<(Epoch, PVTSolution)>::None,
         })
     }
     /// Try to resolve a PVTSolution at desired "t".
@@ -232,33 +239,57 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
 
                 if let Some(mut interpolated) = (self.interpolator)(t_tx, c.sv, interp_order) {
                     let mut c = c.clone();
-                    if modeling.sv_relativistic_effect {
-                        if interpolated.velocity.is_some() { // otherwise, following calaculations would diverge
+
+                    let rot = match modeling.earth_rotation {
+                        true => {
+                            const EARTH_OMEGA_E_WGS84: f64 = 7.2921151467E-5;
+                            let we = EARTH_OMEGA_E_WGS84 * dt_ttx;
+                            let (we_cos, we_sin) = (we.cos(), we.sin());
+                            Matrix3::<f64>::new(
+                                we_cos, we_sin, 0.0_f64, -we_sin, we_cos, 0.0_f64, 0.0_f64,
+                                0.0_f64, 1.0_f64,
+                            )
+                        },
+                        false => Matrix3::<f64>::new(
+                            1.0_f64, 0.0_f64, 0.0_f64, 0.0_f64, 1.0_f64, 0.0_f64, 0.0_f64, 0.0_f64,
+                            1.0_f64,
+                        ),
+                    };
+
+                    interpolated.position = rot * interpolated.position;
+
+                    if modeling.relativistic_clock_bias {
+                        if interpolated.velocity.is_none() {
+                            /*
+                             * we're interested in determining inst. speed vector
+                             * but the user did not provide it, let's eval. it ourselves
+                             */
+                            if let Some((p_ttx, p_position)) = self.prev_sv_state.get(&c.sv) {
+                                let dt = (t_tx - *p_ttx).to_seconds();
+                                interpolated.velocity =
+                                    Some((rot * interpolated.position - rot * p_position) / dt);
+                            }
+                        }
+
+                        if interpolated.velocity.is_some() {
+                            // otherwise, following calaculations would diverge
                             let orbit = interpolated.orbit(t_tx, self.earth_frame);
                             const EARTH_SEMI_MAJOR_AXIS_WGS84: f64 = 6378137.0_f64;
                             const EARTH_GRAVITATIONAL_CONST: f64 = 3986004.418 * 10.0E8;
                             let orbit = interpolated.orbit(t_tx, self.earth_frame);
                             let ea_rad = deg2rad(orbit.ea_deg());
-                            let gm = (EARTH_SEMI_MAJOR_AXIS_WGS84 * EARTH_GRAVITATIONAL_CONST).sqrt();
-                            let bias = -2.0_f64 * orbit.ecc() * ea_rad.sin() * gm / SPEED_OF_LIGHT / SPEED_OF_LIGHT * Unit::Second;
-                            debug!(
-                                "{:?} ({}) : relativistic effect: {}",
-                                t_tx, c.sv, bias
-                            );
+                            let gm =
+                                (EARTH_SEMI_MAJOR_AXIS_WGS84 * EARTH_GRAVITATIONAL_CONST).sqrt();
+                            let bias = -2.0_f64 * orbit.ecc() * ea_rad.sin() * gm
+                                / SPEED_OF_LIGHT
+                                / SPEED_OF_LIGHT
+                                * Unit::Second;
+                            debug!("{:?} ({}) : relativistic clock bias: {}", t_tx, c.sv, bias);
                             c.clock_corr += bias;
-                        } else {
-                            warn!("{:?} ({}) : relativistic effect can't be compensated without velocity interpolation", t_tx, c.sv);
                         }
-                    }
-                    if modeling.earth_rotation {
-                        const EARTH_OMEGA_E_WGS84: f64 = 7.2921151467E-5;
-                        let we = EARTH_OMEGA_E_WGS84 * dt_ttx;
-                        let (we_cos, we_sin) = (we.cos(), we.sin());
-                        let r = Matrix3::<f64>::new(
-                            we_cos, we_sin, 0.0_f64, -we_sin, we_cos, 0.0_f64, 0.0_f64, 0.0_f64,
-                            1.0_f64,
-                        );
-                        interpolated.position = r * interpolated.position;
+
+                        self.prev_sv_state
+                            .insert(c.sv, (t_tx, interpolated.position));
                     }
 
                     debug!(

From c566324588930d03b2ee759ef9f12f66a94a766c Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Fri, 1 Dec 2023 09:52:14 +0100
Subject: [PATCH 03/16] update doc

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 README.md     | 99 ++++++++++++++++++++++++++++++---------------------
 doc/config.md | 10 ++++++
 doc/lsqspp.md |  4 +++
 doc/ppp.md    |  4 +++
 doc/spp.md    | 32 +++++++++++++++++
 5 files changed, 108 insertions(+), 41 deletions(-)
 create mode 100644 doc/config.md
 create mode 100644 doc/lsqspp.md
 create mode 100644 doc/ppp.md
 create mode 100644 doc/spp.md

diff --git a/README.md b/README.md
index 1111a84..c998152 100644
--- a/README.md
+++ b/README.md
@@ -8,46 +8,76 @@ GNSS-RTK
 
 Precise position solver :crab:
 
-Solving method
-==============
+Notes on this ecosystem
+=======================
 
-The solver uses linear algebra implemented in the `nalgebra` library.
+We use `nalgebra` in the solving process (Matrix, Vectors..).  
+We use `nyx` for ephemerides, orbital calculations, navigation filters, etc..  
+`Hifitime` is at the very basis of this work and allows pretty much everything here.  
+[The RINEX Wiki](https://github.com/georust/rinex/wiki) describes extensive application of this framework.  
+[The RNX2CGGTTS and its Wiki](https://github.com/georust/rinex/wiki) is another interesting application of this framework.
 
-The current solver is straightforward and does not require initialization cycles.  
-We can resolve a PVT for every input as long as the criterias for the current setup are met.
+PVT Solutions
+=============
 
-Some criterias are fixed by physics, others are customized and adjusted in the `Cfg` structure.
+The objective is to resolve a PVT solution, implemented in the form of a `prelude::PVTSolution`,
+ideally the most precise as possible, at a given _Epoch_.
 
-The minimum requirements to form a solution :
+Resolving a PVT requires a minimum number of SV (actual measurements)
+in sight and within the predefined criterias:
 
-- the minimal number of SV within the criterias are in sight: 
-  - 4 SV in default mode
-  - 3 SV in fixed altitude mode
-  - 1 SV in time only mode
-- user was able to provide observations that satisfy the resolution strategy
-- user was able to interpolate the SV state vector at the required _Epoch_
+- 4 SV in default mode
+- 3 SV in fixed altitude mode
+- 1 SV in time only mode
 
-PVT Solutions
-=============
+The predefined criterias are manually set in the configuration file, 
+refer to its [own documentation](./doc/config.md). This means the criterias can be
+loosened or tightened, depending on what you want to achieve.
+
+Some constraints on the measurements may apply too, as a rule of thumb:
+
+- `SPP` strategies will only require Pseudo Range and will resolve without Phase observations.
+This makes these strategies capapble to deploy on constrained and degraged environment.
+For example: old receivers, single constellation in sight.. etc..
+
+- `PPP` strategies require not only Pseudo Range but also Phase observations
+and the sampling of two different radio frequencies. No solutions will be formed
+if this predicate (on top of all the others explained here) do not stand.
+
+We support several methods (also sometimes refered to as _strategies_) which
+will have the solver behave differently, due to their very different nature.   
+The method is set by the `solver::Mode`.
 
-The solver tries to resolve a Position Velocity Time (PVT) solution of the receiver.  
-Currently the Velocity term is not evaluated, therefore we only output Positions and Time errors.
-Dilution of precisions are estimated and attached to each solution.
+Dilution of precisions along other meaningful information are attached to each PVT solution.
 
-Single Point Position (SPP)
-===========================
+This solver will eventually be able to resolve along any known GNSS constellation,
+including a mixed pool of spacecrafts (for example a combination of GAL + GPS), 
+and also express the PVT solution against any supported GNSS Time system (for example GST).
 
-The solver supports the SPP resolution method.
+Strategies and behavior
+=======================
 
-SPP can be deployed in degraded conditions where only one carrier signal and/or a unique constellation is in sight.
-In such conditions:
+We support a couple of strategies, only advanced strategies will give the best results
+(most precise solutions):
 
-- you can only hope for a precision of a few meters.    
-- an interpolation order of 9 makes sense, going above will increase the computation load without any benefit.
-- the solver will have to estimate the total bias due to Ionospheric delay. If you can define
-these components accurately yourself, you are encouranged to do it (see related API section and following paragraphs)
+- [spp](./doc/spp.md)
+- [lsqspp](./doc/lsqspp.md)
+- [ppp](./doc/ppp.md)
 
-## Atmospherical and Environmental conditions modeling
+As a rule of thumb: 
+
+- Non recursive strategies (like [spp](./doc/spp.md)) will generate
+a solution as along as enough signal measurements were provided.
+- Any physical phenomena can be accounted for in this framework,
+on any strategy, even though some will not be meaningful depending on the
+configuration setup.  
+- The GDOP criteria can still be used to reject PVT solutions of non recursive strategies
+- Only recursive strategies will take truly use other solver options of the configuration file.
+
+The current API allows changing the strategy from one Epoch to another and this framework will most likely behave fine.   
+But note that this has not been tested and is not the typical use of such tool.
+
+## Atmospherical and Environmental biases
 
 This solver is always capable of modelizing all conditions and form a solution.   
 It is important to understand how our API is designed and operate it the best you can to get the best results.
@@ -90,16 +120,3 @@ It is important to understand how, when and what to customize depending on your
 
 When working with the "cli" application, you can provide an RTKConfiguration
 in the form of JSON, with `--rtk-cfg`.
-
-RINEX Post Processing
-=====================
-
-The [rinex-cli application](https://github.com/georust/rinex) is there
-to post process RINEX data and resolve PVT solutions from it.
-
-CGGTTS: PVT and advanced Timing
-===============================
-
-The [rnx2cggtts application](https://github.com/georust/rinex) uses this
-solver to estimate local clock performances and allow the comparison
-of two remote clocks.
diff --git a/doc/config.md b/doc/config.md
new file mode 100644
index 0000000..6ea7808
--- /dev/null
+++ b/doc/config.md
@@ -0,0 +1,10 @@
+Solver configuration
+====================
+
+The `Cfg` structure is how to tune and adjust this position solver.  
+This solver is smart enough to generate more than decent results even if we rely on default settings
+(quick test / quick setups).
+
+Note that we always prefer precision over computation load.  
+In other words, for each mode the default options are already computation intensive
+and can be simplified.
diff --git a/doc/lsqspp.md b/doc/lsqspp.md
new file mode 100644
index 0000000..34771c2
--- /dev/null
+++ b/doc/lsqspp.md
@@ -0,0 +1,4 @@
+LSQSPP
+======
+
+TODO
diff --git a/doc/ppp.md b/doc/ppp.md
new file mode 100644
index 0000000..b9348de
--- /dev/null
+++ b/doc/ppp.md
@@ -0,0 +1,4 @@
+PPP
+===
+
+TODO
diff --git a/doc/spp.md b/doc/spp.md
new file mode 100644
index 0000000..8655af9
--- /dev/null
+++ b/doc/spp.md
@@ -0,0 +1,32 @@
+Single Point Positioning
+========================
+
+SPP can be deployed in degraded conditions where only one carrier signal and/or a unique constellation is in sight.
+
+Measurements requirements
+=========================
+
+`Mode::SPP` has less constraint on the measurements:
+you only need to gather the approriate amount of pseudo range measurements. 
+
+Several "bonus" may be unlocked if you provide extra measurements
+
+* Providing dopplers has the benefit to validate the interpolated SV state vector
+* Providing phase measurements will unlock fractional pseudo range estimate *in the future*
+* Providing measurements on another signal will unlock the true ionosphere delay to be estimated.
+
+Expected SPP configration
+=========================
+
+In `SPP` an interpolation of 9 makes sense and high orders will not be meaningful.  
+You can use that to reduce the computation load.
+
+Usually, in this mode the user provides measurements from a unique carrier signal.  
+In this context, only a modeling of the ionosphere impact is feasible.  
+As in other 
+
+- you can only hope for a precision of a few meters.    
+- an interpolation order of 9 makes sense, going above will increase the computation load without any benefit.
+- the solver will have to estimate the total bias due to Ionospheric delay. If you can define
+these components accurately yourself, you are encouranged to do it (see related API section and following paragraphs)
+

From e2aa6b78a85b1ababeb11edfaf3c1538c72fc0c5 Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Fri, 1 Dec 2023 17:58:27 +0100
Subject: [PATCH 04/16] working on residual validation

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/cfg.rs       | 86 +++++++++++++++++++++++++++++++-----------------
 src/solutions.rs | 10 ++++--
 src/solver.rs    | 56 +++++++++++++++++++++++--------
 3 files changed, 105 insertions(+), 47 deletions(-)

diff --git a/src/cfg.rs b/src/cfg.rs
index b94a3ba..1903b05 100644
--- a/src/cfg.rs
+++ b/src/cfg.rs
@@ -103,6 +103,35 @@ pub struct InternalDelay {
     pub frequency: f64,
 }
 
+#[derive(Default, Clone, Debug, PartialEq)]
+#[cfg_attr(feature = "serde", derive(Deserialize))]
+pub struct SolverOpts {
+    /// GDOP Threshold (max. limit) to invalidate ongoing GDOP
+    pub gdop_threshold: Option<f64>,
+    /// Weight Matrix in LSQ solving process
+    #[cfg_attr(feature = "serde", serde(default = "default_lsq_weight"))]
+    pub lsq_weight: Option<LSQWeight>,
+}
+   
+impl SolverOpts {
+    /*
+     * form the weight matrix to be used in the solving process
+     */
+    pub(crate) fn lsq_weight_matrix(&self, nb_rows: usize, sv_elev: Vec<f64>) -> DMatrix<f64> {
+        let mut mat = DMatrix::identity(sv_elev.len(), sv_elev.len());
+        match &self.lsq_weight {
+            Some(LSQWeight::LSQWeightCovar) => panic!("not implemented yet"),
+            Some(LSQWeight::LSQWeightMappingFunction(mapf)) => {
+                for i in 0..sv_elev.len() - 1 {
+                    mat[(i, i)] = mapf.a + mapf.b * ((-sv_elev[i]) / mapf.c).exp();
+                }
+            },
+            None => {},
+        }
+        mat
+    }
+}
+
 /// Atmospherical, Physical and Environmental modeling
 #[derive(Copy, Clone, Debug, PartialEq)]
 #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
@@ -171,9 +200,9 @@ pub struct Config {
     /// Internal delays
     #[cfg_attr(feature = "serde", serde(default))]
     pub int_delay: Vec<InternalDelay>,
-    /// Weight Matrix in LSQ solving process
-    #[cfg_attr(feature = "serde", serde(default = "default_lsq_weight"))]
-    pub lsq_weight: Option<LSQWeight>,
+    /// Solver customization
+    #[cfg_attr(feature = "serde", serde(default))]
+    pub solver: SolverOpts,
     /// Time Reference Delay, as defined by BIPM in
     /// "GPS Receivers Accurate Time Comparison" : the time delay
     /// between the GNSS receiver external reference clock and the sampling clock
@@ -203,8 +232,8 @@ pub struct Config {
 }
 
 impl Config {
-    pub fn default(solver: Mode) -> Self {
-        match solver {
+    pub fn default(mode: Mode) -> Self {
+        match mode {
             Mode::SPP => Self {
                 timescale: default_timescale(),
                 fixed_altitude: None,
@@ -217,7 +246,10 @@ impl Config {
                 max_sv: default_max_sv(),
                 int_delay: Default::default(),
                 externalref_delay: Default::default(),
-                lsq_weight: None,
+                solver: SolverOpts {
+                    lsq_weight: None,
+                    gdop_threshold: None,
+                },
             },
             Mode::LSQSPP => Self {
                 timescale: default_timescale(),
@@ -231,13 +263,15 @@ impl Config {
                 max_sv: default_max_sv(),
                 int_delay: Default::default(),
                 externalref_delay: Default::default(),
-                lsq_weight: Some(LSQWeight::LSQWeightMappingFunction(
-                    ElevationMappingFunction {
-                        a: 0.03,
-                        b: 0.03,
-                        c: 100.0,
-                    },
-                )),
+                solver: SolverOpts {
+                    gdop_threshold: Some(30.0),
+                    lsq_weight: Some(LSQWeight::LSQWeightMappingFunction(
+                        ElevationMappingFunction {
+                            a: 0.03,
+                            b: 0.03,
+                            c: 100.0,
+                        })),
+                },
             },
             Mode::PPP => Self {
                 timescale: default_timescale(),
@@ -252,25 +286,17 @@ impl Config {
                 max_sv: default_max_sv(),
                 int_delay: Default::default(),
                 externalref_delay: Default::default(),
-                lsq_weight: default_lsq_weight(),
-            },
-        }
-    }
-    /*
-     * form the weight matrix to be used in the solving process
-     */
-    pub(crate) fn lsq_weight_matrix(&self, nb_rows: usize, sv_elev: Vec<f64>) -> DMatrix<f64> {
-        let mut mat = DMatrix::identity(sv_elev.len(), sv_elev.len());
-        match &self.lsq_weight {
-            Some(LSQWeight::LSQWeightCovar) => panic!("not implemented yet"),
-            Some(LSQWeight::LSQWeightMappingFunction(mapf)) => {
-                for i in 0..sv_elev.len() - 1 {
-                    mat[(i, i)] = mapf.a + mapf.b * ((-sv_elev[i]) / mapf.c).exp();
-                }
+                solver: SolverOpts {
+                    gdop_threshold: None, //TODO
+                    lsq_weight: Some(LSQWeight::LSQWeightMappingFunction(
+                        ElevationMappingFunction {
+                            a: 0.03,
+                            b: 0.03,
+                            c: 100.0,
+                        })),
+                },
             },
-            None => {},
         }
-        mat
     }
 }
 
diff --git a/src/solutions.rs b/src/solutions.rs
index 6eb66c3..2c9a5a1 100644
--- a/src/solutions.rs
+++ b/src/solutions.rs
@@ -167,15 +167,19 @@ impl PVTSolution {
     pub fn sv(&self) -> Vec<SV> {
         self.sv.keys().copied().collect()
     }
-    /// Returns Horizontal Dilution of Precision of this PVT
+    /// Returns Geometrical Dilution of Precision of Self
+    pub fn gdop(&self) -> f64 {
+        (self.estimate.covar[(0, 0)] + self.estimate.covar[(1, 1)] + self.estimate.covar[(2,2)] + self.estimate.covar[(3, 3)]).sqrt()
+    }
+    /// Returns Horizontal Dilution of Precision of Self
     pub fn hdop(&self) -> f64 {
         (self.estimate.covar[(0, 0)] + self.estimate.covar[(1, 1)]).sqrt()
     }
-    /// Returns Vertical Dilution of Precision of this PVT
+    /// Returns Vertical Dilution of Precision of Self
     pub fn vdop(&self) -> f64 {
         self.estimate.covar[(2, 2)].sqrt()
     }
-    /// Returns Time Dilution of Precision of this PVT
+    /// Returns Time Dilution of Precision of Self
     pub fn tdop(&self) -> f64 {
         self.estimate.covar[(3, 3)].sqrt()
     }
diff --git a/src/solver.rs b/src/solver.rs
index 8ca441e..44d9155 100644
--- a/src/solver.rs
+++ b/src/solver.rs
@@ -17,10 +17,10 @@ use gnss::prelude::SV;
 use nalgebra::{DVector, Matrix3, Matrix4, Matrix4x1, MatrixXx4, Vector3};
 
 use crate::{
+    cfg::Config,
+    candidate::Candidate,
     apriori::AprioriPosition,
     bias::{IonosphericBias, TroposphericBias},
-    candidate::Candidate,
-    cfg::Config,
     solutions::{Estimate, PVTSVData, PVTSolution, PVTSolutionType},
 };
 
@@ -91,6 +91,8 @@ pub enum Error {
     PhysicalNonSenseRxPriorTx,
     #[error("physical non sense: t_rx is too late")]
     PhysicalNonSenseRxTooLate,
+    #[error("solution invalidated: gdop too large {0:.3E}")]
+    GDOPInvalidation(f64),
 }
 
 /// Interpolation result (state vector) that needs to be
@@ -154,12 +156,12 @@ where
     earth_frame: Frame,
     /* Sun frame */
     sun_frame: Frame,
-    /* prev. state vector for speed determination */
-    prev_sv_state: HashMap<SV, (Epoch, Vector3<f64>)>,
     /* prev. solution for internal logic */
     prev_pvt: Option<(Epoch, PVTSolution)>,
     /* prev. estimate for recursive impl. */
     prev_state: Option<Estimate>,
+    /* prev. state vector for internal velocity determination */
+    prev_sv_state: HashMap<SV, (Epoch, Vector3<f64>)>,
 }
 
 impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I> {
@@ -228,7 +230,9 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
             self.apriori.geodetic.z,
         );
 
+        let mode = self.mode;
         let modeling = self.cfg.modeling;
+        let solver_opts = &self.cfg.solver;
         let interp_order = self.cfg.interp_order;
 
         /* interpolate positions */
@@ -355,7 +359,7 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
                 };
                 if eclipsed {
                     debug!(
-                        "{:?} ({}): earth eclipsed, dropping",
+                        "{:?} ({}): dropped - eclipsed by earth",
                         pool[idx].t, pool[idx].sv
                     );
                     let _ = pool.swap_remove(idx);
@@ -380,7 +384,7 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
             if let Ok(sv_data) = cd.resolve(
                 t,
                 &self.cfg,
-                self.mode,
+                mode,
                 (x0, y0, z0),
                 (lat_ddeg, lon_ddeg, altitude_above_sea_m),
                 iono_bias,
@@ -393,8 +397,9 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
             }
         }
 
-        let w = self.cfg.lsq_weight_matrix(
-            nb_candidates,
+        let w = self.cfg.solver.lsq_weight_matrix(
+            //nb_candidates,
+            pvt_sv_data.len(), 
             pvt_sv_data
                 .iter()
                 .map(|(sv, sv_data)| sv_data.elevation)
@@ -402,11 +407,34 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         );
 
         let mut pvt_solution =
-            PVTSolution::new(g.clone(), w, y, pvt_sv_data, self.prev_state.clone())?;
+            PVTSolution::new(g.clone(), w.clone(), y.clone(), pvt_sv_data, self.prev_state.clone())?;
 
+        /*
+         * Solution validation : GDOP
+         */
+        if let Some(max_gdop) = solver_opts.gdop_threshold {
+            let gdop = pvt_solution.gdop();
+            if gdop > max_gdop {
+                // invalidate
+                return Err(Error::GDOPInvalidation(gdop));
+            }
+        }
+        
+        if mode.is_recursive() {
+            /*
+             * Solution validation : residual vector
+             */
+            let mut residuals = DVector::<f64>::zeros(nb_candidates);
+            for i in 0..nb_candidates {
+                residuals[i] = y[i] - 0.0_f64 / w[(i , i)];
+            }
+            // let denom = m - n - 1;
+            // residuals.clone().transpose() * residuals / denom;
+        }
+        
         if let Some((prev_t, prev_pvt)) = &self.prev_pvt {
             pvt_solution.v = (pvt_solution.p - prev_pvt.p) / (t - *prev_t).to_seconds();
-            if self.mode.is_recursive() {
+            if mode.is_recursive() {
                 self.prev_state = Some(pvt_solution.estimate.clone());
             }
         }
@@ -418,16 +446,16 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         if let Some(alt) = self.cfg.fixed_altitude {
             pvt_solution.p.z = self.apriori.ecef.z - alt;
             pvt_solution.v.z = 0.0_f64;
-            pvt_solution.estimate.covar[(2, 2)] = 0.0_f64;
+            // pvt_solution.estimate.covar[(2, 2)] = 0.0_f64;
         }
 
         match solution {
             PVTSolutionType::TimeOnly => {
                 pvt_solution.p = Vector3::<f64>::default();
                 pvt_solution.v = Vector3::<f64>::default();
-                pvt_solution.estimate.covar[(0, 0)] = 0.0_f64;
-                pvt_solution.estimate.covar[(1, 1)] = 0.0_f64;
-                pvt_solution.estimate.covar[(2, 2)] = 0.0_f64;
+                // pvt_solution.estimate.covar[(0, 0)] = 0.0_f64;
+                // pvt_solution.estimate.covar[(1, 1)] = 0.0_f64;
+                // pvt_solution.estimate.covar[(2, 2)] = 0.0_f64;
             },
             _ => {},
         }

From e3d033961240a08b7f34e831bb93e3a7887da896 Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Fri, 1 Dec 2023 17:58:55 +0100
Subject: [PATCH 05/16] lib update

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 Cargo.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Cargo.toml b/Cargo.toml
index be794d9..e3b8d4e 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "gnss-rtk"
-version = "0.3.0"
+version = "0.3.1"
 license = "MIT OR Apache-2.0"
 authors = ["Guillaume W. Bres <guillaume.bressaix@gmail.com>"]
 description = "GNSS position solver"

From 2b02d8072ae134644d84ef89f8e17669eec2c332 Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Fri, 1 Dec 2023 18:15:55 +0100
Subject: [PATCH 06/16] make sure we still have enough candidates

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/solver.rs | 32 +++++++++++++++++++++++++-------
 1 file changed, 25 insertions(+), 7 deletions(-)

diff --git a/src/solver.rs b/src/solver.rs
index 44d9155..6130d02 100644
--- a/src/solver.rs
+++ b/src/solver.rs
@@ -332,17 +332,36 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
 
         /* remove observed signals above snr mask (if any) */
         if let Some(min_snr) = self.cfg.min_snr {
+            let mut nb_removed :usize = 0;
             for idx in 0..pool.len() {
-                let (init_code, init_phase) = (pool[idx].code.len(), pool[idx].phase.len());
-                pool[idx].min_snr_mask(min_snr);
-                let delta_code = init_code - pool[idx].code.len();
-                let delta_phase = init_phase - pool[idx].phase.len();
+                let (init_code, init_phase) = (pool[idx - nb_removed].code.len(), pool[idx - nb_removed].phase.len());
+                pool[idx - nb_removed].min_snr_mask(min_snr);
+                let delta_code = init_code - pool[idx - nb_removed].code.len();
+                let delta_phase = init_phase - pool[idx - nb_removed].phase.len();
                 if delta_code > 0 || delta_phase > 0 {
                     debug!(
                         "{:?} ({}) : {} code | {} phase below snr mask",
-                        pool[idx].t, pool[idx].sv, delta_code, delta_phase
+                        pool[idx - nb_removed].t, pool[idx - nb_removed].sv, delta_code, delta_phase
                     );
                 }
+                /* make sure we're still compliant with the strategy */
+                match mode {
+                    Mode::SPP | Mode::LSQSPP => {
+                        if pool[idx - nb_removed].code.len() == 0 {
+                            debug!("{:?} ({}) dropped on bad snr", pool[idx].t, pool[idx].sv);
+                            let _ = pool.swap_remove(idx - nb_removed);
+                            nb_removed += 1;
+                        }
+                    },
+                    Mode::PPP => {
+                        let mut drop = !pool[idx - nb_removed].dual_freq_pseudorange();
+                        drop |= !pool[idx - nb_removed].dual_freq_phase();
+                        if drop {
+                            let _ = pool.swap_remove(idx - nb_removed);
+                            nb_removed += 1;
+                        }
+                    },
+                }
             }
         }
 
@@ -398,8 +417,7 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         }
 
         let w = self.cfg.solver.lsq_weight_matrix(
-            //nb_candidates,
-            pvt_sv_data.len(), 
+            nb_candidates,
             pvt_sv_data
                 .iter()
                 .map(|(sv, sv_data)| sv_data.elevation)

From c175e11816f0fb8c99cf340bdc1da0ca351919d2 Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Fri, 1 Dec 2023 18:18:05 +0100
Subject: [PATCH 07/16] test signal against mode compliance

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 Cargo.toml       |  1 +
 src/candidate.rs | 26 +++++++++++++++++++++++++-
 2 files changed, 26 insertions(+), 1 deletion(-)

diff --git a/Cargo.toml b/Cargo.toml
index e3b8d4e..5ca1f96 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -15,6 +15,7 @@ readme = "README.md"
 log = "0.4"
 thiserror = "1"
 map_3d = "0.1.5"
+itertools = "0.12.0"
 nalgebra = "0.32.3"
 nyx-space = "=2.0.0-beta.0" 
 gnss-rs = { version = "2.1.2", features = ["serde"] }
diff --git a/src/candidate.rs b/src/candidate.rs
index e22c506..f2f49c3 100644
--- a/src/candidate.rs
+++ b/src/candidate.rs
@@ -2,10 +2,13 @@
 
 use gnss::prelude::SV;
 use hifitime::Unit;
+use itertools::Itertools;
 use log::debug;
 use map_3d::deg2rad;
+
 use nyx::cosmic::SPEED_OF_LIGHT;
 use nyx::linalg::{DVector, MatrixXx4};
+use nyx::md::prelude::Frame;
 
 use crate::prelude::{Config, Duration, Epoch, InterpolationResult, Mode, Vector3};
 use crate::solutions::{PVTBias, PVTSVData};
@@ -14,7 +17,6 @@ use crate::{
     bias::{IonosphericBias, TropoModel, TroposphericBias},
     Error,
 };
-use nyx::md::prelude::Frame;
 
 /// Signal observation to attach to each candidate
 #[derive(Debug, Default, Clone)]
@@ -95,6 +97,28 @@ impl Candidate {
             // .map(|pr| pr.value)
             .reduce(|k, _| k)
     }
+    /*
+     * Returns true if we have pseudo range observ on two carriers
+     */
+    pub(crate) fn dual_freq_pseudorange(&self) -> bool {
+        self.code
+            .iter()
+            .map(|c| (c.frequency * 10.0) as u16)
+            .unique()
+            .count()
+            > 0
+    }
+    /*
+     * Returns true if we have phase observ on two carriers
+     */
+    pub(crate) fn dual_freq_phase(&self) -> bool {
+        self.phase
+            .iter()
+            .map(|c| (c.frequency * 10.0) as u16)
+            .unique()
+            .count()
+            > 0
+    }
     /*
      * apply min SNR mask
      */

From 135f09bb781e2c34cd531c3fe414dd89912b8b78 Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Fri, 1 Dec 2023 18:18:12 +0100
Subject: [PATCH 08/16] run linter

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/cfg.rs       |  8 +++++---
 src/solutions.rs |  6 +++++-
 src/solver.rs    | 31 +++++++++++++++++++++----------
 3 files changed, 31 insertions(+), 14 deletions(-)

diff --git a/src/cfg.rs b/src/cfg.rs
index 1903b05..56ee519 100644
--- a/src/cfg.rs
+++ b/src/cfg.rs
@@ -112,7 +112,7 @@ pub struct SolverOpts {
     #[cfg_attr(feature = "serde", serde(default = "default_lsq_weight"))]
     pub lsq_weight: Option<LSQWeight>,
 }
-   
+
 impl SolverOpts {
     /*
      * form the weight matrix to be used in the solving process
@@ -270,7 +270,8 @@ impl Config {
                             a: 0.03,
                             b: 0.03,
                             c: 100.0,
-                        })),
+                        },
+                    )),
                 },
             },
             Mode::PPP => Self {
@@ -293,7 +294,8 @@ impl Config {
                             a: 0.03,
                             b: 0.03,
                             c: 100.0,
-                        })),
+                        },
+                    )),
                 },
             },
         }
diff --git a/src/solutions.rs b/src/solutions.rs
index 2c9a5a1..3c3b8dc 100644
--- a/src/solutions.rs
+++ b/src/solutions.rs
@@ -169,7 +169,11 @@ impl PVTSolution {
     }
     /// Returns Geometrical Dilution of Precision of Self
     pub fn gdop(&self) -> f64 {
-        (self.estimate.covar[(0, 0)] + self.estimate.covar[(1, 1)] + self.estimate.covar[(2,2)] + self.estimate.covar[(3, 3)]).sqrt()
+        (self.estimate.covar[(0, 0)]
+            + self.estimate.covar[(1, 1)]
+            + self.estimate.covar[(2, 2)]
+            + self.estimate.covar[(3, 3)])
+            .sqrt()
     }
     /// Returns Horizontal Dilution of Precision of Self
     pub fn hdop(&self) -> f64 {
diff --git a/src/solver.rs b/src/solver.rs
index 6130d02..587ea56 100644
--- a/src/solver.rs
+++ b/src/solver.rs
@@ -17,10 +17,10 @@ use gnss::prelude::SV;
 use nalgebra::{DVector, Matrix3, Matrix4, Matrix4x1, MatrixXx4, Vector3};
 
 use crate::{
-    cfg::Config,
-    candidate::Candidate,
     apriori::AprioriPosition,
     bias::{IonosphericBias, TroposphericBias},
+    candidate::Candidate,
+    cfg::Config,
     solutions::{Estimate, PVTSVData, PVTSolution, PVTSolutionType},
 };
 
@@ -332,16 +332,22 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
 
         /* remove observed signals above snr mask (if any) */
         if let Some(min_snr) = self.cfg.min_snr {
-            let mut nb_removed :usize = 0;
+            let mut nb_removed: usize = 0;
             for idx in 0..pool.len() {
-                let (init_code, init_phase) = (pool[idx - nb_removed].code.len(), pool[idx - nb_removed].phase.len());
+                let (init_code, init_phase) = (
+                    pool[idx - nb_removed].code.len(),
+                    pool[idx - nb_removed].phase.len(),
+                );
                 pool[idx - nb_removed].min_snr_mask(min_snr);
                 let delta_code = init_code - pool[idx - nb_removed].code.len();
                 let delta_phase = init_phase - pool[idx - nb_removed].phase.len();
                 if delta_code > 0 || delta_phase > 0 {
                     debug!(
                         "{:?} ({}) : {} code | {} phase below snr mask",
-                        pool[idx - nb_removed].t, pool[idx - nb_removed].sv, delta_code, delta_phase
+                        pool[idx - nb_removed].t,
+                        pool[idx - nb_removed].sv,
+                        delta_code,
+                        delta_phase
                     );
                 }
                 /* make sure we're still compliant with the strategy */
@@ -424,8 +430,13 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
                 .collect(),
         );
 
-        let mut pvt_solution =
-            PVTSolution::new(g.clone(), w.clone(), y.clone(), pvt_sv_data, self.prev_state.clone())?;
+        let mut pvt_solution = PVTSolution::new(
+            g.clone(),
+            w.clone(),
+            y.clone(),
+            pvt_sv_data,
+            self.prev_state.clone(),
+        )?;
 
         /*
          * Solution validation : GDOP
@@ -437,19 +448,19 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
                 return Err(Error::GDOPInvalidation(gdop));
             }
         }
-        
+
         if mode.is_recursive() {
             /*
              * Solution validation : residual vector
              */
             let mut residuals = DVector::<f64>::zeros(nb_candidates);
             for i in 0..nb_candidates {
-                residuals[i] = y[i] - 0.0_f64 / w[(i , i)];
+                residuals[i] = y[i] - 0.0_f64 / w[(i, i)];
             }
             // let denom = m - n - 1;
             // residuals.clone().transpose() * residuals / denom;
         }
-        
+
         if let Some((prev_t, prev_pvt)) = &self.prev_pvt {
             pvt_solution.v = (pvt_solution.p - prev_pvt.p) / (t - *prev_t).to_seconds();
             if mode.is_recursive() {

From d8195f883146825e3e1623afebb312b54d606909 Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Sat, 2 Dec 2023 21:53:34 +0100
Subject: [PATCH 09/16] working on code residuals

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/cfg.rs       | 38 +++++++++++++-------------
 src/solutions.rs | 69 +++++++++++++++++++++++++++++++++++-------------
 src/solver.rs    | 53 ++++++++++++++++++++++++++++---------
 3 files changed, 109 insertions(+), 51 deletions(-)

diff --git a/src/cfg.rs b/src/cfg.rs
index 56ee519..f80ce31 100644
--- a/src/cfg.rs
+++ b/src/cfg.rs
@@ -84,7 +84,17 @@ fn default_relativistic_path_range() -> bool {
 }
 
 fn default_lsq_weight() -> Option<LSQWeight> {
-    None
+    Some(LSQWeight::LSQWeightMappingFunction(
+        ElevationMappingFunction {
+            a: 5.0,
+            b: 0.0,
+            c: 10.0,
+        },
+    ))
+}
+
+fn default_gdop_threshold() -> Option<f64> {
+    Some(0.15)
 }
 
 #[derive(Default, Debug, Clone, PartialEq)]
@@ -107,6 +117,7 @@ pub struct InternalDelay {
 #[cfg_attr(feature = "serde", derive(Deserialize))]
 pub struct SolverOpts {
     /// GDOP Threshold (max. limit) to invalidate ongoing GDOP
+    #[cfg_attr(feature = "serde", serde(default = "default_gdop_threshold"))]
     pub gdop_threshold: Option<f64>,
     /// Weight Matrix in LSQ solving process
     #[cfg_attr(feature = "serde", serde(default = "default_lsq_weight"))]
@@ -123,7 +134,8 @@ impl SolverOpts {
             Some(LSQWeight::LSQWeightCovar) => panic!("not implemented yet"),
             Some(LSQWeight::LSQWeightMappingFunction(mapf)) => {
                 for i in 0..sv_elev.len() - 1 {
-                    mat[(i, i)] = mapf.a + mapf.b * ((-sv_elev[i]) / mapf.c).exp();
+                    let sigma = mapf.a + mapf.b * ((-sv_elev[i]) / mapf.c).exp();
+                    mat[(i, i)] = 1.0 / sigma.powi(2);
                 }
             },
             None => {},
@@ -248,7 +260,7 @@ impl Config {
                 externalref_delay: Default::default(),
                 solver: SolverOpts {
                     lsq_weight: None,
-                    gdop_threshold: None,
+                    gdop_threshold: default_gdop_threshold(),
                 },
             },
             Mode::LSQSPP => Self {
@@ -264,14 +276,8 @@ impl Config {
                 int_delay: Default::default(),
                 externalref_delay: Default::default(),
                 solver: SolverOpts {
-                    gdop_threshold: Some(30.0),
-                    lsq_weight: Some(LSQWeight::LSQWeightMappingFunction(
-                        ElevationMappingFunction {
-                            a: 0.03,
-                            b: 0.03,
-                            c: 100.0,
-                        },
-                    )),
+                    gdop_threshold: default_gdop_threshold(),
+                    lsq_weight: default_lsq_weight(),
                 },
             },
             Mode::PPP => Self {
@@ -288,14 +294,8 @@ impl Config {
                 int_delay: Default::default(),
                 externalref_delay: Default::default(),
                 solver: SolverOpts {
-                    gdop_threshold: None, //TODO
-                    lsq_weight: Some(LSQWeight::LSQWeightMappingFunction(
-                        ElevationMappingFunction {
-                            a: 0.03,
-                            b: 0.03,
-                            c: 100.0,
-                        },
-                    )),
+                    gdop_threshold: default_gdop_threshold(),
+                    lsq_weight: default_lsq_weight(),
                 },
             },
         }
diff --git a/src/solutions.rs b/src/solutions.rs
index 3c3b8dc..7a13815 100644
--- a/src/solutions.rs
+++ b/src/solutions.rs
@@ -27,7 +27,7 @@ impl std::fmt::Display for PVTSolutionType {
     }
 }
 
-use nalgebra::base::{DMatrix, DVector, Matrix1x4, Matrix4, Matrix4x1, MatrixXx4};
+use nalgebra::base::{DMatrix, DVector, Matrix1x4, Matrix3, Matrix4, Matrix4x1, MatrixXx4};
 use nyx::cosmic::SPEED_OF_LIGHT;
 
 #[cfg(feature = "serde")]
@@ -87,10 +87,10 @@ pub struct PVTSVData {
 
 #[derive(Default, Clone, Debug)]
 pub struct Estimate {
-    /* x: estimate */
+    /* x estimate */
     pub(crate) x: Matrix4x1<f64>,
-    /* and its covariance matrix */
-    pub(crate) covar: Matrix4<f64>,
+    /* Q matrix */
+    pub(crate) q: Matrix4<f64>,
 }
 
 /// PVT Solution, always expressed as the correction to apply
@@ -132,21 +132,21 @@ impl PVTSolution {
                     .try_inverse()
                     .ok_or(Error::MatrixInversionError)?;
                 let x = (q * g_prime.clone()) * w.clone() * y;
-                Estimate { x, covar: q }
+                Estimate { x, q }
             },
             Some(state) => {
                 let p_1 = state
-                    .covar
+                    .q
                     .try_inverse()
                     .ok_or(Error::CovarMatrixInversionError)?;
 
                 let q = g_prime.clone() * w.clone() * g.clone();
-                let covar = (p_1 + q)
+                let q = (p_1 + q)
                     .try_inverse()
                     .ok_or(Error::CovarMatrixInversionError)?;
 
-                let x = covar * (p_1 * state.x + (g_prime.clone() * w.clone() * y));
-                Estimate { x, covar }
+                let x = q * (p_1 * state.x + (g_prime.clone() * w.clone() * y));
+                Estimate { x, q }
             },
         };
 
@@ -167,24 +167,55 @@ impl PVTSolution {
     pub fn sv(&self) -> Vec<SV> {
         self.sv.keys().copied().collect()
     }
-    /// Returns Geometrical Dilution of Precision of Self
+    /// Returns Geometric Dilution of Precision of Self
     pub fn gdop(&self) -> f64 {
-        (self.estimate.covar[(0, 0)]
-            + self.estimate.covar[(1, 1)]
-            + self.estimate.covar[(2, 2)]
-            + self.estimate.covar[(3, 3)])
+        (self.estimate.q[(0, 0)]
+            + self.estimate.q[(1, 1)]
+            + self.estimate.q[(2, 2)]
+            + self.estimate.q[(3, 3)])
             .sqrt()
     }
+    /// Returns Position Diultion of Precision of Self
+    pub fn pdop(&self) -> f64 {
+        (self.estimate.q[(0, 0)] + self.estimate.q[(1, 1)] + self.estimate.q[(2, 2)]).sqrt()
+    }
+    fn q_enu(&self, lat: f64, lon: f64) -> Matrix3<f64> {
+        let r = Matrix3::<f64>::new(
+            -lon.sin(),
+            -lon.cos() * lat.sin(),
+            lat.cos() * lon.cos(),
+            lon.cos(),
+            -lat.sin() * lon.sin(),
+            lat.cos() * lon.sin(),
+            0.0_f64,
+            lat.cos(),
+            lon.sin(),
+        );
+        let q_3 = Matrix3::<f64>::new(
+            self.estimate.q[(0, 0)],
+            self.estimate.q[(0, 1)],
+            self.estimate.q[(0, 2)],
+            self.estimate.q[(1, 0)],
+            self.estimate.q[(1, 1)],
+            self.estimate.q[(1, 2)],
+            self.estimate.q[(2, 0)],
+            self.estimate.q[(2, 1)],
+            self.estimate.q[(2, 2)],
+        );
+        r.clone().transpose() * q_3 * r.clone()
+    }
     /// Returns Horizontal Dilution of Precision of Self
-    pub fn hdop(&self) -> f64 {
-        (self.estimate.covar[(0, 0)] + self.estimate.covar[(1, 1)]).sqrt()
+    pub fn hdop(&self, lat: f64, lon: f64) -> f64 {
+        let q_enu = self.q_enu(lat, lon);
+        (q_enu[(0, 0)] + q_enu[(1, 1)]).sqrt()
     }
     /// Returns Vertical Dilution of Precision of Self
-    pub fn vdop(&self) -> f64 {
-        self.estimate.covar[(2, 2)].sqrt()
+    pub fn vdop(&self, lat: f64, lon: f64) -> f64 {
+        let q_enu = self.q_enu(lat, lon);
+        q_enu[(1, 1)].sqrt()
     }
     /// Returns Time Dilution of Precision of Self
     pub fn tdop(&self) -> f64 {
-        self.estimate.covar[(3, 3)].sqrt()
+        self.estimate.q[(3, 3)].sqrt()
     }
 }
diff --git a/src/solver.rs b/src/solver.rs
index 587ea56..7cf1521 100644
--- a/src/solver.rs
+++ b/src/solver.rs
@@ -434,10 +434,18 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
             g.clone(),
             w.clone(),
             y.clone(),
-            pvt_sv_data,
+            pvt_sv_data.clone(),
             self.prev_state.clone(),
         )?;
 
+        if mode.is_recursive() {
+            self.prev_state = Some(pvt_solution.estimate.clone());
+        }
+
+        if let Some((prev_t, prev_pvt)) = &self.prev_pvt {
+            pvt_solution.v = (pvt_solution.p - prev_pvt.p) / (t - *prev_t).to_seconds();
+        }
+
         /*
          * Solution validation : GDOP
          */
@@ -451,21 +459,40 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
 
         if mode.is_recursive() {
             /*
-             * Solution validation : residual vector
+             * Solution validation : code residual
              */
-            let mut residuals = DVector::<f64>::zeros(nb_candidates);
-            for i in 0..nb_candidates {
-                residuals[i] = y[i] - 0.0_f64 / w[(i, i)];
+            let mut residuals = DVector::<f64>::zeros(pool.len());
+            for (idx, cd) in pool.iter().enumerate() {
+                let sv = pvt_sv_data
+                    .iter()
+                    .filter_map(|(sv, data)| if *sv == cd.sv { Some(data) } else { None })
+                    .reduce(|k, _| k)
+                    .unwrap();
+
+                let pr = cd.prefered_pseudorange().unwrap().value;
+                let state = cd.state.unwrap().position;
+                let estimate = self.apriori.ecef + pvt_solution.p;
+
+                let (sv_x, sv_y, sv_z) = (state[0], state[1], state[2]);
+                let (x, y, z) = (estimate[0], estimate[1], estimate[2]);
+
+                let rho = ((sv_x - x).powi(2) + (sv_y - y).powi(2) + (sv_z - z).powi(2)).sqrt();
+
+                let c_dt = (cd.clock_corr.to_seconds() - pvt_solution.dt) * SPEED_OF_LIGHT;
+
+                residuals[idx] = pr;
+                residuals[idx] -= rho;
+                residuals[idx] += c_dt;
+                residuals[idx] -= sv.tropo_bias.value().unwrap_or(0.0);
+                residuals[idx] -= sv.iono_bias.value().unwrap_or(0.0);
+                residuals[idx] /= w[(idx, idx)];
+                debug!("{:?}: coderes={}/w={}", t, residuals[idx], w[(idx, idx)]);
             }
-            // let denom = m - n - 1;
-            // residuals.clone().transpose() * residuals / denom;
-        }
 
-        if let Some((prev_t, prev_pvt)) = &self.prev_pvt {
-            pvt_solution.v = (pvt_solution.p - prev_pvt.p) / (t - *prev_t).to_seconds();
-            if mode.is_recursive() {
-                self.prev_state = Some(pvt_solution.estimate.clone());
-            }
+            /*
+             * If stablized: validate new innovation
+             */
+            for i in 0..nb_candidates {}
         }
 
         /*

From f4e98da8a9bbb7bf5dbdca33b10f24625ebf257b Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Sun, 3 Dec 2023 11:39:00 +0100
Subject: [PATCH 10/16] working on solution validation

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/cfg.rs       |  2 +-
 src/solutions.rs |  2 +-
 src/solver.rs    | 99 +++++++++++++++++++++++++++++++-----------------
 3 files changed, 66 insertions(+), 37 deletions(-)

diff --git a/src/cfg.rs b/src/cfg.rs
index f80ce31..95cacb8 100644
--- a/src/cfg.rs
+++ b/src/cfg.rs
@@ -3,8 +3,8 @@ use thiserror::Error;
 #[cfg(feature = "serde")]
 use serde::{Deserialize, Serialize};
 
+use nalgebra::DMatrix; 
 use crate::prelude::{Mode, TimeScale};
-use nalgebra::{DMatrix, DVector, Matrix3, MatrixXx4, Vector3};
 
 /// Configuration Error
 #[derive(Debug, Error)]
diff --git a/src/solutions.rs b/src/solutions.rs
index 7a13815..0313b71 100644
--- a/src/solutions.rs
+++ b/src/solutions.rs
@@ -27,7 +27,7 @@ impl std::fmt::Display for PVTSolutionType {
     }
 }
 
-use nalgebra::base::{DMatrix, DVector, Matrix1x4, Matrix3, Matrix4, Matrix4x1, MatrixXx4};
+use nalgebra::base::{DMatrix, DVector, Matrix3, Matrix4, Matrix4x1, MatrixXx4};
 use nyx::cosmic::SPEED_OF_LIGHT;
 
 #[cfg(feature = "serde")]
diff --git a/src/solver.rs b/src/solver.rs
index 7cf1521..fc10065 100644
--- a/src/solver.rs
+++ b/src/solver.rs
@@ -14,7 +14,7 @@ use nyx::md::prelude::{Bodies, Frame, LightTimeCalc};
 
 use gnss::prelude::SV;
 
-use nalgebra::{DVector, Matrix3, Matrix4, Matrix4x1, MatrixXx4, Vector3};
+use nalgebra::{DVector, Vector3, DMatrix, Matrix3, MatrixXx4};
 
 use crate::{
     apriori::AprioriPosition,
@@ -458,41 +458,15 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         }
 
         if mode.is_recursive() {
-            /*
-             * Solution validation : code residual
-             */
-            let mut residuals = DVector::<f64>::zeros(pool.len());
-            for (idx, cd) in pool.iter().enumerate() {
-                let sv = pvt_sv_data
-                    .iter()
-                    .filter_map(|(sv, data)| if *sv == cd.sv { Some(data) } else { None })
-                    .reduce(|k, _| k)
-                    .unwrap();
-
-                let pr = cd.prefered_pseudorange().unwrap().value;
-                let state = cd.state.unwrap().position;
-                let estimate = self.apriori.ecef + pvt_solution.p;
-
-                let (sv_x, sv_y, sv_z) = (state[0], state[1], state[2]);
-                let (x, y, z) = (estimate[0], estimate[1], estimate[2]);
-
-                let rho = ((sv_x - x).powi(2) + (sv_y - y).powi(2) + (sv_z - z).powi(2)).sqrt();
-
-                let c_dt = (cd.clock_corr.to_seconds() - pvt_solution.dt) * SPEED_OF_LIGHT;
-
-                residuals[idx] = pr;
-                residuals[idx] -= rho;
-                residuals[idx] += c_dt;
-                residuals[idx] -= sv.tropo_bias.value().unwrap_or(0.0);
-                residuals[idx] -= sv.iono_bias.value().unwrap_or(0.0);
-                residuals[idx] /= w[(idx, idx)];
-                debug!("{:?}: coderes={}/w={}", t, residuals[idx], w[(idx, idx)]);
-            }
+            let validator = SolutionValidator::new(
+                &self.apriori.ecef, 
+                &pool, 
+                &w, 
+                &pvt_solution);
+            
+            if !validator.valid().is_ok() {
 
-            /*
-             * If stablized: validate new innovation
-             */
-            for i in 0..nb_candidates {}
+            }
         }
 
         /*
@@ -553,3 +527,58 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         }
     }
 }
+
+pub enum SolutionInvalidation {
+    InnovationOutlier,
+    CodeResidual,
+}
+
+struct SolutionValidator {
+    residuals: DVector::<f64>,
+}
+
+impl SolutionValidator {
+    fn new(
+        apriori_ecef: &Vector3<f64>,
+        pool: &Vec<Candidate>, 
+        w: &DMatrix<f64>,
+        solution: &PVTSolution,
+    ) -> Self {
+            
+        let mut residuals = DVector::<f64>::zeros(pool.len());
+
+        for (idx, cd) in pool.iter().enumerate() {
+            let sv = solution.sv
+                .iter()
+                .filter_map(|(sv, data)| if *sv == cd.sv { Some(data) } else { None })
+                .reduce(|k, _| k)
+                .unwrap();
+
+            let pr = cd.prefered_pseudorange().unwrap().value;
+            let state = cd.state.unwrap().position;
+            let estimate = apriori_ecef + solution.p;
+
+            let (sv_x, sv_y, sv_z) = (state[0], state[1], state[2]);
+            let (x, y, z) = (estimate[0], estimate[1], estimate[2]);
+
+            let rho = ((sv_x - x).powi(2) + (sv_y - y).powi(2) + (sv_z - z).powi(2)).sqrt();
+
+            let dt = cd.clock_corr.to_seconds() - solution.dt;
+
+            residuals[idx] = pr;
+            residuals[idx] -= rho;
+            residuals[idx] += dt * SPEED_OF_LIGHT;
+            residuals[idx] -= sv.tropo_bias.value().unwrap_or(0.0);
+            residuals[idx] -= sv.iono_bias.value().unwrap_or(0.0);
+            residuals[idx] /= w[(idx, idx)];
+            debug!("{:?} ({}): coderes={}/w={}", cd.t, cd.sv, residuals[idx], w[(idx, idx)]);
+        }
+        Self { residuals }
+    }
+    /*
+     * Solution validation process
+     */
+    fn valid(&self) -> Result<(), SolutionInvalidation> {
+        Ok(())
+    }
+}

From 8c49610bdb8ca73bc933f8665ba1874530d1b73d Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Sun, 3 Dec 2023 12:03:19 +0100
Subject: [PATCH 11/16] introduce solution validator

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/candidate.rs                       |  1 -
 src/cfg.rs                             | 14 +++-
 src/{solutions.rs => solutions/mod.rs} |  2 +
 src/solutions/validator.rs             | 90 ++++++++++++++++++++++++
 src/solver.rs                          | 95 ++++----------------------
 5 files changed, 116 insertions(+), 86 deletions(-)
 rename src/{solutions.rs => solutions/mod.rs} (99%)
 create mode 100644 src/solutions/validator.rs

diff --git a/src/candidate.rs b/src/candidate.rs
index f2f49c3..c44ffaf 100644
--- a/src/candidate.rs
+++ b/src/candidate.rs
@@ -8,7 +8,6 @@ use map_3d::deg2rad;
 
 use nyx::cosmic::SPEED_OF_LIGHT;
 use nyx::linalg::{DVector, MatrixXx4};
-use nyx::md::prelude::Frame;
 
 use crate::prelude::{Config, Duration, Epoch, InterpolationResult, Mode, Vector3};
 use crate::solutions::{PVTBias, PVTSVData};
diff --git a/src/cfg.rs b/src/cfg.rs
index 95cacb8..c3987d7 100644
--- a/src/cfg.rs
+++ b/src/cfg.rs
@@ -3,8 +3,8 @@ use thiserror::Error;
 #[cfg(feature = "serde")]
 use serde::{Deserialize, Serialize};
 
-use nalgebra::DMatrix; 
 use crate::prelude::{Mode, TimeScale};
+use nalgebra::DMatrix;
 
 /// Configuration Error
 #[derive(Debug, Error)]
@@ -97,6 +97,10 @@ fn default_gdop_threshold() -> Option<f64> {
     Some(0.15)
 }
 
+fn default_innov_threshold() -> Option<f64> {
+    None
+}
+
 #[derive(Default, Debug, Clone, PartialEq)]
 #[cfg_attr(feature = "serde", derive(Deserialize))]
 /// System Internal Delay as defined by BIPM in
@@ -122,6 +126,9 @@ pub struct SolverOpts {
     /// Weight Matrix in LSQ solving process
     #[cfg_attr(feature = "serde", serde(default = "default_lsq_weight"))]
     pub lsq_weight: Option<LSQWeight>,
+    /// Threshold on new filter innovation
+    #[cfg_attr(feature = "serde", serde(default = "default_innov_threshold"))]
+    pub innovation_threshold: Option<f64>,
 }
 
 impl SolverOpts {
@@ -259,8 +266,9 @@ impl Config {
                 int_delay: Default::default(),
                 externalref_delay: Default::default(),
                 solver: SolverOpts {
-                    lsq_weight: None,
                     gdop_threshold: default_gdop_threshold(),
+                    lsq_weight: None,
+                    innovation_threshold: None,
                 },
             },
             Mode::LSQSPP => Self {
@@ -278,6 +286,7 @@ impl Config {
                 solver: SolverOpts {
                     gdop_threshold: default_gdop_threshold(),
                     lsq_weight: default_lsq_weight(),
+                    innovation_threshold: default_innov_threshold(),
                 },
             },
             Mode::PPP => Self {
@@ -296,6 +305,7 @@ impl Config {
                 solver: SolverOpts {
                     gdop_threshold: default_gdop_threshold(),
                     lsq_weight: default_lsq_weight(),
+                    innovation_threshold: default_innov_threshold(),
                 },
             },
         }
diff --git a/src/solutions.rs b/src/solutions/mod.rs
similarity index 99%
rename from src/solutions.rs
rename to src/solutions/mod.rs
index 0313b71..a7f7b21 100644
--- a/src/solutions.rs
+++ b/src/solutions/mod.rs
@@ -3,6 +3,8 @@ use crate::prelude::{Vector3, SV};
 use crate::Error;
 use std::collections::HashMap;
 
+pub(crate) mod validator;
+
 #[derive(Debug, Copy, Clone, Default)]
 pub enum PVTSolutionType {
     /// Default, complete solution with Position,
diff --git a/src/solutions/validator.rs b/src/solutions/validator.rs
new file mode 100644
index 0000000..5cb933f
--- /dev/null
+++ b/src/solutions/validator.rs
@@ -0,0 +1,90 @@
+use log::debug;
+use nalgebra::{DMatrix, DVector, Vector3};
+use nyx::cosmic::SPEED_OF_LIGHT;
+use thiserror::Error;
+
+use crate::{
+    cfg::SolverOpts,
+    prelude::{Candidate, PVTSolution},
+};
+
+#[derive(Clone, Debug, Error)]
+pub enum SolutionInvalidation {
+    #[error("gdop limit exceeded {0}")]
+    GDOPOutlier(f64),
+    #[error("innovation outlier |{0}|")]
+    InnovationOutlier(f64),
+    #[error("coderes limit exceeded {0}")]
+    CodeResidual(f64),
+}
+
+pub(crate) struct SolutionValidator {
+    gdop: f64,
+    residuals: DVector<f64>,
+}
+
+impl SolutionValidator {
+    pub fn new(
+        apriori_ecef: &Vector3<f64>,
+        pool: &Vec<Candidate>,
+        w: &DMatrix<f64>,
+        solution: &PVTSolution,
+    ) -> Self {
+        let gdop = solution.gdop();
+        let mut residuals = DVector::<f64>::zeros(pool.len());
+
+        for (idx, cd) in pool.iter().enumerate() {
+            let sv = solution
+                .sv
+                .iter()
+                .filter_map(|(sv, data)| if *sv == cd.sv { Some(data) } else { None })
+                .reduce(|k, _| k)
+                .unwrap();
+
+            let pr = cd.prefered_pseudorange().unwrap().value;
+            let state = cd.state.unwrap().position;
+            let estimate = apriori_ecef + solution.p;
+
+            let (sv_x, sv_y, sv_z) = (state[0], state[1], state[2]);
+            let (x, y, z) = (estimate[0], estimate[1], estimate[2]);
+
+            let rho = ((sv_x - x).powi(2) + (sv_y - y).powi(2) + (sv_z - z).powi(2)).sqrt();
+
+            let dt = cd.clock_corr.to_seconds() - solution.dt;
+
+            residuals[idx] = pr;
+            residuals[idx] -= rho;
+            residuals[idx] += dt * SPEED_OF_LIGHT;
+            residuals[idx] -= sv.tropo_bias.value().unwrap_or(0.0);
+            residuals[idx] -= sv.iono_bias.value().unwrap_or(0.0);
+            residuals[idx] /= w[(idx, idx)];
+            debug!(
+                "{:?} ({}): coderes={}/w={}",
+                cd.t,
+                cd.sv,
+                residuals[idx],
+                w[(idx, idx)]
+            );
+        }
+        Self { residuals, gdop }
+    }
+    /*
+     * Solution validation process
+     */
+    pub fn valid(&self, opts: &SolverOpts) -> Result<(), SolutionInvalidation> {
+        if let Some(max_gdop) = opts.gdop_threshold {
+            if self.gdop > max_gdop {
+                return Err(SolutionInvalidation::GDOPOutlier(self.gdop));
+            }
+        }
+        if let Some(threshold) = opts.innovation_threshold {
+            for idx in 0..self.residuals.len() {
+                let res = self.residuals[idx].abs();
+                if res > threshold {
+                    return Err(SolutionInvalidation::InnovationOutlier(res));
+                }
+            }
+        }
+        Ok(())
+    }
+}
diff --git a/src/solver.rs b/src/solver.rs
index fc10065..8fef2e2 100644
--- a/src/solver.rs
+++ b/src/solver.rs
@@ -14,14 +14,17 @@ use nyx::md::prelude::{Bodies, Frame, LightTimeCalc};
 
 use gnss::prelude::SV;
 
-use nalgebra::{DVector, Vector3, DMatrix, Matrix3, MatrixXx4};
+use nalgebra::{DMatrix, DVector, Matrix3, MatrixXx4, Vector3};
 
 use crate::{
     apriori::AprioriPosition,
     bias::{IonosphericBias, TroposphericBias},
     candidate::Candidate,
     cfg::Config,
-    solutions::{Estimate, PVTSVData, PVTSolution, PVTSolutionType},
+    solutions::{
+        validator::{SolutionInvalidation, SolutionValidator},
+        Estimate, PVTSVData, PVTSolution, PVTSolutionType,
+    },
 };
 
 /// Solving mode
@@ -91,8 +94,8 @@ pub enum Error {
     PhysicalNonSenseRxPriorTx,
     #[error("physical non sense: t_rx is too late")]
     PhysicalNonSenseRxTooLate,
-    #[error("solution invalidated: gdop too large {0:.3E}")]
-    GDOPInvalidation(f64),
+    #[error("invalidated solution: {0}")]
+    InvalidatedSolution(SolutionInvalidation),
 }
 
 /// Interpolation result (state vector) that needs to be
@@ -446,27 +449,13 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
             pvt_solution.v = (pvt_solution.p - prev_pvt.p) / (t - *prev_t).to_seconds();
         }
 
-        /*
-         * Solution validation : GDOP
-         */
-        if let Some(max_gdop) = solver_opts.gdop_threshold {
-            let gdop = pvt_solution.gdop();
-            if gdop > max_gdop {
-                // invalidate
-                return Err(Error::GDOPInvalidation(gdop));
-            }
-        }
+        self.prev_pvt = Some((t, pvt_solution.clone()));
 
-        if mode.is_recursive() {
-            let validator = SolutionValidator::new(
-                &self.apriori.ecef, 
-                &pool, 
-                &w, 
-                &pvt_solution);
-            
-            if !validator.valid().is_ok() {
+        let validator = SolutionValidator::new(&self.apriori.ecef, &pool, &w, &pvt_solution);
 
-            }
+        let valid = validator.valid(solver_opts);
+        if !valid.is_ok() {
+            return Err(Error::InvalidatedSolution(valid.err().unwrap()));
         }
 
         /*
@@ -476,21 +465,16 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         if let Some(alt) = self.cfg.fixed_altitude {
             pvt_solution.p.z = self.apriori.ecef.z - alt;
             pvt_solution.v.z = 0.0_f64;
-            // pvt_solution.estimate.covar[(2, 2)] = 0.0_f64;
         }
 
         match solution {
             PVTSolutionType::TimeOnly => {
                 pvt_solution.p = Vector3::<f64>::default();
                 pvt_solution.v = Vector3::<f64>::default();
-                // pvt_solution.estimate.covar[(0, 0)] = 0.0_f64;
-                // pvt_solution.estimate.covar[(1, 1)] = 0.0_f64;
-                // pvt_solution.estimate.covar[(2, 2)] = 0.0_f64;
             },
             _ => {},
         }
 
-        self.prev_pvt = Some((t, pvt_solution.clone()));
         Ok((t, pvt_solution))
     }
     /*
@@ -527,58 +511,3 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         }
     }
 }
-
-pub enum SolutionInvalidation {
-    InnovationOutlier,
-    CodeResidual,
-}
-
-struct SolutionValidator {
-    residuals: DVector::<f64>,
-}
-
-impl SolutionValidator {
-    fn new(
-        apriori_ecef: &Vector3<f64>,
-        pool: &Vec<Candidate>, 
-        w: &DMatrix<f64>,
-        solution: &PVTSolution,
-    ) -> Self {
-            
-        let mut residuals = DVector::<f64>::zeros(pool.len());
-
-        for (idx, cd) in pool.iter().enumerate() {
-            let sv = solution.sv
-                .iter()
-                .filter_map(|(sv, data)| if *sv == cd.sv { Some(data) } else { None })
-                .reduce(|k, _| k)
-                .unwrap();
-
-            let pr = cd.prefered_pseudorange().unwrap().value;
-            let state = cd.state.unwrap().position;
-            let estimate = apriori_ecef + solution.p;
-
-            let (sv_x, sv_y, sv_z) = (state[0], state[1], state[2]);
-            let (x, y, z) = (estimate[0], estimate[1], estimate[2]);
-
-            let rho = ((sv_x - x).powi(2) + (sv_y - y).powi(2) + (sv_z - z).powi(2)).sqrt();
-
-            let dt = cd.clock_corr.to_seconds() - solution.dt;
-
-            residuals[idx] = pr;
-            residuals[idx] -= rho;
-            residuals[idx] += dt * SPEED_OF_LIGHT;
-            residuals[idx] -= sv.tropo_bias.value().unwrap_or(0.0);
-            residuals[idx] -= sv.iono_bias.value().unwrap_or(0.0);
-            residuals[idx] /= w[(idx, idx)];
-            debug!("{:?} ({}): coderes={}/w={}", cd.t, cd.sv, residuals[idx], w[(idx, idx)]);
-        }
-        Self { residuals }
-    }
-    /*
-     * Solution validation process
-     */
-    fn valid(&self) -> Result<(), SolutionInvalidation> {
-        Ok(())
-    }
-}

From 47b4ed1bb96a6af915d443e33051c06278745d87 Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Sun, 3 Dec 2023 13:29:13 +0100
Subject: [PATCH 12/16] renamed positioning mode

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/cfg.rs | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/cfg.rs b/src/cfg.rs
index c3987d7..8ee34e8 100644
--- a/src/cfg.rs
+++ b/src/cfg.rs
@@ -277,7 +277,7 @@ impl Config {
                 interp_order: default_interp(),
                 code_smoothing: default_smoothing(),
                 min_sv_sunlight_rate: None,
-                min_sv_elev: Some(10.0),
+                min_sv_elev: Some(15.0),
                 min_snr: Some(30.0),
                 modeling: Modeling::default(),
                 max_sv: default_max_sv(),
@@ -314,10 +314,10 @@ impl Config {
 
 #[derive(Default, Debug, Clone, Copy, PartialEq)]
 #[cfg_attr(feature = "serde", derive(Deserialize))]
-pub enum SolverMode {
-    /// Receiver is kept at fixed location
+pub enum PositioningMode {
+    /// Receiver is static
     #[default]
     Static,
-    /// Receiver is not static
+    /// Receiver is moving
     Kinematic,
 }

From 58efddddfd4908925f0ff8954241b41bb95a3968 Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Sun, 3 Dec 2023 13:29:27 +0100
Subject: [PATCH 13/16] fixed issue in filter

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/solver.rs | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/solver.rs b/src/solver.rs
index 8fef2e2..ec8d9dd 100644
--- a/src/solver.rs
+++ b/src/solver.rs
@@ -357,7 +357,7 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
                 match mode {
                     Mode::SPP | Mode::LSQSPP => {
                         if pool[idx - nb_removed].code.len() == 0 {
-                            debug!("{:?} ({}) dropped on bad snr", pool[idx].t, pool[idx].sv);
+                            debug!("{:?} ({}) dropped on bad snr", pool[idx - nb_removed].t, pool[idx - nb_removed].sv);
                             let _ = pool.swap_remove(idx - nb_removed);
                             nb_removed += 1;
                         }
@@ -376,9 +376,10 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
 
         /* apply eclipse filter (if need be) */
         if let Some(min_rate) = self.cfg.min_sv_sunlight_rate {
-            for idx in 0..pool.len() - 1 {
-                let state = pool[idx].state.unwrap(); // infaillible
-                let orbit = state.orbit(pool[idx].t, self.earth_frame);
+            let mut nb_removed :usize = 0;
+            for idx in 0..pool.len() {
+                let state = pool[idx - nb_removed].state.unwrap(); // infaillible
+                let orbit = state.orbit(pool[idx - nb_removed].t, self.earth_frame);
                 let state = eclipse_state(&orbit, self.sun_frame, self.earth_frame, &self.cosmic);
                 let eclipsed = match state {
                     EclipseState::Umbra => true,
@@ -388,9 +389,10 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
                 if eclipsed {
                     debug!(
                         "{:?} ({}): dropped - eclipsed by earth",
-                        pool[idx].t, pool[idx].sv
+                        pool[idx - nb_removed].t, pool[idx - nb_removed].sv
                     );
-                    let _ = pool.swap_remove(idx);
+                    let _ = pool.swap_remove(idx - nb_removed);
+                    nb_removed += 1;
                 }
             }
         }

From bc53b0f5578f56d489810ea580a37dc71a3d54dc Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Sun, 3 Dec 2023 15:38:46 +0100
Subject: [PATCH 14/16] working on default cfg

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/cfg.rs | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/src/cfg.rs b/src/cfg.rs
index 8ee34e8..481c63f 100644
--- a/src/cfg.rs
+++ b/src/cfg.rs
@@ -84,17 +84,18 @@ fn default_relativistic_path_range() -> bool {
 }
 
 fn default_lsq_weight() -> Option<LSQWeight> {
-    Some(LSQWeight::LSQWeightMappingFunction(
-        ElevationMappingFunction {
-            a: 5.0,
-            b: 0.0,
-            c: 10.0,
-        },
-    ))
+    //Some(LSQWeight::LSQWeightMappingFunction(
+    //    ElevationMappingFunction {
+    //        a: 5.0,
+    //        b: 0.0,
+    //        c: 10.0,
+    //    },
+    //))
+    None
 }
 
 fn default_gdop_threshold() -> Option<f64> {
-    Some(0.15)
+    Some(1.0)
 }
 
 fn default_innov_threshold() -> Option<f64> {

From b3be4254fbdb0eee4328e98f19623d3e1af6feba Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Sun, 3 Dec 2023 19:30:49 +0100
Subject: [PATCH 15/16] apc correction

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 src/candidate.rs           |   4 +-
 src/cfg.rs                 |   7 +++
 src/lib.rs                 |   2 +-
 src/solutions/validator.rs |   2 +-
 src/solver.rs              | 124 +++++++++++++++++++++++++++----------
 5 files changed, 103 insertions(+), 36 deletions(-)

diff --git a/src/candidate.rs b/src/candidate.rs
index c44ffaf..7c1330d 100644
--- a/src/candidate.rs
+++ b/src/candidate.rs
@@ -4,7 +4,6 @@ use gnss::prelude::SV;
 use hifitime::Unit;
 use itertools::Itertools;
 use log::debug;
-use map_3d::deg2rad;
 
 use nyx::cosmic::SPEED_OF_LIGHT;
 use nyx::linalg::{DVector, MatrixXx4};
@@ -240,7 +239,8 @@ impl Candidate {
         let (x0, y0, z0) = apriori;
         let clock_corr = self.clock_corr.to_seconds();
         let (azimuth, elevation) = (state.azimuth, state.elevation);
-        let (sv_x, sv_y, sv_z) = state.position();
+        let sv_p = state.position();
+        let (sv_x, sv_y, sv_z) = (sv_p[0], sv_p[1], sv_p[2]);
 
         let mut sv_data = PVTSVData::default();
         sv_data.azimuth = azimuth;
diff --git a/src/cfg.rs b/src/cfg.rs
index 481c63f..98d9b36 100644
--- a/src/cfg.rs
+++ b/src/cfg.rs
@@ -83,6 +83,10 @@ fn default_relativistic_path_range() -> bool {
     false
 }
 
+fn default_sv_apc() -> bool {
+    false
+}
+
 fn default_lsq_weight() -> Option<LSQWeight> {
     //Some(LSQWeight::LSQWeightMappingFunction(
     //    ElevationMappingFunction {
@@ -161,6 +165,8 @@ pub struct Modeling {
     #[cfg_attr(feature = "serde", serde(default))]
     pub sv_total_group_delay: bool,
     #[cfg_attr(feature = "serde", serde(default))]
+    pub sv_apc: bool,
+    #[cfg_attr(feature = "serde", serde(default))]
     pub relativistic_clock_bias: bool,
     #[cfg_attr(feature = "serde", serde(default))]
     pub relativistic_path_range: bool,
@@ -176,6 +182,7 @@ impl Default for Modeling {
     fn default() -> Self {
         Self {
             sv_clock_bias: default_sv_clock(),
+            sv_apc: default_sv_apc(),
             iono_delay: default_iono(),
             tropo_delay: default_tropo(),
             sv_total_group_delay: default_sv_tgd(),
diff --git a/src/lib.rs b/src/lib.rs
index e329649..688e832 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -22,7 +22,7 @@ pub mod prelude {
     pub use crate::candidate::{Candidate, Observation};
     pub use crate::cfg::Config;
     pub use crate::solutions::{PVTSolution, PVTSolutionType};
-    pub use crate::solver::{InterpolationResult, Mode, Solver};
+    pub use crate::solver::{InterpolatedPosition, InterpolationResult, Mode, Solver};
     // re-export
     pub use gnss::prelude::{Constellation, SV};
     pub use hifitime::{Duration, Epoch, TimeScale};
diff --git a/src/solutions/validator.rs b/src/solutions/validator.rs
index 5cb933f..0cd17b2 100644
--- a/src/solutions/validator.rs
+++ b/src/solutions/validator.rs
@@ -42,7 +42,7 @@ impl SolutionValidator {
                 .unwrap();
 
             let pr = cd.prefered_pseudorange().unwrap().value;
-            let state = cd.state.unwrap().position;
+            let state = cd.state.unwrap().position();
             let estimate = apriori_ecef + solution.p;
 
             let (sv_x, sv_y, sv_z) = (state[0], state[1], state[2]);
diff --git a/src/solver.rs b/src/solver.rs
index ec8d9dd..fada121 100644
--- a/src/solver.rs
+++ b/src/solver.rs
@@ -14,7 +14,7 @@ use nyx::md::prelude::{Bodies, Frame, LightTimeCalc};
 
 use gnss::prelude::SV;
 
-use nalgebra::{DMatrix, DVector, Matrix3, MatrixXx4, Vector3};
+use nalgebra::{DVector, Matrix3, MatrixXx4, Vector3};
 
 use crate::{
     apriori::AprioriPosition,
@@ -98,12 +98,27 @@ pub enum Error {
     InvalidatedSolution(SolutionInvalidation),
 }
 
+/// Position interpolator helper
+#[derive(Copy, Clone, Debug, PartialEq)]
+pub enum InterpolatedPosition {
+    /// Providing Mass Center position
+    MassCenter(Vector3<f64>),
+    /// Providing APC position
+    AntennaPhaseCenter(Vector3<f64>),
+}
+
+impl Default for InterpolatedPosition {
+    fn default() -> Self {
+        Self::AntennaPhaseCenter(Vector3::<f64>::default())
+    }
+}
+
 /// Interpolation result (state vector) that needs to be
 /// resolved for every single candidate.
 #[derive(Copy, Clone, Debug, Default, PartialEq)]
 pub struct InterpolationResult {
     /// Position vector in [m] ECEF
-    pub position: Vector3<f64>,
+    pub position: InterpolatedPosition,
     /// Velocity vector in [m/s] ECEF
     pub velocity: Option<Vector3<f64>>,
     /// Elevation compared to reference position and horizon in [°]
@@ -113,27 +128,55 @@ pub struct InterpolationResult {
 }
 
 impl InterpolationResult {
-    pub(crate) fn position(&self) -> (f64, f64, f64) {
-        (self.position[0], self.position[1], self.position[2])
+    pub(crate) fn position(&self) -> Vector3<f64> {
+        match self.position {
+            InterpolatedPosition::MassCenter(mc) => mc,
+            InterpolatedPosition::AntennaPhaseCenter(pc) => pc,
+        }
     }
-    pub(crate) fn velocity(&self) -> Option<(f64, f64, f64)> {
-        let velocity = self.velocity?;
-        Some((velocity[0], velocity[1], velocity[2]))
+    pub(crate) fn velocity(&self) -> Option<Vector3<f64>> {
+        self.velocity
     }
     pub(crate) fn orbit(&self, dt: Epoch, frame: Frame) -> Orbit {
-        let (x, y, z) = self.position();
-        let (v_x, v_y, v_z) = self.velocity().unwrap_or((0.0_f64, 0.0_f64, 0.0_f64));
+        let p = self.position();
+        let v = self.velocity().unwrap_or(Vector3::<f64>::default());
         Orbit::cartesian(
-            x / 1000.0,
-            y / 1000.0,
-            z / 1000.0,
-            v_x / 1000.0,
-            v_y / 1000.0,
-            v_z / 1000.0,
+            p[0] / 1000.0,
+            p[1] / 1000.0,
+            p[2] / 1000.0,
+            v[0] / 1000.0,
+            v[1] / 1000.0,
+            v[2] / 1000.0,
             dt,
             frame,
         )
     }
+    /*
+     * Applies the APC conversion, if need be
+     */
+    pub(crate) fn mc_apc_correction(&mut self, r_sun: Vector3<f64>, delta_apc: f64) {
+        match self.position {
+            InterpolatedPosition::MassCenter(r_sat) => {
+                let k = -r_sat / (r_sat[0].powi(2) + r_sat[1].powi(2) + r_sat[2].powi(2)).sqrt();
+                let norm = ((r_sun[0] - r_sat[0]).powi(2)
+                    + (r_sun[1] - r_sat[1]).powi(2)
+                    + (r_sun[2] - r_sat[2]).powi(2))
+                .sqrt();
+                let e = (r_sun - r_sat) / norm;
+                let j = Vector3::<f64>::new(k[0] * e[0], k[1] * e[1], k[2] * e[2]);
+                let i = Vector3::<f64>::new(j[0] * k[0], j[1] * k[1], j[2] * k[2]);
+                let r = Matrix3::<f64>::new(i[0], j[0], k[0], i[1], j[1], k[1], i[2], j[2], k[2]);
+                let r_dot = Vector3::<f64>::new(
+                    (i[0] + j[0] + k[0]) * delta_apc,
+                    (i[1] + j[1] + k[1]) * delta_apc,
+                    (i[2] + j[2] + k[2]) * delta_apc,
+                );
+                let r_apc = r_sat + r_dot;
+                self.position = InterpolatedPosition::AntennaPhaseCenter(r_apc);
+            },
+            _ => {},
+        }
+    }
 }
 
 /// PVT Solver
@@ -155,9 +198,17 @@ where
     pub interpolator: I,
     /* Cosmic model */
     cosmic: Arc<Cosm>,
-    /* Earth frame */
+    /*
+     * (Reference) Earth frame.
+     * Could be relevant to rename this the day we want to
+     * resolve solutions on other Planets...  ; )
+     */
     earth_frame: Frame,
-    /* Sun frame */
+    /*
+     * Sun Body frame
+     * Could be relevant to rename this the day we want to resolve
+     * solutions in other Star systems....   o___O
+     */
     sun_frame: Frame,
     /* prev. solution for internal logic */
     prev_pvt: Option<(Epoch, PVTSolution)>,
@@ -238,6 +289,9 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         let solver_opts = &self.cfg.solver;
         let interp_order = self.cfg.interp_order;
 
+        let cosmic = &self.cosmic;
+        let earth_frame = self.earth_frame;
+
         /* interpolate positions */
         let mut pool: Vec<Candidate> = pool
             .iter()
@@ -245,6 +299,11 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
                 let (t_tx, dt_ttx) = c.transmission_time(&self.cfg).ok()?;
 
                 if let Some(mut interpolated) = (self.interpolator)(t_tx, c.sv, interp_order) {
+                    if modeling.sv_apc {
+                        let r_sun = Self::sun_unit_vector(&earth_frame, cosmic, t_tx);
+                        interpolated.mc_apc_correction(r_sun, 0.0_f64);
+                    }
+
                     let mut c = c.clone();
 
                     let rot = match modeling.earth_rotation {
@@ -263,7 +322,8 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
                         ),
                     };
 
-                    interpolated.position = rot * interpolated.position;
+                    interpolated.position =
+                        InterpolatedPosition::AntennaPhaseCenter(rot * interpolated.position());
 
                     if modeling.relativistic_clock_bias {
                         if interpolated.velocity.is_none() {
@@ -274,7 +334,7 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
                             if let Some((p_ttx, p_position)) = self.prev_sv_state.get(&c.sv) {
                                 let dt = (t_tx - *p_ttx).to_seconds();
                                 interpolated.velocity =
-                                    Some((rot * interpolated.position - rot * p_position) / dt);
+                                    Some((rot * interpolated.position() - rot * p_position) / dt);
                             }
                         }
 
@@ -296,7 +356,7 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
                         }
 
                         self.prev_sv_state
-                            .insert(c.sv, (t_tx, interpolated.position));
+                            .insert(c.sv, (t_tx, interpolated.position()));
                     }
 
                     debug!(
@@ -357,7 +417,11 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
                 match mode {
                     Mode::SPP | Mode::LSQSPP => {
                         if pool[idx - nb_removed].code.len() == 0 {
-                            debug!("{:?} ({}) dropped on bad snr", pool[idx - nb_removed].t, pool[idx - nb_removed].sv);
+                            debug!(
+                                "{:?} ({}) dropped on bad snr",
+                                pool[idx - nb_removed].t,
+                                pool[idx - nb_removed].sv
+                            );
                             let _ = pool.swap_remove(idx - nb_removed);
                             nb_removed += 1;
                         }
@@ -376,7 +440,7 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
 
         /* apply eclipse filter (if need be) */
         if let Some(min_rate) = self.cfg.min_sv_sunlight_rate {
-            let mut nb_removed :usize = 0;
+            let mut nb_removed: usize = 0;
             for idx in 0..pool.len() {
                 let state = pool[idx - nb_removed].state.unwrap(); // infaillible
                 let orbit = state.orbit(pool[idx - nb_removed].t, self.earth_frame);
@@ -389,7 +453,8 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
                 if eclipsed {
                     debug!(
                         "{:?} ({}): dropped - eclipsed by earth",
-                        pool[idx - nb_removed].t, pool[idx - nb_removed].sv
+                        pool[idx - nb_removed].t,
+                        pool[idx - nb_removed].sv
                     );
                     let _ = pool.swap_remove(idx - nb_removed);
                     nb_removed += 1;
@@ -480,17 +545,12 @@ impl<I: std::ops::Fn(Epoch, SV, usize) -> Option<InterpolationResult>> Solver<I>
         Ok((t, pvt_solution))
     }
     /*
-     * Evaluates Sun/Earth vector, <!> expressed in Km <!>
-     * for all SV NAV Epochs in provided context
+     * Evaluates Sun/Earth vector in meter ECEF at given Epoch
      */
-    fn sun_earth_vector(&mut self, t: Epoch) -> Vector3<f64> {
+    fn sun_unit_vector(ref_frame: &Frame, cosmic: &Arc<Cosm>, t: Epoch) -> Vector3<f64> {
         let sun_body = Bodies::Sun;
-        let orbit = self.cosmic.celestial_state(
-            sun_body.ephem_path(),
-            t,
-            self.earth_frame,
-            LightTimeCalc::None,
-        );
+        let orbit =
+            cosmic.celestial_state(sun_body.ephem_path(), t, *ref_frame, LightTimeCalc::None);
         Vector3::new(
             orbit.x_km * 1000.0,
             orbit.y_km * 1000.0,

From 9efda7149da342dc975f70d5ff7c42c4ea6d9ad0 Mon Sep 17 00:00:00 2001
From: "Guillaume W. Bres" <guillaume.bressaix@gmail.com>
Date: Sun, 3 Dec 2023 20:09:16 +0100
Subject: [PATCH 16/16] v0.4.0

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>
---
 Cargo.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Cargo.toml b/Cargo.toml
index 5ca1f96..6990c02 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "gnss-rtk"
-version = "0.3.1"
+version = "0.4.0"
 license = "MIT OR Apache-2.0"
 authors = ["Guillaume W. Bres <guillaume.bressaix@gmail.com>"]
 description = "GNSS position solver"