Repository: paarnes/GNSS_Multipath_Analysis_Software (language: python)

Issue:
Create run-tests.yml



Discussion / hints:

Fix patch:
diff --git a/.devcontainer/Dockerfile b/.devcontainer/Dockerfile
new file mode 100644
index 0000000..3ab5ba3
--- /dev/null
+++ b/.devcontainer/Dockerfile
@@ -0,0 +1,23 @@
+# [Choice] Python version (use -bullseye variants on local arm64/Apple Silicon): 3, 3.10, 3.9, 3.8, 3.7, 3.6, 3-bullseye, 3.10-bullseye, 3.9-bullseye, 3.8-bullseye, 3.7-bullseye, 3.6-bullseye, 3-buster, 3.10-buster, 3.9-buster, 3.8-buster, 3.7-buster, 3.6-buster
+ARG VARIANT="3.10-bullseye"
+FROM mcr.microsoft.com/vscode/devcontainers/python:0-${VARIANT}
+
+# [Choice] Node.js version: none, lts/*, 16, 14, 12, 10
+ARG NODE_VERSION="none"
+RUN if [ "${NODE_VERSION}" != "none" ]; then su vscode -c "umask 0002 && . /usr/local/share/nvm/nvm.sh && nvm install ${NODE_VERSION} 2>&1"; fi
+
+# [Optional] If your pip requirements rarely change, uncomment this section to add them to the image.
+ COPY .devcontainer/requirements.txt /tmp/pip-tmp/
+  RUN pip3 --disable-pip-version-check --no-cache-dir install -r /tmp/pip-tmp/requirements.txt \
+    && rm -rf /tmp/pip-tmp
+
+
+
+# [Optional] Uncomment this section to install additional OS packages.
+RUN apt-get update && export DEBIAN_FRONTEND=noninteractive \
+     && apt-get -y install --no-install-recommends libpq-dev gdal-bin libgdal-dev  
+     # <your-package-list-here> 
+
+# [Optional] Uncomment this line to install global node packages.
+# RUN su vscode -c "source /usr/local/share/nvm/nvm.sh && npm install -g <your-package-here>" 2>&1
+
diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json
new file mode 100644
index 0000000..0e85335
--- /dev/null
+++ b/.devcontainer/devcontainer.json
@@ -0,0 +1,55 @@
+// For format details, see https://aka.ms/devcontainer.json. For config options, see the README at:
+// https://github.com/microsoft/vscode-dev-containers/tree/v0.238.1/containers/python-3
+{
+	"name": "Python 3",
+	"build": {
+		"dockerfile": "Dockerfile",
+		"context": "..",
+		"args": { 
+			// Update 'VARIANT' to pick a Python version: 3, 3.10, 3.9, 3.8, 3.7, 3.6
+			// Append -bullseye or -buster to pin to an OS version.
+			// Use -bullseye variants on local on arm64/Apple Silicon.
+			"VARIANT": "3.10-bullseye",
+			// Options
+			"NODE_VERSION": "lts/*"
+		}
+	},
+	
+	
+	// Configure tool-specific properties.
+	"customizations": {
+		// Configure properties specific to VS Code.
+		"vscode": {
+			// Set *default* container specific settings.json values on container create.
+			"settings": { 
+				"python.defaultInterpreterPath": "/usr/local/bin/python",
+				"python.linting.enabled": true,
+				"python.linting.pylintEnabled": true,
+				"python.formatting.autopep8Path": "/usr/local/py-utils/bin/autopep8",
+				"python.formatting.blackPath": "/usr/local/py-utils/bin/black",
+				"python.formatting.yapfPath": "/usr/local/py-utils/bin/yapf",
+				"python.linting.banditPath": "/usr/local/py-utils/bin/bandit",
+				"python.linting.flake8Path": "/usr/local/py-utils/bin/flake8",
+				"python.linting.mypyPath": "/usr/local/py-utils/bin/mypy",
+				"python.linting.pycodestylePath": "/usr/local/py-utils/bin/pycodestyle",
+				"python.linting.pydocstylePath": "/usr/local/py-utils/bin/pydocstyle",
+				"python.linting.pylintPath": "/usr/local/py-utils/bin/pylint"
+			},
+			
+			// Add the IDs of extensions you want installed when the container is created.
+			"extensions": [
+				"ms-python.python",
+				"ms-python.vscode-pylance"
+			]
+		}
+	},
+
+	// Use 'forwardPorts' to make a list of ports inside the container available locally.
+	// "forwardPorts": [],
+
+	// Use 'postCreateCommand' to run commands after the container is created.
+	// "postCreateCommand": "pip3 install --user -r requirements.txt",
+
+	// Comment out to connect as root instead. More info: https://aka.ms/vscode-remote/containers/non-root.
+	"remoteUser": "vscode"
+}
diff --git a/.devcontainer/requirements.txt b/.devcontainer/requirements.txt
new file mode 100644
index 0000000..558f978
--- /dev/null
+++ b/.devcontainer/requirements.txt
@@ -0,0 +1,4 @@
+matplotlib
+numpy
+pandas
+tqdm
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 1508cd2..93d3514 100644
--- a/.gitignore
+++ b/.gitignore
@@ -129,3 +129,4 @@ dmypy.json
 .pyre/
 
 Admin 
+launch.json
\ No newline at end of file
diff --git a/Multipath_analysis/Running_multipath_analysis_example.py b/Multipath_analysis/Running_multipath_analysis_example.py
deleted file mode 100644
index 70a12dd..0000000
--- a/Multipath_analysis/Running_multipath_analysis_example.py
+++ /dev/null
@@ -1,124 +0,0 @@
-import os, pickle,numpy as np
-from GNSS_MultipathAnalysis import GNSS_MultipathAnalysis
-from readRinexObs import readRinexObs
-
-abs_path = os.path.abspath("GNSS_MultipathAnalysis.py") #find path to the main file
-base_path = os.path.normpath(os.getcwd() + os.sep + os.pardir)
-## Path to TestData
-relpath_to_testdata = 'TestData'
-full_path_testdata = os.path.join(base_path, relpath_to_testdata) 
-
-#Path to outputdir 
-relpath_to_outputdir = 'Results'
-full_path_ouputdir = os.path.join(base_path, relpath_to_outputdir) 
-
-## -----  Defining input data --------
-## Rinex observation file
-rinObsFilename1 = full_path_testdata  + '/ObservationFiles/' + 'OPEC00NOR_S_20220010000_01D_30S_MO_3.04.rnx'
-rinObsFilename2 = full_path_testdata  + '/ObservationFiles/' + 'NMBUS_SAMSUNG_S20.20o'
-
-## SP3 files
-sp3NavFilename_1      = full_path_testdata  + '/SP3/' + 'NMBUS_2020 10 30.SP3'
-sp3NavFilename_1_opec = full_path_testdata  + '/SP3/' + 'Testfile_20220101.eph'
-
-
-## Broadcast ephemerides
-broadcastNav1 =  full_path_testdata  + '/NavigationFiles/' + 'OPEC00NOR_S_20220010000_01D_GN.rnx'
-broadcastNav2 =  full_path_testdata  + '/NavigationFiles/' + 'OPEC00NOR_S_20220010000_01D_RN.rnx'
-broadcastNav3 =  full_path_testdata  + '/NavigationFiles/' + 'BCEmerge_30_10_2020.20p'
-broadcastNav4 =  full_path_testdata  + '/NavigationFiles/' + 'BRDC00IGS_R_20220010000_01D_MN.rnx'
-
-
-
-#%% #####-- Simple examples for running analysis (few userdefined settings) 
-# Using SP3  OPEC
-
-analysisResults = GNSS_MultipathAnalysis(rinObsFilename1,
-                                         desiredGNSSsystems=["G"],
-                                         sp3NavFilename_1 = sp3NavFilename_1_opec)
-
-
-#%% Using broadcasted eph OPEC
-analysisResults = GNSS_MultipathAnalysis(rinObsFilename1,
-                                         desiredGNSSsystems=['G'], # Run analysis on GPS only ['G','R','E'] for GPS,GLONASS and Galileo 
-                                         outputDir = full_path_ouputdir,
-                                         broadcastNav1=broadcastNav4
-                                         )
-
-#%% Using broadcasted eph NMBUS
-analysisResults = GNSS_MultipathAnalysis(rinObsFilename2,
-                                          desiredGNSSsystems=['E'],
-                                          outputDir = full_path_ouputdir,
-                                          broadcastNav1=broadcastNav3
-                                          )
-#%% Using SP3  NMBUS
-analysisResults = GNSS_MultipathAnalysis(rinObsFilename2,
-                                         desiredGNSSsystems=['G'],
-                                         outputDir = full_path_ouputdir,
-                                         sp3NavFilename_1=sp3NavFilename_1,
-                                         plotEstimates=True
-                                         )
-
-#%% Using SP3  NMBUS without plotting
-analysisResults = GNSS_MultipathAnalysis(rinObsFilename2,
-                                         desiredGNSSsystems=['G'],
-                                         outputDir = full_path_ouputdir,
-                                         broadcastNav1=broadcastNav3,
-                                         plotEstimates=False,
-                                         plot_polarplot=False
-                                         )
-
-#%% -- Advanced example (more user defined settings)
-
-## Parameters
-GNSSsystems                 = ["R"] # run analysis in GLONASS only
-phaseCodeLimit              = 6.667
-ionLimit                    = 0.0667
-cutoff_elevation_angle      = 10 # 10 degree elevation cutoff
-outputDir                   = full_path_ouputdir
-plotEstimates               = True
-plot_polarplot              = True
-includeResultSummary        = True
-includeCompactSummary       = True
-includeObservationOverview  = True
-includeLLIOverview          = True
-
-
-
-analysisResults = GNSS_MultipathAnalysis(rinObsFilename1,
-                       desiredGNSSsystems=GNSSsystems,
-                       sp3NavFilename_1 = sp3NavFilename_1_opec,
-                       phaseCodeLimit = phaseCodeLimit,
-                       ionLimit = ionLimit,
-                       cutoff_elevation_angle = cutoff_elevation_angle,
-                       outputDir = outputDir, 
-                       plotEstimates = plotEstimates,
-                       plot_polarplot=plot_polarplot,
-                       includeResultSummary = includeResultSummary, 
-                       includeCompactSummary = includeCompactSummary, 
-                       includeObservationOverview = includeObservationOverview,
-                       includeLLIOverview = includeLLIOverview
-                       )
- 
-  
-
-#%% ------------ How to read in the result file from a analysis
-path_to_resFile = os.path.join(full_path_ouputdir, 'analysisResults.pkl')
-file_to_read = open(path_to_resFile, "rb")
-loaded_dictionary = pickle.load(file_to_read)
-
-
-#%% ----------- Examplecode for using readrinex to read in GPS and Galielo observation
-rinObsFilename = r'C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\TestData\ObservationFiles\OPEC00NOR_S_20220010000_01D_30S_MO_3.04.rnx'
-desiredGNSSsystems = ["G","E"]          # GPS and Galileo
-desiredObsCodes = ["C", "L"]       # code and phase observations
-desiredObsBands = list(np.arange(1,10)) # all carrier bands. 
-
-## --- Read RINEX 3.0x observation file
-[GNSS_obs, GNSS_LLI, GNSS_SS, GNSS_SVs, time_epochs, nepochs, GNSSsystems,\
-    obsCodes, approxPosition, max_sat, tInterval, markerName, rinexVersion, recType, timeSystem, leapSec, gnssType,\
-    rinexProgr, rinexDate, antDelta, tFirstObs, tLastObs, clockOffsetsON, GLO_Slot2ChannelMap, success] = \
-    readRinexObs(rinObsFilename,  
-                 desiredGNSSsystems = desiredGNSSsystems,
-                 desiredObsCodes = desiredObsCodes,
-                 desiredObsBands = desiredObsBands)
diff --git a/README.md b/README.md
index ed3e438..bfcc9b7 100644
--- a/README.md
+++ b/README.md
@@ -22,7 +22,7 @@ To install the required packages, run:
 `pip install -r requirements.txt`   
 where the *requirements.txt* is located [here](https://github.com/paarnes/GNSS/blob/master/Multipath_analysis/requirements.txt).
 
-Note: In the example plots, TEX is used to get prettier text formatting. However, this requires TEX/LaTex to be installed on your computer. The program will first try to use TEX, and if it's not possible, standard text formatting will be used. So TEX/LaTex is not required to run the program and make plots anymore. The software will just not use TEX as a interpreter. In the previous release, the software required TEX to be installed to successfully create plots.
+Note: In the example plots, TEX is used to get prettier text formatting. However, this requires TEX/LaTex to be installed on your computer. The program will first try to use TEX, and if it's not possible, standard text formatting will be used. So TEX/LaTex is not required to run the program and make plots.
 
 ## The steps are:
 1. Reads in the RINEX observation file 
diff --git a/src/Examples_on_how_to_run_it.ipynb b/src/Examples_on_how_to_run_it.ipynb
new file mode 100644
index 0000000..2a5b0be
--- /dev/null
+++ b/src/Examples_on_how_to_run_it.ipynb
@@ -0,0 +1,149 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Example on how to run the software\n",
+    "\n",
+    "Run the first cell first to import modules and set paths\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "from GNSS_MultipathAnalysis import GNSS_MultipathAnalysis\n",
+    "\n",
+    "base_path = os.getcwd()\n",
+    "parent_dir = os.path.abspath(os.path.join(base_path, os.pardir))\n",
+    "\n",
+    "## Path to TestData\n",
+    "path_to_testdata = os.path.join(parent_dir, 'TestData')\n",
+    "ouputdir = os.path.join(parent_dir, 'Results')\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Define the path to the RINEX observation file and the SP3 or navigation file\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rinObsFilename1 = os.path.join(path_to_testdata, 'ObservationFiles', 'OPEC00NOR_S_20220010000_01D_30S_MO_3.04.rnx')\n",
+    "broadcastNav4 = os.path.join(path_to_testdata, 'NavigationFiles', 'BRDC00IGS_R_20220010000_01D_MN.rnx')\n",
+    "analysisResults = GNSS_MultipathAnalysis(rinObsFilename1, broadcastNav1=broadcastNav4)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "####  Advanced example (more user defined settings)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Parameters\n",
+    "GNSSsystems                 = [\"R\"] # run analysis in GLONASS only\n",
+    "phaseCodeLimit              = 6.667\n",
+    "ionLimit                    = 0.0667\n",
+    "cutoff_elevation_angle      = 10 # 10 degree elevation cutoff\n",
+    "outputDir                   = ouputdir\n",
+    "plotEstimates               = True\n",
+    "plot_polarplot              = True\n",
+    "includeResultSummary        = True\n",
+    "includeCompactSummary       = True\n",
+    "includeObservationOverview  = True\n",
+    "includeLLIOverview          = True\n",
+    "\n",
+    "\n",
+    "## Rinex observation file\n",
+    "rinObsFilename1 = path_to_testdata  + '/ObservationFiles/' + 'OPEC00NOR_S_20220010000_01D_30S_MO_3.04.rnx'\n",
+    "## SP3 files\n",
+    "sp3NavFilename_1_opec = path_to_testdata  + '/SP3/' + 'Testfile_20220101.eph'\n",
+    "\n",
+    "analysisResults = GNSS_MultipathAnalysis(rinObsFilename1,\n",
+    "                       desiredGNSSsystems=GNSSsystems,\n",
+    "                       sp3NavFilename_1 = sp3NavFilename_1_opec,\n",
+    "                       phaseCodeLimit = phaseCodeLimit,\n",
+    "                       ionLimit = ionLimit,\n",
+    "                       cutoff_elevation_angle = cutoff_elevation_angle,\n",
+    "                       outputDir = outputDir, \n",
+    "                       plotEstimates = plotEstimates,\n",
+    "                       plot_polarplot=plot_polarplot,\n",
+    "                       includeResultSummary = includeResultSummary, \n",
+    "                       includeCompactSummary = includeCompactSummary, \n",
+    "                       includeObservationOverview = includeObservationOverview,\n",
+    "                       includeLLIOverview = includeLLIOverview\n",
+    "                       )"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "How to read in the result file from a analysis"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#%% ------------ How to read in the result file from a analysis\n",
+    "import pickle\n",
+    "path_to_resFile = os.path.join(ouputdir, 'analysisResults.pkl')\n",
+    "file_to_read = open(path_to_resFile, \"rb\")\n",
+    "loaded_dictionary = pickle.load(file_to_read)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "base",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.7"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "477348ab2256e79b25884e91eb7bf3223e9c8ae3ede6ff068233d5704eb1b027"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/Multipath_analysis/GNSS_MultipathAnalysis.py b/src/GNSS_MultipathAnalysis.py
similarity index 98%
rename from Multipath_analysis/GNSS_MultipathAnalysis.py
rename to src/GNSS_MultipathAnalysis.py
index ecaf909..b4dedc1 100644
--- a/Multipath_analysis/GNSS_MultipathAnalysis.py
+++ b/src/GNSS_MultipathAnalysis.py
@@ -1,4 +1,6 @@
-import os, numpy as np, pickle
+import os
+import numpy as np
+import pickle
 from readRinexObs import *
 from Geodetic_functions import *
 from computeSatElevations import computeSatElevations
@@ -199,33 +201,33 @@ def GNSS_MultipathAnalysis(rinObsFilename,
     
     ## ---  Control of the user input arguments 
     if type(sp3NavFilename_1) != str:
-        print('ERROR(GNSS_Receiver_QC_2020): The input variable sp3NavFilename_1 must be a string\n' \
+        print('ERROR(GNSS_MultipathAnalysis): The input variable sp3NavFilename_1 must be a string\n' \
             'Argument is now of type %s\n' %  (type(sp3NavFilename_1)))
         analysisResults = np.nan
         return
     
     if type(sp3NavFilename_2) != str:
-        print('ERROR(GNSS_Receiver_QC_2020): The input variable sp3NavFilename_2 must be a string\n' \
+        print('ERROR(GNSS_MultipathAnalysis): The input variable sp3NavFilename_2 must be a string\n' \
             'Argument is now of type %s\n' %  (type(sp3NavFilename_2)))
         analysisResults = np.nan
         return
     
     
     if type(sp3NavFilename_3) != str:
-        print('ERROR(GNSS_Receiver_QC_2020): The input variable sp3NavFilename_3must be a string\n' \
+        print('ERROR(GNSS_MultipathAnalysis): The input variable sp3NavFilename_3must be a string\n' \
             'Argument is now of type %s\n' %  (type(sp3NavFilename_3)))
         analysisResults = np.nan
         return
     
     if type(rinObsFilename) != str:
-        print('ERROR(GNSS_Receiver_QC_2020): The input variable rinObsFilename must be a string\n' \
+        print('ERROR(GNSS_MultipathAnalysis): The input variable rinObsFilename must be a string\n' \
             'Argument is now of type %s\n' %  (type(rinObsFilename)))
         analysisResults = np.nan
         return
  
     
     if not os.path.isfile(rinObsFilename):
-        print('ERROR(GNSS_Receiver_QC_2020): RINEX observation file can not be found. Please check that the path is correct.\n') 
+        print('ERROR(GNSS_MultipathAnalysis): RINEX observation file can not be found. Please check that the path is correct.\n') 
         analysisResults = np.nan
         return 
     
diff --git a/Multipath_analysis/Geodetic_functions.py b/src/Geodetic_functions.py
similarity index 100%
rename from Multipath_analysis/Geodetic_functions.py
rename to src/Geodetic_functions.py
diff --git a/Multipath_analysis/Rinex_Frequency_Overview.txt b/src/Rinex_Frequency_Overview.txt
similarity index 100%
rename from Multipath_analysis/Rinex_Frequency_Overview.txt
rename to src/Rinex_Frequency_Overview.txt
diff --git a/src/Tests/analysisResults_NMBUS.pkl b/src/Tests/analysisResults_NMBUS.pkl
new file mode 100644
index 0000000..de92237
Binary files /dev/null and b/src/Tests/analysisResults_NMBUS.pkl differ
diff --git a/src/Tests/analysisResults_OPEC00NOR_S_20220010000_01D_30S_MO.pkl b/src/Tests/analysisResults_OPEC00NOR_S_20220010000_01D_30S_MO.pkl
new file mode 100644
index 0000000..3f58c21
Binary files /dev/null and b/src/Tests/analysisResults_OPEC00NOR_S_20220010000_01D_30S_MO.pkl differ
diff --git a/Multipath_analysis/barylag.py b/src/barylag.py
similarity index 100%
rename from Multipath_analysis/barylag.py
rename to src/barylag.py
diff --git a/Multipath_analysis/computeDelayStats.py b/src/computeDelayStats.py
similarity index 100%
rename from Multipath_analysis/computeDelayStats.py
rename to src/computeDelayStats.py
diff --git a/Multipath_analysis/computeSatElevAzimuth_fromNav.py b/src/computeSatElevAzimuth_fromNav.py
similarity index 100%
rename from Multipath_analysis/computeSatElevAzimuth_fromNav.py
rename to src/computeSatElevAzimuth_fromNav.py
diff --git a/Multipath_analysis/computeSatElevations.py b/src/computeSatElevations.py
similarity index 98%
rename from Multipath_analysis/computeSatElevations.py
rename to src/computeSatElevations.py
index 38c5d69..076aae2 100644
--- a/Multipath_analysis/computeSatElevations.py
+++ b/src/computeSatElevations.py
@@ -255,6 +255,8 @@ def get_elevation_angle(sys, PRN, week, tow, sat_positions, nEpochs, epoch_dates
         ## -- Calculate elevation and azimut angle from receiver to satellite
         if not np.isnan(u) and not np.isnan(e) and not np.isnan(n):
             elevation_angle = np.rad2deg(atanc(u, sqrt(e**2 + n**2)))
+            if not 0 < elevation_angle < 90: # if the satellite is below the horizon (elevation angle is below zero)
+                elevation_angle = np.nan
             
             # Azimut computation and quadrant correction 
             if (e> 0 and n< 0) or (e < 0 and n < 0):
diff --git a/Multipath_analysis/detectClockJumps.py b/src/detectClockJumps.py
similarity index 100%
rename from Multipath_analysis/detectClockJumps.py
rename to src/detectClockJumps.py
diff --git a/Multipath_analysis/detectCycleSlips.py b/src/detectCycleSlips.py
similarity index 100%
rename from Multipath_analysis/detectCycleSlips.py
rename to src/detectCycleSlips.py
diff --git a/Multipath_analysis/estimateSignalDelays.py b/src/estimateSignalDelays.py
similarity index 98%
rename from Multipath_analysis/estimateSignalDelays.py
rename to src/estimateSignalDelays.py
index 118de03..9e5da6a 100644
--- a/Multipath_analysis/estimateSignalDelays.py
+++ b/src/estimateSignalDelays.py
@@ -1,343 +1,343 @@
-import numpy as np
-from detectCycleSlips import detectCycleSlips, orgSlipEpochs
-import warnings
-warnings.filterwarnings(action='ignore', message='Mean of empty slice')
-
-def estimateSignalDelays(range1_Code, range2_Code,phase1_Code, phase2_Code, carrier_freq1, \
-                         carrier_freq2, nepochs, max_sat, GNSS_SVs, obsCodes, GNSS_obs, \
-                             currentGNSSsystem, tInterval, phaseCodeLimit, ionLimit):
-    """
-     Function that takes observations of from the observation period and
-     estimates the following delays on the signal:
-                                                   ionospheric delay on the
-                                                   first phase signal
-    
-                                                   multipath delay on the
-                                                   first range signal
-    
-                                                   multipath delay on the
-                                                   second range signal
-    
-     Also estimates epochs with ambiguity slip for the first signal.
-    
-     The estimates are all relative estimates. The multipath estimates are
-     relative to a mean values. The ionosphere delay estimates are relative to
-     the first estimate. These estimates must be reduced again by the
-     releative value at ambiguity slips. The function calls on another
-     function to estimate periods of ambiguity slips. It then corrects the
-     delay estimates for these slips.
-    --------------------------------------------------------------------------------------------------------------------------
-     INPUTS
-    
-     range1_Code:          string, defines the observation type that will be 
-                           used as the first range observation. ex. "C1X", "C5X"    
-    
-     range2_Code:          string, defines the observation type that will be 
-                           used as the second range observation. ex. "C1X", "C5X"
-    
-     phase1_Code:          string, defines the observation type that will be 
-                           used as the first phase observation. ex. "L1X", "L5X"
-    
-     phase2_Code:          string, defines the observation type that will be 
-                           used as the second phase observation. ex. "L1X", "L5X"
-    
-     carrier_freq1:        carrier frequency of first phase signal. unit Hz
-    
-     carrier_freq2:        carrier frequency of second phase signal. unit Hz
-    
-     nepochs:              number of epochs with observations in rinex observation file.
-    
-     max_sat:              max PRN number of current GNSS system.
-    
-     GNSS_SVs:             matrix containing number of satellites with 
-                           obsevations for each epoch, and PRN for those satellites
-    
-                           GNSS_SVs(epoch, j)  j=1: number of observed satellites
-                                               j>1: PRN of observed satellites
-    
-     obsCodes:             cell contaning strings. Cell defines the observation
-                           codes available from the current GNSS system. Each
-                           string is a three-character string, the first 
-                           character (a capital letter) is an observation code 
-                           ex. "L" or "C". The second character (a digit)
-                           is a frequency code. The third character(a Capital letter)  
-                           is the attribute, ex. "P" or "X"
-    
-     GNSS_obs:             3D matrix containing all observation of current 
-                           GNSS system for all epochs. Order of obsType index
-                           is same order as in obsCodes cell
-    
-                           GNSS_obs(PRN, obsType, epoch)
-                                               PRN: double
-                                               ObsType: double: 1,2,...,numObsTypes
-                                               epoch: double
-    
-     currentGNSSsystem:    string, code to indicate which GNSS system is the
-                           current observations are coming from. Only used to
-                           decide if system uses CDMA or FDMA.
-    
-     tInterval:            observations interval seconds. 
-    
-     phaseCodeLimit:       critical limit that indicates cycle slip for
-                           phase-code combination. Unit: m/s. If set to 0,
-                           default value will be used
-    
-     ionLimit:             critical limit that indicates cycle slip for
-                           the rate of change of the ionopheric delay. 
-                           Unit: m/s. If set to 0, default value will be used
-    --------------------------------------------------------------------------------------------------------------------------
-     OUTPUTS
-    
-     ion_delay_phase1:     matrix containing estimates of ionospheric delays 
-                           on the first phase signal for each PRN, at each epoch.
-    
-                           ion_delay_phase1(epoch, PRN)
-    
-     multipath_range1:     matrix containing estimates of multipath delays 
-                           on the first range signal for each PRN, at each epoch.
-    
-                           multipath_range2(epoch, PRN)
-    
-     multipath_range2:     matrix containing estimates of multipath delays 
-                           on the second range signal for each PRN, at each epoch.
-    
-                           multipath_range2(epoch, PRN)
-    
-     range1_slip_periods:  cell, each cell element contains one matrix for
-                           every PRN. Each matrix contains epochs of range1
-                           slip starts in first column and slip ends in second
-                           column.
-    
-     range1_observations:  matrix. Contains all range1 observations for all
-                           epochs and for all SVs. 
-    
-                           range1_observations(epoch, PRN)
-    
-     phase1_observations:  matrix. Contains all phase1 observations for all
-                           epochs and for all SVs. 
-    
-                           range1_observations(epoch, PRN)
-    
-     success:              boolean, 1 if no error occurs, 0 otherwise
-    
-    --------------------------------------------------------------------------------------------------------------------------
-    """
-    
-    FDMA_used = 0
-    success = 1
-    
-    if 'R' in currentGNSSsystem:
-        FDMA_used = 1
-        carrier_freq1_list = carrier_freq1
-        carrier_freq2_list = carrier_freq2
-    else:
-        alpha = carrier_freq1**2/carrier_freq2**2 # amplfication factor
-    
-    # Define parameters
-    c = 299792458 # speed of light
-    
-    if ionLimit ==0:
-        ionLimit = 4/60   # critical rate of change of ionosphere delay  to indicate ambiguity slip on either 
-                          # range1/phase1 signal, range2/phase2 signal, or both
-    
-    
-    if phaseCodeLimit == 0:
-        phaseCodeLimit = 4/60*100   # critical rate of change of 
-                                    # N1_pseudo_estimate to indicate ambiguity slip on
-                                    # the range1/phase1 signal
- 
-    ## -- Initialize data matrices
-    ion_delay_phase1        = np.zeros([nepochs, max_sat+1]) # Addin 1 to since python is nullindexed
-    multipath_range1        = np.zeros([nepochs, max_sat+1])
-    # multipath_range2        = np.zeros([nepochs, max_sat+1]) # kommenterte bort 23.01.2023
-    N1_pseudo_estimate      = np.zeros([nepochs, max_sat+1])
-    missing_obs_overview    = np.zeros([nepochs, max_sat+1])
-    missing_range1_overview = np.zeros([nepochs, max_sat+1])
-    
-    ## -- Initialize cell for storing phase slip periods
-    ambiguity_slip_periods = {}
-    ## -- Initialize cell for storing slip periods for only range1/phase1
-
-    range1_slip_periods = {}
-    
-    ## -- Make estimates of delays
-    for epoch in np.arange(0,nepochs):
-        ## -- Calculate how many satelites with observations in current epoch
-        n_sat = int(GNSS_SVs[epoch,0])
-        ## - Iterate over all satellites of current epoch
-        for i in np.arange(0,n_sat):  
-            PRN = int(GNSS_SVs[epoch, 1+i]) # Get PRN ## KANSJE FJERNE +1??? for å få med PRN 11
-            
-            if FDMA_used:
-                carrier_freq1 = carrier_freq1_list[PRN]
-                carrier_freq2 = carrier_freq2_list[PRN]
-                alpha = carrier_freq1**2/carrier_freq2**2 # amplfication factor
-            
-            
-            ## -- Get observations            
-            range1 = GNSS_obs[epoch+1][PRN, ismember(obsCodes[currentGNSSsystem],range1_Code)] 
-            range2 = GNSS_obs[epoch+1][PRN, ismember(obsCodes[currentGNSSsystem],range2_Code)]
-            
-            phase1 = GNSS_obs[epoch+1][PRN, ismember(obsCodes[currentGNSSsystem],phase1_Code)]*c/carrier_freq1
-            phase2 = GNSS_obs[epoch+1][PRN, ismember(obsCodes[currentGNSSsystem],phase2_Code)]*c/carrier_freq2
-            
-            if any([str(range1), str(range2), str(phase1), str(phase2)]) == '[]':
-                print('ERROR(estimateSignalDelays): There is no observation type #s. Check for missing data in RINEX observation file!' % (phase1_Code))
-                success = 0
-                return success 
-            
-            
-            ## -- If any of the four observations are missing, ie 0, estimate remains 0 for that epoch and satellite
-            if [range1, range2, phase1, phase2].count(0) == 0:
-                ## -- Calculate estimate for Ionospheric delay on phase 1 signal
-                ion_delay_phase1[epoch, PRN] = 1/(alpha-1)*(phase1-phase2)
-                
-                ## -- Calculate multipath on both code signals
-                multipath_range1[epoch,PRN] = range1 - (1 + 2/(alpha-1))*phase1 + (2/(alpha-1))*phase2
-                # multipath_range2[epoch,PRN] = range2 - (2*alpha/(alpha-1))*phase1 + ((2*alpha)/(alpha-1) - 1)*phase2 ## kommenterte bort 23.01.2023
-                N1_pseudo_estimate[epoch, PRN] = phase1 - range1
-            else:
-                ## Flag epoch and PRN as missing obs
-                missing_obs_overview[epoch, PRN] = 1
-                if range1 == 0 or phase1 == 0:
-                    missing_range1_overview[epoch, PRN] = 1 
-                else:
-                    N1_pseudo_estimate[epoch, PRN] = phase1 - range1
-
-    
-    ## -- Detect and correct for ambiguity slips
-    for PRN in np.arange(0,max_sat):
-        PRN = PRN + 1
-        ## -- Get first and last epoch with observations for current PRN
-        if len(np.nonzero(ion_delay_phase1[:,PRN])[0]) != 0:
-            epoch_first_obs = np.nonzero(ion_delay_phase1[:,PRN])[0][0]
-            epoch_last_obs = np.nonzero(ion_delay_phase1[:,PRN])[0][-1]
-        else:
-            epoch_first_obs = np.nonzero(ion_delay_phase1[:,PRN])[0]
-            epoch_last_obs = np.nonzero(ion_delay_phase1[:,PRN])[0]
-                
-               
-        ## -- Get first and last epoch with range 1 observations for current PRN
-        if len(np.nonzero(N1_pseudo_estimate[:,PRN])[0]) != 0:
-            epoch_first_range1obs =  np.nonzero(N1_pseudo_estimate[:,PRN])[0][0]
-            epoch_last_range1obs =  np.nonzero(N1_pseudo_estimate[:,PRN])[0][-1]
-        else:
-            epoch_first_range1obs =  np.nonzero(N1_pseudo_estimate[:,PRN])[0]
-            epoch_last_range1obs =  np.nonzero(N1_pseudo_estimate[:,PRN])[0]
-        
-        
-        ## -- Run function to detect cycle slips for current epoch for range1/phase1 only 
-        range1_slip_epochs = detectCycleSlips(N1_pseudo_estimate[:, PRN], missing_range1_overview[:, PRN],epoch_first_range1obs, epoch_last_range1obs, tInterval, phaseCodeLimit)
-       
-        # Run function to detect cycle slips for current epoch for either range1/phase1 signal, range2/phase2 signal, or both
-        ionosphere_slip_epochs = detectCycleSlips(ion_delay_phase1[:, PRN], missing_obs_overview[:, PRN],epoch_first_obs, epoch_last_obs, tInterval, ionLimit)
-           
-        ## -- Make combined array of slip epochs from both lin. combinations used to detects slips
-        ambiguity_slip_epochs = np.union1d(range1_slip_epochs, ionosphere_slip_epochs) 
-        # range1_slip_epochs = np.intersect1d(range1_slip_epochs, ionosphere_slip_epochs) # hvorfor intersect her?
-        range1_slip_epochs = range1_slip_epochs #tester om det gir stor forskjell uten intersect
-
-        ## -- Organize slips detected on range1/phase1 signal only
-        range1_slip_periods[PRN],_ = orgSlipEpochs(range1_slip_epochs)
-        
-        ## -- Orginize combined slips detected on range1/phase1 signal only
-        ambiguity_slip_periods[PRN], n_slip_periods = orgSlipEpochs(ambiguity_slip_epochs)
-        
-        # Set all zero estimates to NaN so that epochs with missing observations are not corrected        
-        ion_delay_phase1[ion_delay_phase1[:, PRN]==0, PRN] = np.nan
-        multipath_range1[multipath_range1[:, PRN]==0, PRN] = np.nan
-        # multipath_range2[multipath_range2[:, PRN]==0, PRN] = np.nan ## 23.01.2023 kommenterer bort
-        
-        ## If there are no slips then there is only one "ambiguity period". All estimates are therefore reduced by the same relative value
-        if len(ambiguity_slip_periods[PRN]) == 0:
-            # if len(list(epoch_first_obs)) == 0:
-            if epoch_first_obs.size == 0: # added 18.02.2023 because of error when running on RINEX v2 (should be like this either way!)
-                pass
-            else:    
-                ion_delay_phase1[epoch_first_obs::, PRN] = ion_delay_phase1[epoch_first_obs::, PRN] - ion_delay_phase1[epoch_first_obs, PRN]  
-                # multipath_range1[epoch_first_obs::, PRN] = multipath_range1[epoch_first_obs::, PRN] - np.mean(np.nonzero(multipath_range1[epoch_first_obs::, PRN]))
-                ### removed np.nonzero and use np.nanmean instead since every zero is replaced by nan in line 247 and 248
-                multipath_range1[epoch_first_obs::, PRN] = multipath_range1[epoch_first_obs::, PRN] - np.nanmean(multipath_range1[epoch_first_obs::, PRN])  
-
-                # multipath_range2[epoch_first_obs::, PRN] = multipath_range2[epoch_first_obs::, PRN] - np.mean(np.nonzero(multipath_range2[epoch_first_obs::, PRN])) 
-        else:
-            ## -- Set all estimates of epochs with cycle slips to nan
-            for slip_period in np.arange(0,n_slip_periods):
-                slip_start   = int(ambiguity_slip_periods[PRN][slip_period,0])
-                slip_end     = int(ambiguity_slip_periods[PRN][slip_period,1])
-               
-                if slip_start == slip_end:  # need a if test bacause if there equal, python dont set to nan
-                    ion_delay_phase1[slip_start, PRN] = np.nan
-                    multipath_range1[slip_start, PRN] = np.nan
-                    # multipath_range2[slip_start, PRN] = np.nan
-                else: 
-                    ion_delay_phase1[slip_start:slip_end+1, PRN] = np.nan # + 1 because a[2:3] gives one element. Matlab a(2:3) gives 2 element.
-                    multipath_range1[slip_start:slip_end+1, PRN] = np.nan # if error msg here, try add "if slip_end != epoch_last_obs else epoch_last_obs" in the slicing (oneliner)
-                    # multipath_range2[slip_start:slip_end+1, PRN] = np.nan
-            
-              
-            ## Extract start and end of each segment and correct multipath and ionosphere estimates for each segment
-            for ambiguity_period in np.arange(0,n_slip_periods+1): # removed + 1 cause of indexproblem 29.11
-                if ambiguity_period == 0: 
-                    ambiguity_period_start  = epoch_first_obs
-                    ambiguity_period_end    = int(ambiguity_slip_periods[PRN][0,0])   #INK -1 igjen???
-                ## -- If last ambiguity period
-                elif ambiguity_period == n_slip_periods: # removed + 1 
-                    ambiguity_period_start       = int(ambiguity_slip_periods[PRN][-1,1] + 1) 
-                    ambiguity_period_end         = epoch_last_obs +1
-                   
-                    ## If last epoch with observation is a slip, then there is no  last ambiguity period
-                    if ambiguity_period_start > epoch_last_obs:
-                        ambiguity_period_start = []
-                        ambiguity_period_end = []
-                    
-                else:
-                    ambiguity_period_start = int(ambiguity_slip_periods[PRN][ambiguity_period-1, 1] + 1)
-                    ambiguity_period_end   = int(ambiguity_slip_periods[PRN][ambiguity_period, 0])
-                
-               
-                ## -- Ionosphere delay estimates of current ambiguity period is reduced by first estimate of ambiguity period
-                if ambiguity_period_start != ambiguity_period_end:
-                    ion_delay_phase1[ambiguity_period_start:ambiguity_period_end, PRN] = ion_delay_phase1[ambiguity_period_start:ambiguity_period_end, PRN] - \
-                        ion_delay_phase1[ambiguity_period_start, PRN]
-                   
-                    ## -- Multipath delays of current ambiguity period are reduced by mean of estimates in ambiguity period, excluding NaN and
-                    multipath_range1[ambiguity_period_start:ambiguity_period_end, PRN] = multipath_range1[ambiguity_period_start:ambiguity_period_end, PRN] - \
-                        np.nanmean(multipath_range1[ambiguity_period_start:ambiguity_period_end, PRN]) # added nanmean 30.11
-                        
-                    # multipath_range2[ambiguity_period_start:ambiguity_period_end, PRN] = multipath_range2[ambiguity_period_start:ambiguity_period_end, PRN] -\
-                        # np.nanmean(multipath_range2[ambiguity_period_start:ambiguity_period_end, PRN])
-                else:
-                    
-                    ion_delay_phase1[ambiguity_period_start, PRN] = ion_delay_phase1[ambiguity_period_start, PRN] - \
-                        ion_delay_phase1[ambiguity_period_start, PRN]
-                    ## -- Multipath delays of current ambiguity period are reduced by mean of estimates in ambiguity period, excluding NaN and
-                    multipath_range1[ambiguity_period_start, PRN] = multipath_range1[ambiguity_period_start, PRN] - \
-                        np.nanmean(multipath_range1[ambiguity_period_start, PRN])
-                       
-                    # multipath_range2[ambiguity_period_start, PRN] = multipath_range2[ambiguity_period_start, PRN] -\
-                        # np.nanmean(multipath_range2[ambiguity_period_start, PRN])
-                        
-        ## -- Get range1 and phase 1 observations for all epochs and PRN
-        range1_observations =  np.zeros([nepochs, max_sat+1]) 
-        phase1_observations =  np.zeros([nepochs, max_sat+1]) 
-        for ep in np.arange(0, len(GNSS_obs)):
-            for PRN in np.arange(0,max_sat):
-                range1_observations[ep,PRN] = GNSS_obs[ep+1][PRN,ismember(obsCodes[currentGNSSsystem],range1_Code)] 
-                phase1_observations[ep,PRN] = GNSS_obs[ep+1][PRN,ismember(obsCodes[currentGNSSsystem],phase1_Code)]
-
-    # return ion_delay_phase1, multipath_range1, multipath_range2, range1_slip_periods, range1_observations, phase1_observations, success
-    # return ion_delay_phase1, multipath_range1, multipath_range2, ambiguity_slip_periods, range1_observations, phase1_observations, success # changeing from range1slip to amgiguity
-    return ion_delay_phase1, multipath_range1, range1_slip_periods,ambiguity_slip_periods, range1_observations, phase1_observations, success #removed multipath_range2
-
-
-def ismember(list_,code):
-    """
-    The function takes in a string and a list, and finds the index of 
-    """
-    indx = [idx for idx, val in enumerate(list_) if val == code]
-    if indx != []:
-        indx = indx[0]
-    return indx
-
-
+import numpy as np
+from detectCycleSlips import detectCycleSlips, orgSlipEpochs
+import warnings
+warnings.filterwarnings(action='ignore', message='Mean of empty slice')
+
+def estimateSignalDelays(range1_Code, range2_Code,phase1_Code, phase2_Code, carrier_freq1, \
+                         carrier_freq2, nepochs, max_sat, GNSS_SVs, obsCodes, GNSS_obs, \
+                             currentGNSSsystem, tInterval, phaseCodeLimit, ionLimit):
+    """
+     Function that takes observations of from the observation period and
+     estimates the following delays on the signal:
+                                                   ionospheric delay on the
+                                                   first phase signal
+    
+                                                   multipath delay on the
+                                                   first range signal
+    
+                                                   multipath delay on the
+                                                   second range signal
+    
+     Also estimates epochs with ambiguity slip for the first signal.
+    
+     The estimates are all relative estimates. The multipath estimates are
+     relative to a mean values. The ionosphere delay estimates are relative to
+     the first estimate. These estimates must be reduced again by the
+     releative value at ambiguity slips. The function calls on another
+     function to estimate periods of ambiguity slips. It then corrects the
+     delay estimates for these slips.
+    --------------------------------------------------------------------------------------------------------------------------
+     INPUTS
+    
+     range1_Code:          string, defines the observation type that will be 
+                           used as the first range observation. ex. "C1X", "C5X"    
+    
+     range2_Code:          string, defines the observation type that will be 
+                           used as the second range observation. ex. "C1X", "C5X"
+    
+     phase1_Code:          string, defines the observation type that will be 
+                           used as the first phase observation. ex. "L1X", "L5X"
+    
+     phase2_Code:          string, defines the observation type that will be 
+                           used as the second phase observation. ex. "L1X", "L5X"
+    
+     carrier_freq1:        carrier frequency of first phase signal. unit Hz
+    
+     carrier_freq2:        carrier frequency of second phase signal. unit Hz
+    
+     nepochs:              number of epochs with observations in rinex observation file.
+    
+     max_sat:              max PRN number of current GNSS system.
+    
+     GNSS_SVs:             matrix containing number of satellites with 
+                           obsevations for each epoch, and PRN for those satellites
+    
+                           GNSS_SVs(epoch, j)  j=1: number of observed satellites
+                                               j>1: PRN of observed satellites
+    
+     obsCodes:             cell contaning strings. Cell defines the observation
+                           codes available from the current GNSS system. Each
+                           string is a three-character string, the first 
+                           character (a capital letter) is an observation code 
+                           ex. "L" or "C". The second character (a digit)
+                           is a frequency code. The third character(a Capital letter)  
+                           is the attribute, ex. "P" or "X"
+    
+     GNSS_obs:             3D matrix containing all observation of current 
+                           GNSS system for all epochs. Order of obsType index
+                           is same order as in obsCodes cell
+    
+                           GNSS_obs(PRN, obsType, epoch)
+                                               PRN: double
+                                               ObsType: double: 1,2,...,numObsTypes
+                                               epoch: double
+    
+     currentGNSSsystem:    string, code to indicate which GNSS system is the
+                           current observations are coming from. Only used to
+                           decide if system uses CDMA or FDMA.
+    
+     tInterval:            observations interval seconds. 
+    
+     phaseCodeLimit:       critical limit that indicates cycle slip for
+                           phase-code combination. Unit: m/s. If set to 0,
+                           default value will be used
+    
+     ionLimit:             critical limit that indicates cycle slip for
+                           the rate of change of the ionopheric delay. 
+                           Unit: m/s. If set to 0, default value will be used
+    --------------------------------------------------------------------------------------------------------------------------
+     OUTPUTS
+    
+     ion_delay_phase1:     matrix containing estimates of ionospheric delays 
+                           on the first phase signal for each PRN, at each epoch.
+    
+                           ion_delay_phase1(epoch, PRN)
+    
+     multipath_range1:     matrix containing estimates of multipath delays 
+                           on the first range signal for each PRN, at each epoch.
+    
+                           multipath_range2(epoch, PRN)
+    
+     multipath_range2:     matrix containing estimates of multipath delays 
+                           on the second range signal for each PRN, at each epoch.
+    
+                           multipath_range2(epoch, PRN)
+    
+     range1_slip_periods:  cell, each cell element contains one matrix for
+                           every PRN. Each matrix contains epochs of range1
+                           slip starts in first column and slip ends in second
+                           column.
+    
+     range1_observations:  matrix. Contains all range1 observations for all
+                           epochs and for all SVs. 
+    
+                           range1_observations(epoch, PRN)
+    
+     phase1_observations:  matrix. Contains all phase1 observations for all
+                           epochs and for all SVs. 
+    
+                           range1_observations(epoch, PRN)
+    
+     success:              boolean, 1 if no error occurs, 0 otherwise
+    
+    --------------------------------------------------------------------------------------------------------------------------
+    """
+    
+    FDMA_used = 0
+    success = 1
+    
+    if 'R' in currentGNSSsystem:
+        FDMA_used = 1
+        carrier_freq1_list = carrier_freq1
+        carrier_freq2_list = carrier_freq2
+    else:
+        alpha = carrier_freq1**2/carrier_freq2**2 # amplfication factor
+    
+    # Define parameters
+    c = 299792458 # speed of light
+    
+    if ionLimit ==0:
+        ionLimit = 4/60   # critical rate of change of ionosphere delay  to indicate ambiguity slip on either 
+                          # range1/phase1 signal, range2/phase2 signal, or both
+    
+    
+    if phaseCodeLimit == 0:
+        phaseCodeLimit = 4/60*100   # critical rate of change of 
+                                    # N1_pseudo_estimate to indicate ambiguity slip on
+                                    # the range1/phase1 signal
+ 
+    ## -- Initialize data matrices
+    ion_delay_phase1        = np.zeros([nepochs, max_sat+1]) # Addin 1 to since python is nullindexed
+    multipath_range1        = np.zeros([nepochs, max_sat+1])
+    # multipath_range2        = np.zeros([nepochs, max_sat+1]) # kommenterte bort 23.01.2023
+    N1_pseudo_estimate      = np.zeros([nepochs, max_sat+1])
+    missing_obs_overview    = np.zeros([nepochs, max_sat+1])
+    missing_range1_overview = np.zeros([nepochs, max_sat+1])
+    
+    ## -- Initialize cell for storing phase slip periods
+    ambiguity_slip_periods = {}
+    ## -- Initialize cell for storing slip periods for only range1/phase1
+
+    range1_slip_periods = {}
+    
+    ## -- Make estimates of delays
+    for epoch in np.arange(0,nepochs):
+        ## -- Calculate how many satelites with observations in current epoch
+        n_sat = int(GNSS_SVs[epoch,0])
+        ## - Iterate over all satellites of current epoch
+        for i in np.arange(0,n_sat):  
+            PRN = int(GNSS_SVs[epoch, 1+i]) # Get PRN ## KANSJE FJERNE +1??? for å få med PRN 11
+            
+            if FDMA_used:
+                carrier_freq1 = carrier_freq1_list[PRN]
+                carrier_freq2 = carrier_freq2_list[PRN]
+                alpha = carrier_freq1**2/carrier_freq2**2 # amplfication factor
+            
+            
+            ## -- Get observations            
+            range1 = GNSS_obs[epoch+1][PRN, ismember(obsCodes[currentGNSSsystem],range1_Code)] 
+            range2 = GNSS_obs[epoch+1][PRN, ismember(obsCodes[currentGNSSsystem],range2_Code)]
+            
+            phase1 = GNSS_obs[epoch+1][PRN, ismember(obsCodes[currentGNSSsystem],phase1_Code)]*c/carrier_freq1
+            phase2 = GNSS_obs[epoch+1][PRN, ismember(obsCodes[currentGNSSsystem],phase2_Code)]*c/carrier_freq2
+            
+            if any([str(range1), str(range2), str(phase1), str(phase2)]) == '[]':
+                print('ERROR(estimateSignalDelays): There is no observation type #s. Check for missing data in RINEX observation file!' % (phase1_Code))
+                success = 0
+                return success 
+            
+            
+            ## -- If any of the four observations are missing, ie 0, estimate remains 0 for that epoch and satellite
+            if [range1, range2, phase1, phase2].count(0) == 0:
+                ## -- Calculate estimate for Ionospheric delay on phase 1 signal
+                ion_delay_phase1[epoch, PRN] = 1/(alpha-1)*(phase1-phase2)
+                
+                ## -- Calculate multipath on both code signals
+                multipath_range1[epoch,PRN] = range1 - (1 + 2/(alpha-1))*phase1 + (2/(alpha-1))*phase2
+                # multipath_range2[epoch,PRN] = range2 - (2*alpha/(alpha-1))*phase1 + ((2*alpha)/(alpha-1) - 1)*phase2 ## kommenterte bort 23.01.2023
+                N1_pseudo_estimate[epoch, PRN] = phase1 - range1
+            else:
+                ## Flag epoch and PRN as missing obs
+                missing_obs_overview[epoch, PRN] = 1
+                if range1 == 0 or phase1 == 0:
+                    missing_range1_overview[epoch, PRN] = 1 
+                else:
+                    N1_pseudo_estimate[epoch, PRN] = phase1 - range1
+
+    
+    ## -- Detect and correct for ambiguity slips
+    for PRN in np.arange(0,max_sat):
+        PRN = PRN + 1
+        ## -- Get first and last epoch with observations for current PRN
+        if len(np.nonzero(ion_delay_phase1[:,PRN])[0]) != 0:
+            epoch_first_obs = np.nonzero(ion_delay_phase1[:,PRN])[0][0]
+            epoch_last_obs = np.nonzero(ion_delay_phase1[:,PRN])[0][-1]
+        else:
+            epoch_first_obs = np.nonzero(ion_delay_phase1[:,PRN])[0]
+            epoch_last_obs = np.nonzero(ion_delay_phase1[:,PRN])[0]
+                
+               
+        ## -- Get first and last epoch with range 1 observations for current PRN
+        if len(np.nonzero(N1_pseudo_estimate[:,PRN])[0]) != 0:
+            epoch_first_range1obs =  np.nonzero(N1_pseudo_estimate[:,PRN])[0][0]
+            epoch_last_range1obs =  np.nonzero(N1_pseudo_estimate[:,PRN])[0][-1]
+        else:
+            epoch_first_range1obs =  np.nonzero(N1_pseudo_estimate[:,PRN])[0]
+            epoch_last_range1obs =  np.nonzero(N1_pseudo_estimate[:,PRN])[0]
+        
+        
+        ## -- Run function to detect cycle slips for current epoch for range1/phase1 only 
+        range1_slip_epochs = detectCycleSlips(N1_pseudo_estimate[:, PRN], missing_range1_overview[:, PRN],epoch_first_range1obs, epoch_last_range1obs, tInterval, phaseCodeLimit)
+       
+        # Run function to detect cycle slips for current epoch for either range1/phase1 signal, range2/phase2 signal, or both
+        ionosphere_slip_epochs = detectCycleSlips(ion_delay_phase1[:, PRN], missing_obs_overview[:, PRN],epoch_first_obs, epoch_last_obs, tInterval, ionLimit)
+           
+        ## -- Make combined array of slip epochs from both lin. combinations used to detects slips
+        ambiguity_slip_epochs = np.union1d(range1_slip_epochs, ionosphere_slip_epochs) 
+        # range1_slip_epochs = np.intersect1d(range1_slip_epochs, ionosphere_slip_epochs) # hvorfor intersect her?
+        range1_slip_epochs = range1_slip_epochs #tester om det gir stor forskjell uten intersect
+
+        ## -- Organize slips detected on range1/phase1 signal only
+        range1_slip_periods[PRN],_ = orgSlipEpochs(range1_slip_epochs)
+        
+        ## -- Orginize combined slips detected on range1/phase1 signal only
+        ambiguity_slip_periods[PRN], n_slip_periods = orgSlipEpochs(ambiguity_slip_epochs)
+        
+        # Set all zero estimates to NaN so that epochs with missing observations are not corrected        
+        ion_delay_phase1[ion_delay_phase1[:, PRN]==0, PRN] = np.nan
+        multipath_range1[multipath_range1[:, PRN]==0, PRN] = np.nan
+        # multipath_range2[multipath_range2[:, PRN]==0, PRN] = np.nan ## 23.01.2023 kommenterer bort
+        
+        ## If there are no slips then there is only one "ambiguity period". All estimates are therefore reduced by the same relative value
+        if len(ambiguity_slip_periods[PRN]) == 0:
+            # if len(list(epoch_first_obs)) == 0:
+            if epoch_first_obs.size == 0: # added 18.02.2023 because of error when running on RINEX v2 (should be like this either way!)
+                pass
+            else:    
+                ion_delay_phase1[epoch_first_obs::, PRN] = ion_delay_phase1[epoch_first_obs::, PRN] - ion_delay_phase1[epoch_first_obs, PRN]  
+                # multipath_range1[epoch_first_obs::, PRN] = multipath_range1[epoch_first_obs::, PRN] - np.mean(np.nonzero(multipath_range1[epoch_first_obs::, PRN]))
+                ### removed np.nonzero and use np.nanmean instead since every zero is replaced by nan in line 247 and 248
+                multipath_range1[epoch_first_obs::, PRN] = multipath_range1[epoch_first_obs::, PRN] - np.nanmean(multipath_range1[epoch_first_obs::, PRN])  
+
+                # multipath_range2[epoch_first_obs::, PRN] = multipath_range2[epoch_first_obs::, PRN] - np.mean(np.nonzero(multipath_range2[epoch_first_obs::, PRN])) 
+        else:
+            ## -- Set all estimates of epochs with cycle slips to nan
+            for slip_period in np.arange(0,n_slip_periods):
+                slip_start   = int(ambiguity_slip_periods[PRN][slip_period,0])
+                slip_end     = int(ambiguity_slip_periods[PRN][slip_period,1])
+               
+                if slip_start == slip_end:  # need a if test bacause if there equal, python dont set to nan
+                    ion_delay_phase1[slip_start, PRN] = np.nan
+                    multipath_range1[slip_start, PRN] = np.nan
+                    # multipath_range2[slip_start, PRN] = np.nan
+                else: 
+                    ion_delay_phase1[slip_start:slip_end+1, PRN] = np.nan # + 1 because a[2:3] gives one element. Matlab a(2:3) gives 2 element.
+                    multipath_range1[slip_start:slip_end+1, PRN] = np.nan # if error msg here, try add "if slip_end != epoch_last_obs else epoch_last_obs" in the slicing (oneliner)
+                    # multipath_range2[slip_start:slip_end+1, PRN] = np.nan
+            
+              
+            ## Extract start and end of each segment and correct multipath and ionosphere estimates for each segment
+            for ambiguity_period in np.arange(0,n_slip_periods+1): # removed + 1 cause of indexproblem 29.11
+                if ambiguity_period == 0: 
+                    ambiguity_period_start  = epoch_first_obs
+                    ambiguity_period_end    = int(ambiguity_slip_periods[PRN][0,0])   #INK -1 igjen???
+                ## -- If last ambiguity period
+                elif ambiguity_period == n_slip_periods: # removed + 1 
+                    ambiguity_period_start       = int(ambiguity_slip_periods[PRN][-1,1] + 1) 
+                    ambiguity_period_end         = epoch_last_obs +1
+                   
+                    ## If last epoch with observation is a slip, then there is no  last ambiguity period
+                    if ambiguity_period_start > epoch_last_obs:
+                        ambiguity_period_start = []
+                        ambiguity_period_end = []
+                    
+                else:
+                    ambiguity_period_start = int(ambiguity_slip_periods[PRN][ambiguity_period-1, 1] + 1)
+                    ambiguity_period_end   = int(ambiguity_slip_periods[PRN][ambiguity_period, 0])
+                
+               
+                ## -- Ionosphere delay estimates of current ambiguity period is reduced by first estimate of ambiguity period
+                if ambiguity_period_start != ambiguity_period_end:
+                    ion_delay_phase1[ambiguity_period_start:ambiguity_period_end, PRN] = ion_delay_phase1[ambiguity_period_start:ambiguity_period_end, PRN] - \
+                        ion_delay_phase1[ambiguity_period_start, PRN]
+                   
+                    ## -- Multipath delays of current ambiguity period are reduced by mean of estimates in ambiguity period, excluding NaN and
+                    multipath_range1[ambiguity_period_start:ambiguity_period_end, PRN] = multipath_range1[ambiguity_period_start:ambiguity_period_end, PRN] - \
+                        np.nanmean(multipath_range1[ambiguity_period_start:ambiguity_period_end, PRN]) # added nanmean 30.11
+                        
+                    # multipath_range2[ambiguity_period_start:ambiguity_period_end, PRN] = multipath_range2[ambiguity_period_start:ambiguity_period_end, PRN] -\
+                        # np.nanmean(multipath_range2[ambiguity_period_start:ambiguity_period_end, PRN])
+                else:
+                    
+                    ion_delay_phase1[ambiguity_period_start, PRN] = ion_delay_phase1[ambiguity_period_start, PRN] - \
+                        ion_delay_phase1[ambiguity_period_start, PRN]
+                    ## -- Multipath delays of current ambiguity period are reduced by mean of estimates in ambiguity period, excluding NaN and
+                    multipath_range1[ambiguity_period_start, PRN] = multipath_range1[ambiguity_period_start, PRN] - \
+                        np.nanmean(multipath_range1[ambiguity_period_start, PRN])
+                       
+                    # multipath_range2[ambiguity_period_start, PRN] = multipath_range2[ambiguity_period_start, PRN] -\
+                        # np.nanmean(multipath_range2[ambiguity_period_start, PRN])
+                        
+        ## -- Get range1 and phase 1 observations for all epochs and PRN
+        range1_observations =  np.zeros([nepochs, max_sat+1]) 
+        phase1_observations =  np.zeros([nepochs, max_sat+1]) 
+        for ep in np.arange(0, len(GNSS_obs)):
+            for PRN in np.arange(0,max_sat):
+                range1_observations[ep,PRN] = GNSS_obs[ep+1][PRN,ismember(obsCodes[currentGNSSsystem],range1_Code)] 
+                phase1_observations[ep,PRN] = GNSS_obs[ep+1][PRN,ismember(obsCodes[currentGNSSsystem],phase1_Code)]
+
+    # return ion_delay_phase1, multipath_range1, multipath_range2, range1_slip_periods, range1_observations, phase1_observations, success
+    # return ion_delay_phase1, multipath_range1, multipath_range2, ambiguity_slip_periods, range1_observations, phase1_observations, success # changeing from range1slip to amgiguity
+    return ion_delay_phase1, multipath_range1, range1_slip_periods,ambiguity_slip_periods, range1_observations, phase1_observations, success #removed multipath_range2
+
+
+def ismember(list_,code):
+    """
+    The function takes in a string and a list, and finds the index of 
+    """
+    indx = [idx for idx, val in enumerate(list_) if val == code]
+    if indx != []:
+        indx = indx[0]
+    return indx
+
+
diff --git a/Multipath_analysis/getLLISlipPeriods.py b/src/getLLISlipPeriods.py
similarity index 100%
rename from Multipath_analysis/getLLISlipPeriods.py
rename to src/getLLISlipPeriods.py
diff --git a/Multipath_analysis/gpstime2date.py b/src/gpstime2date.py
similarity index 100%
rename from Multipath_analysis/gpstime2date.py
rename to src/gpstime2date.py
diff --git a/Multipath_analysis/make_polarplot.py b/src/make_polarplot.py
similarity index 100%
rename from Multipath_analysis/make_polarplot.py
rename to src/make_polarplot.py
diff --git a/Multipath_analysis/make_polarplot_dont_use_TEX.py b/src/make_polarplot_dont_use_TEX.py
similarity index 100%
rename from Multipath_analysis/make_polarplot_dont_use_TEX.py
rename to src/make_polarplot_dont_use_TEX.py
diff --git a/Multipath_analysis/plotResults.py b/src/plotResults.py
similarity index 100%
rename from Multipath_analysis/plotResults.py
rename to src/plotResults.py
diff --git a/Multipath_analysis/preciseOrbits2ECEF.py b/src/preciseOrbits2ECEF.py
similarity index 100%
rename from Multipath_analysis/preciseOrbits2ECEF.py
rename to src/preciseOrbits2ECEF.py
diff --git a/Multipath_analysis/readFrequencyOverview.py b/src/readFrequencyOverview.py
similarity index 100%
rename from Multipath_analysis/readFrequencyOverview.py
rename to src/readFrequencyOverview.py
diff --git a/Multipath_analysis/readRinexNav.py b/src/readRinexNav.py
similarity index 100%
rename from Multipath_analysis/readRinexNav.py
rename to src/readRinexNav.py
diff --git a/Multipath_analysis/readRinexObs.py b/src/readRinexObs.py
similarity index 99%
rename from Multipath_analysis/readRinexObs.py
rename to src/readRinexObs.py
index 03ace1b..262596b 100644
--- a/Multipath_analysis/readRinexObs.py
+++ b/src/readRinexObs.py
@@ -3372,13 +3372,14 @@ def rinexReadObsBlockHead211(fid):
         print(msg2)
 
 
-    SVlist = []
-    sat_overview = line.split(" ")[-1][len(str(numSV))::]
-    pattern = re.compile(r'[A-Z][0-9]{2}')
-    sat_list = re.findall(pattern, sat_overview)
-    for sat in sat_list: ## droppe denne forloopen ?? SVlist og sat_list inneholder det samme??
-        SVlist.append(sat)
-    
+    # SVlist = []
+    # sat_overview = line.split(" ")[-1][len(str(numSV))::]
+    # pattern = re.compile(r'[A-Z][0-9]{2}')
+    # sat_list = re.findall(pattern, sat_overview)
+    # for sat in sat_list: ## droppe denne forloopen ?? SVlist og sat_list inneholder det samme??
+    #     SVlist.append(sat)    
+
+    SVlist = re.findall(r'[A-Z][0-9]{2}', line)   
 
     return success, epochflag, clockOffset, date, numSV, SVlist, eof
 
diff --git a/Multipath_analysis/readRinexObs2_11.py b/src/readRinexObs2_11.py
similarity index 100%
rename from Multipath_analysis/readRinexObs2_11.py
rename to src/readRinexObs2_11.py
diff --git a/Multipath_analysis/read_SP3Nav.py b/src/read_SP3Nav.py
similarity index 100%
rename from Multipath_analysis/read_SP3Nav.py
rename to src/read_SP3Nav.py
diff --git a/Multipath_analysis/requirements.txt b/src/requirements.txt
similarity index 100%
rename from Multipath_analysis/requirements.txt
rename to src/requirements.txt
diff --git a/Multipath_analysis/signalAnalysis.py b/src/signalAnalysis.py
similarity index 100%
rename from Multipath_analysis/signalAnalysis.py
rename to src/signalAnalysis.py
diff --git a/Multipath_analysis/writeOutputFile.py b/src/writeOutputFile.py
similarity index 100%
rename from Multipath_analysis/writeOutputFile.py
rename to src/writeOutputFile.py


Test patch:
diff --git a/.github/workflows/run-tests.yml b/.github/workflows/run-tests.yml
new file mode 100644
index 0000000..ad1a46a
--- /dev/null
+++ b/.github/workflows/run-tests.yml
@@ -0,0 +1,28 @@
+name: Python application
+
+on:
+  pull_request:
+    branches: [ "master","dev"]
+
+permissions:
+  contents: read
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v3
+    - name: Set up Python 3.10
+      uses: actions/setup-python@v3
+      with:
+        python-version: "3.10"
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install pytest
+        pip install -r .devcontainer/requirements.txt
+    - name: Test with pytest
+      run: |
+        pytest src/Tests/test_GNSS_MultipathAnalysis.py
diff --git a/src/Tests/test_GNSS_MultipathAnalysis.py b/src/Tests/test_GNSS_MultipathAnalysis.py
new file mode 100644
index 0000000..747c548
--- /dev/null
+++ b/src/Tests/test_GNSS_MultipathAnalysis.py
@@ -0,0 +1,105 @@
+import pytest
+import pickle
+import sys
+import os
+import numpy as np
+from numpy.testing import assert_almost_equal
+
+# os.chdir(r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\Multipath_analysis")
+sys.path.append(r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\src")
+from GNSS_MultipathAnalysis import GNSS_MultipathAnalysis
+os.chdir(r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\src")
+
+
+
+def read_pickle(file_path):
+    with open(file_path, 'rb') as file:
+        data = pickle.load(file)
+    return data
+
+
+
+def test_GNSS_MultipathAnalysis_sp3_file():
+    """
+    Test the results from OPEC2022 and sp3 files
+    """
+    rinObs_file =  r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\TestData\ObservationFiles\NMBUS_SAMSUNG_S20.20o"
+    sp3Nav_file = r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\TestData\SP3\NMBUS_2020 10 30.SP3"
+    expected_res = r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\Multipath_analysis\Tests\analysisResults_NMBUS.pkl"
+    result = GNSS_MultipathAnalysis(rinObsFilename=rinObs_file, sp3NavFilename_1=sp3Nav_file,
+                                    plotEstimates=False,
+                                    plot_polarplot=False)
+    
+    expected_result = read_pickle(expected_res)
+    # Compare the result with the expected result
+    assert_almost_equal(expected_result['Sat_position']['G']['Position'][2],  result['Sat_position']['G']['Position'][2],decimal=4)
+    assert_almost_equal(expected_result['Sat_position']['G']['Position'][4],  result['Sat_position']['G']['Position'][4],decimal=4)
+    assert_almost_equal(expected_result['Sat_position']['G']['Position'][30],  result['Sat_position']['G']['Position'][30],decimal=4)
+
+    assert_almost_equal(expected_result['Sat_position']['G']['Azimut'],  result['Sat_position']['G']['Azimut'],decimal=4)
+    assert_almost_equal(expected_result['Sat_position']['E']['Azimut'],  result['Sat_position']['E']['Azimut'],decimal=4)
+
+    assert_almost_equal(expected_result['Sat_position']['G']['Elevation'],  result['Sat_position']['G']['Elevation'],decimal=4)
+    assert_almost_equal(expected_result['Sat_position']['E']['Elevation'],  result['Sat_position']['E']['Elevation'],decimal=4)
+
+
+# def test_GNSS_MultipathAnalysis_broadcast_navfile():
+#     """
+#     Test the results from OPEC2022 and use of broadcast ephemerids
+#     """
+#     rinObs_file =  r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\TestData\ObservationFiles\OPEC00NOR_S_20220010000_01D_30S_MO_3.04.rnx"
+#     broadNav_file = r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\TestData\NavigationFiles\BRDC00IGS_R_20220010000_01D_MN.rnx"
+#     expected_res = r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\Multipath_analysis\Tests\analysisResults_OPEC00NOR_S_20220010000_01D_30S_MO.pkl"
+#     result = GNSS_MultipathAnalysis(rinObsFilename=rinObs_file, broadcastNav1=broadNav_file,
+#                                     plotEstimates=False,
+#                                     plot_polarplot=False)
+#     expected_result = read_pickle(expected_res)
+#     # Compare the result with the expected result
+#     assert_almost_equal(expected_result['Sat_position']['G']['Position'][1],  result['Sat_position']['G']['Position'][1],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['G']['Position'][8],  result['Sat_position']['G']['Position'][8],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['G']['Position'][10], result['Sat_position']['G']['Position'][10],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['R']['Position'][1],  result['Sat_position']['R']['Position'][1],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['R']['Position'][7],  result['Sat_position']['R']['Position'][7],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['R']['Position'][8],  result['Sat_position']['R']['Position'][8],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['E']['Position'][1],  result['Sat_position']['E']['Position'][1],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['E']['Position'][3],  result['Sat_position']['E']['Position'][3],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['E']['Position'][8],  result['Sat_position']['E']['Position'][8],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['C']['Position'][5],  result['Sat_position']['C']['Position'][5],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['C']['Position'][6],  result['Sat_position']['C']['Position'][6],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['C']['Position'][9],  result['Sat_position']['C']['Position'][9],decimal=4)
+
+#     assert_almost_equal(expected_result['Sat_position']['G']['Azimut'],  result['Sat_position']['G']['Azimut'],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['R']['Azimut'],  result['Sat_position']['R']['Azimut'],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['E']['Azimut'],  result['Sat_position']['E']['Azimut'],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['C']['Azimut'],  result['Sat_position']['C']['Azimut'],decimal=4)
+
+#     assert_almost_equal(expected_result['Sat_position']['G']['Elevation'],  result['Sat_position']['G']['Elevation'],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['R']['Elevation'],  result['Sat_position']['R']['Elevation'],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['E']['Elevation'],  result['Sat_position']['E']['Elevation'],decimal=4)
+#     assert_almost_equal(expected_result['Sat_position']['C']['Elevation'],  result['Sat_position']['C']['Elevation'],decimal=4)
+
+#     assert_almost_equal(expected_result['GPS']['Band_1']['C1C']['rms_multipath_range1_averaged'], result['GPS']['Band_1']['C1C']['rms_multipath_range1_averaged'], decimal=4)
+#     assert_almost_equal(expected_result['GPS']['Band_1']['C1C']['elevation_weighted_average_rms_multipath_range1'], result['GPS']['Band_1']['C1C']['elevation_weighted_average_rms_multipath_range1'], decimal=4)
+#     assert_almost_equal(expected_result['GLONASS']['Band_1']['C1C']['rms_multipath_range1_averaged'], result['GLONASS']['Band_1']['C1C']['rms_multipath_range1_averaged'], decimal=4)
+#     assert_almost_equal(expected_result['GLONASS']['Band_1']['C1C']['elevation_weighted_average_rms_multipath_range1'], result['GLONASS']['Band_1']['C1C']['elevation_weighted_average_rms_multipath_range1'], decimal=4)
+#     assert_almost_equal(expected_result['Galileo']['Band_1']['C1X']['rms_multipath_range1_averaged'], result['Galileo']['Band_1']['C1X']['rms_multipath_range1_averaged'], decimal=4)
+#     assert_almost_equal(expected_result['Galileo']['Band_1']['C1X']['elevation_weighted_average_rms_multipath_range1'], result['Galileo']['Band_1']['C1X']['elevation_weighted_average_rms_multipath_range1'], decimal=4)
+#     assert_almost_equal(expected_result['BeiDou']['Band_2']['C2X']['rms_multipath_range1_averaged'], result['BeiDou']['Band_2']['C2X']['rms_multipath_range1_averaged'], decimal=4)
+#     assert_almost_equal(expected_result['BeiDou']['Band_2']['C2X']['elevation_weighted_average_rms_multipath_range1'], result['BeiDou']['Band_2']['C2X']['elevation_weighted_average_rms_multipath_range1'], decimal=4)
+
+#     assert_almost_equal(expected_result['GPS']['Band_1']['C1C']['nEstimates'], result['GPS']['Band_1']['C1C']['nEstimates'], decimal=4)
+#     assert_almost_equal(expected_result['GLONASS']['Band_1']['C1C']['nEstimates'], result['GLONASS']['Band_1']['C1C']['nEstimates'], decimal=4)
+#     assert_almost_equal(expected_result['Galileo']['Band_1']['C1X']['nEstimates'], result['Galileo']['Band_1']['C1X']['nEstimates'], decimal=4)
+#     assert_almost_equal(expected_result['BeiDou']['Band_2']['C2X']['nEstimates'], result['BeiDou']['Band_2']['C2X']['nEstimates'], decimal=4)
+
+#     assert expected_result['GPS']['Band_1']['C1C']['cycle_slip_distribution']     == result['GPS']['Band_1']['C1C']['cycle_slip_distribution'] 
+#     assert expected_result['GLONASS']['Band_1']['C1C']['cycle_slip_distribution'] == result['GLONASS']['Band_1']['C1C']['cycle_slip_distribution'] 
+#     assert expected_result['Galileo']['Band_1']['C1X']['cycle_slip_distribution'] == result['Galileo']['Band_1']['C1X']['cycle_slip_distribution'] 
+#     assert expected_result['BeiDou']['Band_2']['C2X']['cycle_slip_distribution']  == result['BeiDou']['Band_2']['C2X']['cycle_slip_distribution'] 
+
+# Run the tests
+if __name__ == '__main__':
+    pytest.main()
+
+
+