Develop (#1)

* sv relativistic effects compensation
* relativistic effect compensation
* update doc
* working on residual validation
* lib update
* make sure we still have enough candidates
* test signal against mode compliance
* run linter
---------

Signed-off-by: Guillaume W. Bres <guillaume.bressaix@gmail.com>

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"
@@ -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"] }

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.

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.

doc/lsqspp.md

@@ -0,0 +1,4 @@
+LSQSPP
+======
+
+TODO

doc/ppp.md

@@ -0,0 +1,4 @@
+PPP
+===
+
+TODO

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)
+

src/candidate.rs

@@ -2,9 +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};
@@ -93,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
      */
@@ -153,26 +179,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);
@@ -241,7 +253,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