Repository: paarnes/GNSS_Multipath_Analysis_Software (language: python)

Issue:
Feature/vectorization
#### The biggest change in this PR is speed improvements. A 24-hour RINEX observation file (30s datarate) that took around 12 min to process before, now takes around 2 minutes due to the use of vectorization. 

1. **Utilize Vectorization for Signal Delay Estimations**
   - Implement vectorization for estimating signal delays to achieve significant speed improvements.

2. **Enhancements in Cycle Slip Detection and Multipath Resetting**
   - Improve cycle slip detection and the resetting of multipath for each continuous arc.

3. **Utilize Vectorization for Conversion/Interpolation Operations**
   - Implement vectorization when converting/interpolating from Kepler elements to Earth-Centered, Earth-Fixed (ECEF) coordinates and when converting GLONASS state vectors to ECEF. This optimization results in major speed improvements.

4. **Utilize Vectorization for Satellite Elevation and Azimuth Computations**
   - Implement vectorization when computing satellite elevation and azimuth angles to achieve better performance.

5. **Fix Bug in RinexNav Class**
   - Address a bug in the RinexNav class where the first epoch for the first satellite was being skipped.

6. **Introduce Typing Hints in Main Function**
   - Enhance code readability and maintainability by incorporating typing hints in the main function.

7. **Add Zstd Compression as User Option**
   - Provide users with the option to choose if one want to use Zstd compression for improved data storage.

8. **Default SNR Analysis (Optional to Turn Off)**
   - Conduct SNR analysis by default, with the flexibility for users to turn it off if needed.

9. **Code Refactoring**
   - Conduct necessary refactoring to enhance code structure and maintainability.





Discussion / hints:

Fix patch:
diff --git a/.github/workflows/pypi-package-upload-pypi.yml b/.github/workflows/pypi-package-upload-pypi.yml
new file mode 100644
index 0000000..6a4205a
--- /dev/null
+++ b/.github/workflows/pypi-package-upload-pypi.yml
@@ -0,0 +1,52 @@
+name: Upload Python Package to PYPI
+
+on:
+  release:
+    types:
+      - created
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v2
+
+    - name: Set up Python
+      uses: actions/setup-python@v2
+      with:
+        python-version: 3.8
+
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install --upgrade setuptools wheel build twine
+
+    - name: Get current version from pyproject.toml
+      id: get_version
+      run: echo "::set-output name=version::$(grep -oP '(?<=version = ")[^"]*' pyproject.toml)"
+
+    - name: Increment version
+      id: increment_version
+      run: echo "::set-output name=version::$(python -c "version='${{ steps.get_version.outputs.version }}'.split('.'); version[-1] = str(int(version[-1]) + 1); print('.'.join(version))")"
+
+    - name: Update version in pyproject.toml
+      run: sed -i "s/version = \".*\"/version = \"${{ steps.increment_version.outputs.version }}\"/" pyproject.toml
+
+    - name: Commit changes to pyproject.toml
+      run: |
+        git config --local user.email "per.helge.aarnes@gmail.com"
+        git config --local user.name "Per Helge Aarnes"
+        git add pyproject.toml
+        git commit -m "Bump version to ${{ steps.increment_version.outputs.version }}"
+
+    - name: Push changes to the master branch
+      run: git push origin HEAD:${{ github.ref }}
+      if: github.ref == 'refs/heads/master'
+
+    - name: Build and upload to PyPI (master branch)
+      if: github.ref == 'refs/heads/master' && github.event_name == 'release' && github.event.action == 'created'
+      run: |
+        python -m build
+        twine upload --repository-url https://upload.pypi.org/legacy/ --username __token__ --password ${{ secrets.PYPI_TOKEN }} dist/*
diff --git a/.gitignore b/.gitignore
index 5a72953..4181f9d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -128,8 +128,9 @@ dmypy.json
 # Pyre type checker
 .pyre/
 
-Admin 
+Admin
 launch.json
+SLETT_*
 
 Output_Files/
 .vscode/
\ No newline at end of file
diff --git a/README.md b/README.md
index b0d98a4..7be4be1 100644
--- a/README.md
+++ b/README.md
@@ -110,11 +110,11 @@ cutoff_elevation_angle = 10 # drop satellites lower than 10 degrees
 nav_data_rate = 60 # desired datarate for ephemerides (to improve speed)
 
 analysisResults = GNSS_MultipathAnalysis(rinex_obs_file,
-										 broadcastNav1=rinNav_file,
-										 include_SNR = True,
-										 outputDir = output_folder,
-										 nav_data_rate = nav_data_rate,
-										 cutoff_elevation_angle = cutoff_elevation_angle)
+					broadcastNav1=rinNav_file,
+					include_SNR = True,
+					outputDir = output_folder,
+					nav_data_rate = nav_data_rate,
+					cutoff_elevation_angle = cutoff_elevation_angle)
 ```
 
 #### Read a RINEX observation file
@@ -122,7 +122,6 @@ analysisResults = GNSS_MultipathAnalysis(rinex_obs_file,
 from gnssmultipath import readRinexObs
 
 rinObs_file = 'OPEC00NOR_S_20220010000_01D_30S_MO_3.04'
-
 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 = \
@@ -136,7 +135,6 @@ GNSS_obs, GNSS_LLI, GNSS_SS, GNSS_SVs, time_epochs, nepochs, GNSSsystems,\
 from gnssmultipath import Rinex_v3_Reader
 
 rinNav_file = 'BRDC00IGS_R_20220010000_01D_MN.rnx'
-
 navdata = Rinex_v3_Reader().read_rinex_nav(rinNav_file, data_rate=60)
 
 ```
@@ -147,8 +145,7 @@ navdata = Rinex_v3_Reader().read_rinex_nav(rinNav_file, data_rate=60)
 from gnssmultipath import PickleHandler
 
 path_to_picklefile = 'analysisResults.pkl'
-
-result_dict = PickleHandler.read_pickle(pickl_file)
+result_dict = PickleHandler.read_pickle(path_to_picklefile)
 
 
 ```
@@ -158,8 +155,7 @@ result_dict = PickleHandler.read_pickle(pickl_file)
 from gnssmultipath import PickleHandler
 
 path_to_picklefile = 'analysisResults.pkl'
-
-result_dict = PickleHandler.read_zstd_pickle(pickl_file)
+result_dict = PickleHandler.read_zstd_pickle(path_to_picklefile)
 
 ```
 
diff --git a/build_and_upload_to_pypi.sh b/build_and_upload_to_pypi.sh
deleted file mode 100644
index a9ce776..0000000
--- a/build_and_upload_to_pypi.sh
+++ /dev/null
@@ -1,53 +0,0 @@
-#!/bin/bash
-
-# Check if .env file exists
-env_file=".env"
-
-if [ ! -f "$env_file" ]; then
-  echo "Error: .env file not found in the current directory."
-  exit 1
-fi
-
-# Read and set environment variables from .env file
-while IFS= read -r line || [[ -n "$line" ]]; do
-  # Ignore lines starting with '#' (comments) or empty lines
-  if [[ "$line" =~ ^\s*# || -z "$line" ]]; then
-    continue
-  fi
-
-  # Split the line into key and value using the '=' delimiter
-  key="${line%%=*}"
-  value="${line#*=}"
-
-  # Remove leading and trailing whitespace from the key and value
-  key="${key// /}"
-  value="${value// /}"
-
-  # Set the environment variable
-  export "$key=$value"
-#   echo "Set environment variable: $key=$value"
-done < "$env_file"
-
-# # Upgrade pip
-py -m pip install --upgrade pip
-
-# Upgrade build
-py -m pip install --upgrade build
-
-# Build the project
-py -m build
-
-# Upgrade twine
-py -m pip install --upgrade twine
-
-# echo "USERNAME: $USERNAME"
-# echo "API_TOKEN_TEST: $API_TOKEN_TEST"
-
-# Upload the distribution files using twine
-
-# Upload to test
-twine upload --repository testpypi dist/*  -u "$USERNAME" -p "$API_TOKEN_TEST"
-
-# Upload to PYPI prod
-# twine upload dist/* -u "$USERNAME" -p "$API_TOKEN"
-
diff --git a/src/gnssmultipath/GNSS_MultipathAnalysis.py b/src/gnssmultipath/GNSS_MultipathAnalysis.py
index dd38213..5ac2c00 100644
--- a/src/gnssmultipath/GNSS_MultipathAnalysis.py
+++ b/src/gnssmultipath/GNSS_MultipathAnalysis.py
@@ -9,6 +9,7 @@
 import warnings
 import logging
 import time
+from typing import Union, List
 import numpy as np
 from tqdm import tqdm
 from gnssmultipath.readRinexObs import readRinexObs
@@ -26,31 +27,30 @@
 
 
 
-def GNSS_MultipathAnalysis(rinObsFilename,
-                          broadcastNav1=None,
-                          broadcastNav2=None,
-                          broadcastNav3=None,
-                          broadcastNav4=None,
-                          sp3NavFilename_1=None,
-                          sp3NavFilename_2=None,
-                          sp3NavFilename_3=None,
-                          desiredGNSSsystems=None,
-                          phaseCodeLimit= None,
-                          ionLimit=None,
-                          cutoff_elevation_angle=None,
-                          outputDir=None,
-                          plotEstimates= None,
-                          plot_polarplot=None,
-                          include_SNR=None,
-                          save_results_as_pickle = True,
-                          nav_data_rate=60,
-                          tLim_R   = None,
-                          tLim_GEC = None,
-                          includeResultSummary= None,
-                          includeCompactSummary=None,
-                          includeObservationOverview=None,
-                          includeLLIOverview= None,
-                          use_LaTex =None
+def GNSS_MultipathAnalysis(rinObsFilename: str,
+                          broadcastNav1: Union[str, None] = None,
+                          broadcastNav2: Union[str, None] = None,
+                          broadcastNav3: Union[str, None] = None,
+                          broadcastNav4: Union[str, None] = None,
+                          sp3NavFilename_1: Union[str, None] = None,
+                          sp3NavFilename_2: Union[str, None] = None,
+                          sp3NavFilename_3: Union[str, None] = None,
+                          desiredGNSSsystems: Union[List[str], None] = None,
+                          phaseCodeLimit: Union[float, int, None] = None,
+                          ionLimit: Union[float, None] = None,
+                          cutoff_elevation_angle: Union[int, None] = None,
+                          outputDir: Union[str, None] = None,
+                          plotEstimates: bool = True,
+                          plot_polarplot: bool = True,
+                          include_SNR: bool = True,
+                          save_results_as_pickle: bool = True,
+                          save_results_as_compressed_pickle: bool = False,
+                          nav_data_rate: int = 60,
+                          includeResultSummary: Union[bool, None] = None,
+                          includeCompactSummary: Union[bool, None] = None,
+                          includeObservationOverview: Union[bool, None] = None,
+                          includeLLIOverview: Union[bool, None] = None,
+                          use_LaTex: bool = True
                           ):
 
     """
@@ -59,16 +59,17 @@ def GNSS_MultipathAnalysis(rinObsFilename,
     Made by: Per Helge Aarnes
     E-mail: per.helge.aarnes@gmail.com
 
-    Based on the Matlab software GNSS_Reciever_QC_2020 made by Bjørn-Eirik Roald.
+    GNSS_MultipathAnalysis is a software for analyzing the multipath effect on Global Navigation Satellite Systems (GNSS) and
+    is based on the MATLAB software "GNSS_Receiver_QC_2020" made by Bjørn-Eirik Roald.
+    This is the main function of the software that, through the help of other functions:
 
-    The main function of the software that, through the help of other functions:
-
-      - reads RINEX 3 observation files
-      - reads SP3 satellite navigation files
+      - reads RINEX observation files
+      - reads SP3 satellite navigation files (if a SP3 file is fed in)
+      - reads rinex navigation files (if a rinex navigation file is fed in)
       - computes elevation angles of satellites for every epoch
-      - makes estimates of multipath, ionospheric delay, and ambiguity slips
-          for all signals in RINEX 3 observation file
-      - plots estimates in graphs, if user desired it to
+      - makes estimates of multipath, ionospheric delay, and cycle slips
+          for all signals in RINEX observation file
+      - plots estimates if user choose to do it
       - computes and stores statistics on estimates
       - writes an output files containing results
 
@@ -76,9 +77,9 @@ def GNSS_MultipathAnalysis(rinObsFilename,
     further more functions. The function called upon directly in
     GNSS_MultipathAnalysis are:
 
-      - readRinexObs304.py
+      - readRinexObs.py
       - computeSatElevations.py
-      - readFrequencyOverview.py
+      - computeSatElevAzimuth_fromNav.py
       - signalAnalysis.py
       - plotEstimates.py
       - make_barplot.py
@@ -86,7 +87,9 @@ def GNSS_MultipathAnalysis(rinObsFilename,
       - writeOutputFile.py
 
     --------------------------------------------------------------------------------------------------------------------------
+    
     INPUTS:
+    ------
 
     rinObsFilename:           string. Path to RINEX 3 observation file
 
@@ -128,12 +131,15 @@ def GNSS_MultipathAnalysis(rinObsFilename,
     include_SNR:              boolean. If not defined, SNR from Rinex obs file will NOT be used (optional)
 
     save_results_as_pickle:   boolean. If True, the results will be stored as dictionary in form of a binary pickle file. Default set to True.
+    
+
+    save_results_as_compressed_pickle : boolean. If True, the results will be stored as dictionary in form of a binary compressed pickle file (zstd compression). Default set to False.
 
 
     nav_data_rate:            integer. The desired data rate of ephemerides given in minutes. Default is 60 min. The purpose with this
-                                parameter is to speed up processing time. Both reading the RINEX navigation file and looping through the
-                                ephemerides aftwerward will be significatnly faster by increasing this value. Note: A too high value will
-                                degrade the accuracy of the interploated satellite coordinates.
+                              parameter is to speed up processing time. Both reading the RINEX navigation file and looping through the
+                              ephemerides aftwerward will be significatnly faster by increasing this value. Note: A too high value will
+                              degrade the accuracy of the interploated satellite coordinates.
 
     includeResultSummary:     boolean. 1 if user desires output file to
                               include more detailed overview of statistics,
@@ -143,9 +149,12 @@ def GNSS_MultipathAnalysis(rinObsFilename,
     includeCompactSummary:    boolean. 1 if user desired output file to
                               include more compact overview og statistics. (optional)
 
-    includeObservationOverview:     boolean. 1 if user desires output file to
-                                      include overview of obseration types observed
-                                      by each satellite. 0 otherwise (optional)
+    includeObservationOverview: boolean. 1 if user desires output file to
+                                include overview of obseration types observed
+                                by each satellite. 0 otherwise (optional)
+                                      
+    use_LaTex:                 boolean. Will use TeX as an interpreter in plots. Default set to true. "Requires TeX installed on computer".
+   
     --------------------------------------------------------------------------------------------------------------------------
     OUTPUTS:
 
@@ -178,36 +187,26 @@ def GNSS_MultipathAnalysis(rinObsFilename,
     if sp3NavFilename_3 is None:
         sp3NavFilename_3 = ""
     if phaseCodeLimit is None:
-        phaseCodeLimit = 0
+        phaseCodeLimit =  4/60 *100
 
     if ionLimit is None:
-        ionLimit = 0
+        ionLimit = 4/60 
 
     if cutoff_elevation_angle is None:
         cutoff_elevation_angle = 0
 
-    ## -- Create output file
+    #  Create output file
     if outputDir is None:
         outputDir = 'Output_Files'
 
     if not os.path.isdir(outputDir):
         os.mkdir(outputDir)
 
-    ## --Unless graph directory already exists, create directory
+    # Unless graph directory already exists, create directory
     graphDir = outputDir + '/Graphs'
     if not os.path.isdir(graphDir):
         os.mkdir(graphDir)
 
-    if plotEstimates is None:
-        plotEstimates = 1
-
-    if plot_polarplot is None:
-        plot_polarplot = 1
-
-    if tLim_R is None:
-        tLim_R = 1800 # 30 min
-    if tLim_GEC is None:
-        tLim_GEC = 7200 # 2 hours
 
     if includeResultSummary is None:
         includeResultSummary = 1
@@ -227,30 +226,28 @@ def GNSS_MultipathAnalysis(rinObsFilename,
     else:
         includeAllGNSSsystems   = 0
 
-    if use_LaTex is None:
-        use_LaTex = True
 
     ## ---  Control of the user input arguments
-    if type(sp3NavFilename_1) != str:
+    if not isinstance(sp3NavFilename_1, str):
         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:
+    if not isinstance(sp3NavFilename_2, str):
         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:
+    if not isinstance(sp3NavFilename_3, str):
         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:
+    if not isinstance(rinObsFilename, str):
         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
@@ -272,6 +269,11 @@ def GNSS_MultipathAnalysis(rinObsFilename,
 
     if not os.path.isfile(sp3NavFilename_3) and len(sp3NavFilename_3) != 0:
         print('WARNING: Third SP3 Navigation file can not be found.\n')
+        
+    
+    # Check for conflicting save options
+    if save_results_as_pickle and save_results_as_compressed_pickle:
+        save_results_as_pickle = False
 
 
     latex_installed = True
@@ -300,14 +302,13 @@ def GNSS_MultipathAnalysis(rinObsFilename,
     ## --- Read observation file
     includeAllObsCodes  = 0
 
-    if include_SNR is None:
-        desiredObsCodes = ["C", "L"] # only code and phase observations
-    elif include_SNR is True:
+    if include_SNR:
         desiredObsCodes = ["C", "L", "S"]
-
-    desiredObsBands = list(np.arange(1,10)) # all carrier bands. Tot 9, but arange stops at 8 -> 10
+    else: 
+        desiredObsCodes = ["C", "L"] # only code and phase observations
 
 
+    desiredObsBands = list(np.arange(1,10)) # all carrier bands. Tot 9, but arange stops at 8 -> 10
     readSS = 1
     readLLI = 1
 
@@ -326,7 +327,7 @@ def GNSS_MultipathAnalysis(rinObsFilename,
             nepochs, time_epochs, max_sat, sp3NavFilename_1, sp3NavFilename_2, sp3NavFilename_3)
     else:
         nav_files = [broadcastNav1,broadcastNav2,broadcastNav3,broadcastNav4]
-        sat_pos, glo_fcn = computeSatElevAzimuth_fromNav(nav_files, approxPosition, GNSS_SVs, GNSS_obs, time_epochs, nav_data_rate, tLim_GEC,tLim_R)
+        sat_pos, glo_fcn = computeSatElevAzimuth_fromNav(nav_files, approxPosition, GNSS_SVs, time_epochs, nav_data_rate)
 
         ## -- Build same struture for satellit elevation angles if broadcast nav defined
         sat_elevation_angles = {}
@@ -334,15 +335,15 @@ def GNSS_MultipathAnalysis(rinObsFilename,
         for sys in np.arange(0,len(GNSSsystems)):
             currentGNSSsystem = GNSSsystems[sys+1]
             if currentGNSSsystem != 'C':
-                sat_elevation_angles[sys] = sat_pos_dummy[currentGNSSsystem]['Elevation'][:,0:37]
+                sat_elevation_angles[sys] = sat_pos_dummy[currentGNSSsystem]['elevation'][:,0:37]
             else:
-                sat_elevation_angles[sys] = sat_pos_dummy[currentGNSSsystem]['Elevation']
+                sat_elevation_angles[sys] = sat_pos_dummy[currentGNSSsystem]['elevation']
 
         ## - Check for missing systems in navigation file, and remove if found
         missing_sys = []
         dummy_GNSSsystems = GNSSsystems.copy()
         for key,sys in dummy_GNSSsystems.items():
-            if len(sat_pos[sys]['Position']) == 0:
+            if len(sat_pos[sys]['position']) == 0:
                 del sat_pos[sys]
                 del GNSSsystems[key]
                 missing_sys.append(sys)
@@ -412,22 +413,18 @@ def GNSS_MultipathAnalysis(rinObsFilename,
 
     ### --- Build the dicture of the dict used for storing results ----
 
-    ## --initialize variable storing total number of observation codes processed
+    ## --Initialize variable storing total number of observation codes processed
     nCodes_Total = 0
-    ## -- initialize results dictionary
+    ## -- Initialize results dictionary
     analysisResults = {}
     analysisResults['nGNSSsystem'] = nGNSSsystems
     analysisResults['GNSSsystems'] = list(GNSSsystems.values())
     # for sys in range(0,nGNSSsystems):
-    for sys in np.arange(0,nGNSSsystems):
-        ## -- Full name of current GNSS system
-        GNSSsystemName = GNSSsystemCode2Fullname[GNSSsystems[sys+1]]
-        ## -- Include full name of current GNSS system
-        analysisResults[GNSSsystemName] =  {}
-        ## -- Initialize dict for current GNSS system
-        current_sys_dict = {}
-        ## -- Initialize observationOverview field as a dict
-        current_sys_dict['observationOverview'] = {}
+    for sys in np.arange(0, nGNSSsystems):
+        GNSSsystemName = GNSSsystemCode2Fullname[GNSSsystems[sys+1]] # Full name of current GNSS system
+        analysisResults[GNSSsystemName] = {}  # Include full name of current GNSS system
+        current_sys_dict = {} # Initialize dict for current GNSS system
+        current_sys_dict['observationOverview'] = {} # Initialize observationOverview field as a dict
         ## -- Extract the possible bands of current GNSS system, example GPS: 1,2,5
         GNSSsystemPossibleBands = GNSSsystem2BandsMap[GNSSsystemName]
         nPossibleBands = len(GNSSsystemPossibleBands)
@@ -620,9 +617,9 @@ def GNSS_MultipathAnalysis(rinObsFilename,
         if sp3NavFilename_1 != '':
             try:
                 sat_pos[currentGNSSsystem] = {}
-                sat_pos[currentGNSSsystem]['Position']  = sat_coordinates[currentGNSSsystem]
-                sat_pos[currentGNSSsystem]['Azimut']    = sat_azimut_angles[sys]
-                sat_pos[currentGNSSsystem]['Elevation'] = sat_elevation_angles[sys]
+                sat_pos[currentGNSSsystem]['position']  = sat_coordinates[currentGNSSsystem]
+                sat_pos[currentGNSSsystem]['azimuth']    = sat_azimut_angles[sys]
+                sat_pos[currentGNSSsystem]['elevation'] = sat_elevation_angles[sys]
             except:
                 pass
 
@@ -699,8 +696,8 @@ def GNSS_MultipathAnalysis(rinObsFilename,
         for sys in analysisResults['GNSSsystems']:
             curr_sys = GNSS_Name2Code[sys]
             try:
-                azimut_currentSys = analysisResults['Sat_position'][curr_sys]['Azimut']
-                elevation_currentSys = analysisResults['Sat_position'][curr_sys]['Elevation']
+                azimut_currentSys = analysisResults['Sat_position'][curr_sys]['azimuth']
+                elevation_currentSys = analysisResults['Sat_position'][curr_sys]['elevation']
                 print('INFO: Making a regular polar plot for showing azimut and elevation angle for each satellite. Please wait...')
                 if use_LaTex:
                     try:
@@ -764,6 +761,11 @@ def GNSS_MultipathAnalysis(rinObsFilename,
         results_name = os.path.join(outputDir, pickle_filename)
         PickleHandler.write_zstd_pickle(analysisResults, results_name)
         print(f'INFO: The analysis results has been written to the file {pickle_filename}.\n')
+    elif save_results_as_compressed_pickle:
+        pickle_filename = 'analysisResults.pkl.zst'
+        results_name = os.path.join(outputDir, pickle_filename)
+        PickleHandler.write_zstd_pickle(analysisResults, results_name)
+        print(f'INFO: The analysis results has been written to the file {pickle_filename}.\n')
 
     end_time = time.time()
     compute_processing_time(start_time, end_time)
diff --git a/src/gnssmultipath/Geodetic_functions.py b/src/gnssmultipath/Geodetic_functions.py
index c7bbbec..6c0a43d 100644
--- a/src/gnssmultipath/Geodetic_functions.py
+++ b/src/gnssmultipath/Geodetic_functions.py
@@ -2,6 +2,8 @@
 import numpy as np
 from datetime import datetime,timedelta
 from datetime import date
+from typing import List, Union, Literal
+from numpy import ndarray
 import warnings
 warnings.filterwarnings("ignore")
 
@@ -95,6 +97,62 @@ def Nrad(a,b,lat):
 
 
 
+def compute_satellite_azimut_and_elevation_angle(X, Y, Z, xm, ym, zm):
+    """
+    Computes the satellites azimute and elevation angle based on satellitte and
+    receiver ECEF-coordinates. Utilizes vectorization (no for loops) for
+    better performance.
+
+    Unit: Degree.
+
+    Parameters
+    ----------
+    X : Satellite X-coordinate (numpy array)
+    Y : Satellite Y-coordinate (numpy array)
+    Z : Satellite Z-coordinate (numpy array)
+    xm : Reciever X-coordinate (numpy array)
+    ym : Reciever Y-coordinate (numpy array)
+    zm : Reciever X-coordinate (numpy array)
+
+    Returns
+    -------
+    az : Azimut angle in degrees  (numpy array)
+    elev: Elevation angle in degrees  (numpy array)
+
+
+    NOT IN USE AT THE MOMENT
+    """
+
+    ## -- WGS 84 ellipsoid:
+    a   =  6378137.0         # semi-major ax
+    b   =  6356752.314245    # semi minor ax
+
+    # Compute latitude and longitude for the receiver
+    lat,lon,h = ECEF2geodb(a,b,xm,ym,zm)
+
+    # Find coordinate difference between satellite and receiver
+    dX = (X - xm)
+    dY = (Y - ym)
+    dZ = (Z - zm)
+
+    # Convert from ECEF to ENU (east,north, up)
+    east, north, up = np.vectorize(ECEF2enu)(lat,lon,dX,dY,dZ)
+
+    # Calculate azimuth angle and correct for quadrants
+    azimuth = np.rad2deg(np.arctan(east/north))
+    azimuth = np.where((east > 0) & (north < 0) | ((east < 0) & (north < 0)), azimuth + 180, azimuth)
+    azimuth = np.where((east < 0) & (north > 0), azimuth + 360, azimuth)
+
+    # Calculate elevation angle
+    elevation = np.rad2deg(np.arctan(up / np.sqrt(east**2 + north**2)))
+    elevation = np.where((elevation <= 0) | (elevation >= 90), np.nan, elevation) # Set elevation angle to NaN if not in the range (0, 90)
+
+    return azimuth, elevation
+
+
+
+
+
 def compute_azimut_elev(X,Y,Z,xm,ym,zm):
     """
     Computes the satellites azimute and elevation angel based on satellitte and
@@ -188,6 +246,14 @@ def date2gpstime(year,month,day,hour,minute,seconds):
     return week, tow
 
 
+def filter_array_on_system(eph_data,system):
+    """
+    Filtering the array with ephemerids based on system code "G", "E" "C" etc.
+    """
+    mask = np.char.startswith(eph_data[:,0],system)
+    return eph_data[mask]
+
+
 def extract_nav_message(data,PRN,tidspunkt):
     """
     Funksjonen samler de satellittene som har de spesifiserte PRN nummrene,
@@ -540,6 +606,44 @@ def gpstime2date(week, tow):
 
 
 
+def gpstime2date_arrays(week: Union[List[int], np.ndarray], tow: Union[List[float], np.ndarray]) -> np.ndarray:
+    """
+    Calculates date from GPS-week number and "time-of-week" to Gregorian calendar.
+    NOT IN USE AT THE MOMENT
+
+    Parameters
+    ----------
+    week : array/list, GPS-week numbers.
+    tow  : array/list, "Time of week" values.
+
+    Returns
+    -------
+    dates : array, An array of dates given in the Gregorian calendar ([year, month, day, hour, min, sec]).
+    """
+    # Convert week and tow to regular Python lists
+    week = week.tolist()
+    tow = tow.tolist()
+    total_seconds = [w * 7 * 24 * 3600 + t for w, t in zip(week, tow)] # Calculate the time since the GPS epoch (6th January 1980) for each week and tow
+    delta = [timedelta(seconds=s) for s in total_seconds] # Calculate the timedelta from the GPS epoch for each week and tow
+    gps_epoch = datetime(1980, 1, 6) # The GPS reference epoch
+    dates = [gps_epoch + d for d in delta] # Calculate the final date by adding the timedelta to the GPS epoch
+    date_array = np.array([[date.year, date.month, date.day, date.hour, date.minute, date.second] for date in dates]) # Create an array by extracting date components
+
+    return date_array
+
+
+
+def date2gpstime_vectorized(gregorian_date_array):
+    """
+    Convert from gregorian dates to GPS time by using
+    np.vectorization.
+    """
+    years, months, days, hours, minutes, seconds = gregorian_date_array.astype(float).astype(int).T
+    weeks, tows = np.vectorize(date2gpstime)(years, months, days, hours, minutes, seconds)
+    return np.round(weeks), np.round(tows)
+
+
+
 def utc_to_gpst(t_utc):
     """
     Convert from UTC to GPST by adding leap seconds.
diff --git a/src/gnssmultipath/RinexNav.py b/src/gnssmultipath/RinexNav.py
index bef04ea..75fbfb6 100644
--- a/src/gnssmultipath/RinexNav.py
+++ b/src/gnssmultipath/RinexNav.py
@@ -69,8 +69,9 @@ def filter_ephemeris_data_on_time(self, ephemeris_list, time_difference_minutes=
                 prev_sat = current_sat
 
             # Calculate the time difference between the current epoch and the previous one
+            # Appends the data if the time diff is greater than the set limit, or if it is the first epoch (filtered_data is empty)
             time_diff = epoch - previous_epoch
-            if time_diff >= time_difference:
+            if time_diff >= time_difference or not filtered_data: 
                 filtered_data.append(ephemeris_sublist)
                 previous_epoch = epoch
 
@@ -413,5 +414,5 @@ def read_rinex_nav(self, filename, desired_GNSS: list = ['G','R','E','C'], dataf
 
 if __name__=="__main__":
     # brod1 = r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\TestData\NavigationFiles\BRDC00IGS_R_20220010000_01D_MN.rnx"
-    # data = Rinex_v3_Reader().read_rinex_nav(brod1, dataframe=True, data_rate=1200000000000)
+    # data = Rinex_v3_Reader().read_rinex_nav(brod1, dataframe=True, data_rate=60)
     pass
diff --git a/src/gnssmultipath/SatelliteEphemerisToECEF.py b/src/gnssmultipath/SatelliteEphemerisToECEF.py
new file mode 100644
index 0000000..1dff983
--- /dev/null
+++ b/src/gnssmultipath/SatelliteEphemerisToECEF.py
@@ -0,0 +1,591 @@
+"""
+Module for converting from broadcasted ephemerides to ECEF and interpolated to the desired epoch.
+
+Made by: Per Helge Aarnes
+E-mail: per.helge.aarnes@gmail.com
+"""
+
+
+from typing import Literal, Optional, List, Union
+import numpy as np
+from numpy import ndarray
+from tqdm import tqdm
+from gnssmultipath.Geodetic_functions import date2gpstime_vectorized, get_leap_seconds, gpstime2date_arrays, ECEF2enu, ECEF2geodb
+from gnssmultipath.RinexNav import Rinex_v3_Reader
+
+
+
+class GLOStateVec2ECEF:
+
+    """
+    Class for interpolating a GLONASS state vector using 4th order Runge Kutta
+    """
+    def __init__(self):
+        self.sys_code = "R" # glonass system code
+
+
+    def interpolate_glonass_coord_runge_kutta(self, filtered_eph_data, time_epochs):
+        """
+        Function that use broadcast ephemerides (from rinex nav file) and interpolates to the
+        desired time defined by "time_epochs" using 4th order Runge-kutta.
+
+        INPUT:
+        ------
+        ephemerides: A array with ephemerides for current satellite
+        time_epochs: An n_obs X 2 sized array with weeks and tow read from the rinex observation file.
+
+        Return:
+        -------
+        pos            : The interpolated satellite coordiantes for the desired time defined by "time_epochs"
+        vel            : The satellite velocity for the desired time defined by "time_epochs"
+        clock_err      : The satellite clock error for the desired time defined by "time_epochs"
+        clock_rate_err : The satellite clock rate error for the desired time defined by "time_epochs"
+
+        """
+        filtered_eph_data[:,0] = np.nan #reming char from string (ex G01 -> 9999)
+        filtered_eph_data = filtered_eph_data.astype(float)
+
+        ## Read in data:
+        week,toc = date2gpstime_vectorized(filtered_eph_data[:, 1:7].astype(int))
+        tauN = filtered_eph_data[:,6]   # SV clock bias (sec) (-TauN)
+        gammaN = filtered_eph_data[:,7] # SV relative frequency bias (+GammaN)
+
+        tow_rec = np.round(time_epochs, 6)  # extracting tow
+        x_te = filtered_eph_data[:,10]  # X-coordinates at t_e in PZ-90 [km]
+        y_te = filtered_eph_data[:,14]  # Y-coordinates at t_e in PZ-90 [km]
+        z_te = filtered_eph_data[:,18]  # Z-coordinates at t_e in PZ-90 [km]
+
+        vx_te = filtered_eph_data[:,11] # Velocity component at t_e in PZ-90 (v_x) [km/s]
+        vy_te = filtered_eph_data[:,15] # Velocity component at t_e in PZ-90 (v_y) [km/s]
+        vz_te = filtered_eph_data[:,19] # Velocity component at t_e in PZ-90 (v_z) [km/s]
+
+        J_x = filtered_eph_data[:,12]   # Moon and sun acceleration at t_e [km/sec**2]
+        J_y = filtered_eph_data[:,16]   # Moon and sun acceleration at t_e [km/sec**2]
+        J_z = filtered_eph_data[:,20]   # Moon and sun acceleration at t_e [km/sec**2]
+
+        ## -- Convert from UTC to GPST by adding leap seconds
+        leap_sec = get_leap_seconds(week[0],toc[0]) # Get correct leap sec based on date (using first epoch)
+        toc_gps_time = toc + leap_sec # convert from UTC to GPST
+
+        ## -- Find time difference
+        tdiff = tow_rec - toc_gps_time
+
+        ## -- Clock correction (except for general relativity which is applied later)
+        clock_err = tauN + tdiff * (gammaN)
+        clock_rate_err = gammaN
+
+        # Vectorized initialization
+        init_state = np.array([x_te, y_te, z_te, vx_te, vy_te, vz_te]).T  # Transpose to make it a 2D array
+        init_state = 1000*init_state        # converting from km to meters
+        acc = 1000*np.array([J_x, J_y,J_z]) # converting from km to meters
+        state = init_state
+        tstep = 90
+        tt = np.where(tdiff < 0, -tstep, tstep)
+        while np.any(np.abs(tdiff) > 1e-9):
+            tt = np.where(np.abs(tdiff) < tstep, tdiff, tt)
+            k1 = self.glonass_diff_eq(state, acc)
+            k2 = self.glonass_diff_eq(state + k1 * tt[:, None] / 2, -acc)
+            k3 = self.glonass_diff_eq(state + k2 * tt[:, None] / 2, -acc)
+            k4 = self.glonass_diff_eq(state + k3 * tt[:, None], -acc)
+            state += (k1 + 2 * k2 + 2 * k3 + k4) * tt[:, None] / 6.0
+            tdiff -= tt
+
+        pos = state[:, :3]
+        vel = state[:, 3:6]
+        return pos, vel, clock_err, clock_rate_err
+
+
+    def glonass_diff_eq(self, state:ndarray , acc:ndarray):
+        """
+        Calculate the derivatives of Glonass satellite coordinates using the Runge-Kutta interpolation method.
+
+        INPUT:
+        ------
+        state : numpy.ndarray. State vector containing the coordinates and velocities of Glonass satellites.
+                Shape: (n, 6), where n is the number of epochs.
+
+        acc   : numpy.ndarray. Acceleration vector. Shape: (3, n), where n is the number of epochs.
+
+        Returns:
+        --------
+        der_state : numpy.ndarray. Derivatives of the state vector representing the change in satellite coordinates and velocities.
+                    Shape: (n, 6), where n is the number of epochs.
+
+        Constants:
+        -----------
+        J2    : float. Second zonal coefficient of the spherical harmonic expression.
+        mu    : float. Gravitational constant [m^3/s^2] (product of the mass of the Earth and the gravitational constant).
+        omega : float. Earth rotation rate [rad/s].
+        ae    : float. Semi-major axis of the PZ-90 ellipsoid [m].
+        """
+
+        J2 = 1.0826257e-3     # Second zonal coefficient of spherical harmonic expression.
+        mu = 3.9860044e14     # Gravitational constant [m3/s2]   (product of the mass of the earth and and gravity constant)
+        omega = 7.292115e-5   # Earth rotation rate    [rad/sek]
+        ae = 6378136.0        # Semi-major axis PZ-90   [m]
+
+        r = np.linalg.norm(state[:, :3], axis=1)  # Euclidean norm for the radius
+        der_state = np.zeros((state.shape[0], 6))
+
+        zero_indices = r**2 < 0
+        der_state[zero_indices] = 0  # Set derivatives to zero for cases where r^2 < 0
+
+        a = 1.5 * J2 * mu * (ae**2) / (r**5)
+        b = 5 * (state[:, 2]**2) / (r**2)
+        c = -mu / (r**3) - a * (1 - b)
+
+        der_state[:, :3] = state[:, 3:]
+        der_state[:, 3] = (c + omega**2) * state[:, 0] + 2 * omega * state[:, 4] + acc[0, :]
+        der_state[:, 4] = (c + omega**2) * state[:, 1] - 2 * omega * state[:, 3] + acc[1, :]
+        der_state[:, 5] = (c - 2 * a) * state[:, 2] + acc[2, :]
+
+        return der_state
+
+
+    def get_leap_seconds(self, week,tow):
+        """
+        Add leap seconds based on date. Input is week and tow for current obs.
+        """
+        year,month,day,_,_,_ = gpstime2date_arrays(week, tow) # convert to gregorian date
+        ep_time = (year,month,day)
+        if ep_time <= (2006, 1, 1):
+            raise ValueError("Have no history on leap seconds before 2006")
+        elif ep_time <= (2009, 1, 1):
+            return 14
+        elif ep_time <= (2012, 7, 1):
+            return 15
+        elif ep_time <= (2015, 7, 1):
+            return 16
+        elif ep_time <= (2017, 1, 1):
+            return 17
+        else:
+            return 18
+
+
+
+
+
+
+
+
+
+
+class Kepler2ECEF:
+    """
+    Class for converting from Kepler elements from a
+    RINEX navigation file to Earth Centered Earth Fixed (ECEF).
+
+    Initialize the class with a RINEX navigation file, the reciecer coordinates in ECEF
+    (x_rec, y_rec, z_rec). data_rate is optional
+    and set to 60 min as default. (60 min between each epoch)
+
+    """
+
+    def __init__(self, x_rec, y_rec, z_rec, data_rate=60):
+        self.x_rec = x_rec
+        self.y_rec = y_rec
+        self.z_rec = z_rec
+        self.data_rate = data_rate
+
+    def kepler2ecef(self, filtered_eph_data, tow_rec):
+        """
+        The function calculates satellite coordinates and corrects for Earth's rotation. It transforms from Kepler elements to Earth-Centered, Earth-Fixed (ECEF) coordinates.
+        Supports GPS,Galileo and BeiDou
+        """
+        GM         = 3.986005e14      # Product of Earth's mass and the gravitational constant
+        omega_e    = 7.2921151467e-5  # Earth's angular velocity [rad/second]
+        c          = 299792458        # Speed of light [m/s]
+
+        filtered_eph_data[:,0] = np.nan  # Remove cell containing a string representing PRN with system code (to be able to convert array to float)
+        filtered_eph_data = filtered_eph_data.astype(float)
+
+        # Extract data:
+        M0         = filtered_eph_data[:,13] # Mean anomaly at reference epoch
+        delta_n    = filtered_eph_data[:,12] # Mean motion difference from computed value [rad/second]
+        e          = filtered_eph_data[:,15] # Eccentricity (e)
+        A          = filtered_eph_data[:,17]**2  # Square root of semimajor axis in (sqrt(meter))
+        OMEGA      = filtered_eph_data[:,20] # Longitude of ascending node (OMEGA0) in radians.
+        i0         = filtered_eph_data[:,22] # Inclination angle at reference time (i0) [rad]
+        omega      = filtered_eph_data[:,24] # Argument of perigee (omega) in [rad]
+        OMEGA_dot  = filtered_eph_data[:,25] # Rate of right ascension (OMEGA_dot) [rad/second]
+        i_dot      = filtered_eph_data[:,26] # Rate of inclination angle (i_dot) [rad/second]
+        C_uc       = filtered_eph_data[:,14] # cosine harmonic correction coefficients (argument of perigee)
+        C_us       = filtered_eph_data[:,16] # sine harminic Correction coefficients (argument of perigee)
+        C_rc       = filtered_eph_data[:,23] # cosine harmonic correction term to the orbit radius (C_rc)
+        C_rs       = filtered_eph_data[:,11] # sine harmonic correction term to the orbit radius (C_rs)
+        C_ic       = filtered_eph_data[:,19] # cosine harmonic correction term to the angle of inclination (C_ic)
+        C_is       = filtered_eph_data[:,21] # sine harmonic correction term to the angle of inclination (C_is)
+        toe        = filtered_eph_data[:,18] # Reference "time of ephemeris" data (toe) (sec of GPS week)
+
+        n0  = np.sqrt(GM/A**3) # Computed mean motion (n0) [rad/second]
+        t_k = tow_rec - toe # Time difference between reception time and ephemeris reference time
+
+        n_k = n0 + delta_n   # Corrected mean motion (n_k)
+        M_k = M0 + n_k*t_k   # Mean anomaly (M_k) [rad]
+
+        # Calculate eccentric anomaly
+        E = M_k.copy()
+        for _ in range(10):
+            E = M_k + e * np.sin(E)
+            dE = np.fmod(E - M_k, 2 * np.pi)
+            mask = abs(dE) < 1.e-12
+            if np.all(mask):
+                break
+
+        cosv = (np.cos(E) - e)/(1 - e*np.cos(E))
+        sinv = (np.sqrt(1 - e**2)*np.sin(E))/(1-e*np.cos(E))
+        tanv = sinv/cosv
+
+        # Quadrant correction
+        v = np.arctan(tanv)
+        mask = (sinv > 0) & (cosv < 0) | (sinv < 0) & (cosv < 0)
+        v[mask] += np.pi
+        mask = sinv < 0
+        v[mask] += 2 * np.pi
+
+        theta   = v + omega
+        # Orbital parameter corrections
+        du_k    = C_uc*np.cos(2*theta) + C_us*np.sin(2*theta)
+        dr_k    = C_rc*np.cos(2*theta) + C_rs*np.sin(2*theta)
+        di_k    = C_ic*np.cos(2*theta) + C_is*np.sin(2*theta)
+
+        # Corrected orbital parameters
+        u_k     = theta + du_k
+        r_k     = A*(1 - e*np.cos(E)) + dr_k
+        i_k     = i0 + i_dot*t_k + di_k
+
+        # Corrected longitude of ascending node
+        OMEGA_k = OMEGA + (OMEGA_dot - omega_e)*t_k - omega_e*toe
+
+        # Satellite position in the orbit
+        x = r_k*np.cos(u_k)
+        y = r_k*np.sin(u_k)
+
+        # ECEF coordinates for the satellite
+        X = x*np.cos(OMEGA_k) - y*np.sin(OMEGA_k)*np.cos(i_k)
+        Y = x*np.sin(OMEGA_k) + y*np.cos(OMEGA_k)*np.cos(i_k)
+        Z = y*np.sin(i_k)
+
+        # Relativistic clock correction
+        dT_rel = (-2/c**2)*np.sqrt(A*GM)*e*np.sin(E)
+
+        if (abs(self.x_rec)) > 1.0 and (abs(self.y_rec)) > 1.0 and (abs(self.z_rec)) > 1.0:
+            TRANS = 0
+            TRANS0 = 0.075  # approximate signal travel time
+            j = 0
+            while np.any(abs(TRANS0 - TRANS) > 1e-10):
+                j += 1
+                if (j > 20):
+                    print('Error: The travel time-rotation does not converge!')
+                    break
+
+                TRANS = TRANS0
+                OMEGA_k = OMEGA + (OMEGA_dot - omega_e) * t_k - omega_e * (toe + TRANS)
+                X = x * np.cos(OMEGA_k) - y * np.sin(OMEGA_k) * np.cos(i_k)
+                Y = x * np.sin(OMEGA_k) + y * np.cos(OMEGA_k) * np.cos(i_k)
+                Z = y * np.sin(i_k)
+                dX = (X - self.x_rec)
+                dY = (Y - self.y_rec)
+                dZ = (Z - self.z_rec)
+                DS = np.sqrt(dX**2 + dY**2 + dZ**2)
+                TRANS0 = DS / c
+        else:
+            # Do not correct for the earth rotation
+            OMEGA_k = OMEGA + (OMEGA_dot - omega_e)*t_k - omega_e*toe
+            X = x*np.cos(OMEGA_k) - y*np.sin(OMEGA_k)*np.cos(i_k)
+            Y = x*np.sin(OMEGA_k) + y*np.cos(OMEGA_k)*np.cos(i_k)
+            Z = y*np.sin(i_k)
+
+        return X, Y, Z, dT_rel
+
+
+
+
+
+
+
+class SatelliteEphemerisToECEF:
+    """
+    Parent class for converting from broadcasted ephemerides to Earth Centered Earth Fixed (ECEF).
+    Make use of the classes "GLOStateVec2ECEF" and "Kepler2ECEF".
+
+    Input:
+    -----
+
+    rinex_nav_file  : List or string. Takes in both a single RINEX navfile or a list of RINEX navigation file. If a list is provided, the data will
+                      be merged in a single array (merged on to one file).
+    x_rec           : Receiver X-coordinates in ECEF
+    y_rec           : Receiver Y-coordinates in ECEF
+    z_rec           : Receiver Z-coordinates in ECEF
+    desired_systems : List of desired system codes ["G","R","E"]
+    data_rate       : Update rate of the ephemerides. Higher values correspond to less data and faster file reading (in minutes)
+
+
+    Examples on how to use the class:
+
+    - Compute satellite coordinates for all systems and satellites where the time is defined GPS time "Time-of-week" (seconds):
+            sat_pos = SatelliteEphemerisToECEF(rin_nav_file, x_rec, y_rec, z_rec).get_sat_ecef_coordinates(tow_epochs) # where "tow_epochs" is an 1D array with time-of-week
+
+
+    - Compute satellite coordinates for a specific satellite (in this case "R10") and where the time is defined in GPS time and "Time-of-week" (seconds):
+            sat_pos = SatelliteEphemerisToECEF(rin_nav_file, x_rec, y_rec, z_rec).get_sat_ecef_coordinates(tow_epochs, PRN = 'R10')
+
+
+    - Compute satellite coordinates for a specific satellite (in this case "G20") and where the time is defined in Gregorian time (year, month,day,hour,minute,seconds):
+            time_epochs = np.array([[2020,   10,   30,   13,   22,   14],
+                                    [2020,   10,   30,   13,   22,   15],
+                                    [2020,   10,   30,   13,   22,   16]])
+
+            sat_pos = SatelliteEphemerisToECEF(rin_nav_file, x_rec, y_rec, z_rec).get_sat_ecef_coordinates(time_epochs, time_fmt='GREGORIAN', PRN = 'G20')
+
+
+    - Compute satellite azimuth and elevation angles
+            CONVERTER = SatelliteEphemerisToECEF(rinnav, x_rec, y_rec, z_rec, data_rate=120) # create an object
+            CONVERTER.get_sat_ecef_coordinates(time_epochs)                                  # compute satellite coordinates
+            sat_pos = CONVERTER.compute_satellite_azimut_and_elevation_angle()               # compute satellite elevation and azimuth angles
+
+
+    """
+
+    def __init__(self, rinex_nav_file:Union[str, List[str]], x_rec, y_rec, z_rec, desired_systems: Optional[List[str]] = None, data_rate=60):
+
+        if desired_systems is None:
+            desired_systems = ["G", "R", "E", "C"]
+
+        if isinstance(rinex_nav_file, list):
+            self.ephemerides, self.glo_fcn = self.read_a_list_of_nav_files(rinex_nav_file, data_rate=data_rate)
+        else:
+            self.nav_data = Rinex_v3_Reader().read_rinex_nav(rinex_nav_file, desired_systems, data_rate=data_rate)
+            self.ephemerides = self.nav_data['ephemerides']
+            self.glo_fcn = self.nav_data['glonass_fcn']
+
+        self.x_rec = x_rec
+        self.y_rec = y_rec
+        self.z_rec = z_rec
+        self.max_sat_per_sys = {"G" : 36, "R" : 36, "E" : 36, "C" : 60}
+        self.available_systems = self.find_available_systems_in_eph_data()
+        self.available_systems = list(set(self.available_systems).intersection(desired_systems))
+        self.sat_coord = {sys: {"position": {str(PRN): None for PRN in range(1, max_sat + 1)}} for sys, max_sat in self.max_sat_per_sys.items() if sys in self.available_systems}
+        self.sat_coord_computed = False
+        self.prn_overview = self.get_availible_satellites_for_a_system()
+        self.system_code_mapper = {"G" : "GPS", "R" : "GLONASS", "C" : "BeiDou", "E" : "Galileo"}
+        self.sys_names = [self.system_code_mapper[sys_code] for sys_code in self.available_systems]
+        self.total_sats = sum(len(self.prn_overview[sys])for sys in self.available_systems)
+
+
+    def read_a_list_of_nav_files(self, rinex_nav_file, data_rate):
+        """
+        Reads in a list of RINEX navigation files and merge the data
+        into one data array.
+        """
+        nav_files = [nav_file for nav_file in rinex_nav_file if nav_file != ""]
+        first_nav_data = Rinex_v3_Reader().read_rinex_nav(nav_files[0], data_rate=data_rate)
+        data = np.concatenate([first_nav_data['ephemerides']] + [Rinex_v3_Reader().read_rinex_nav(nav_file, data_rate=data_rate)['ephemerides'] for nav_file in nav_files[1:]], axis=0)
+        glonass_fcn = first_nav_data.get('glonass_fcn', None)
+        data = np.unique(data, axis=0) # ensure no duplicates
+        return data, glonass_fcn
+
+
+    def filter_array_on_PRN(self, PRN:str) -> ndarray:
+        """
+        Filtering the array with ephemerides based on system code "G", "E" "C" or "R"
+        based on PRN number. PRN needs to be in the form "G04".
+        """
+        mask = np.char.startswith(self.ephemerides[:,0], PRN)
+        return self.ephemerides[mask]
+
+
+    def filter_array_on_system(self, sys_code:str) -> ndarray:
+        """ Filtering the array with ephemerides based on system code "G", "E" "C" or "R" """
+        mask = np.char.startswith(self.ephemerides[:,0], sys_code)
+        return self.ephemerides[mask]
+
+    def get_availible_satellites_for_a_system(self) -> dict:
+        """Creates a dict of list of available satellites for each availible system """
+        aviliable_sats = {sys:None for sys in self.available_systems}
+        for sys in self.available_systems:
+            curr_sys = self.filter_array_on_system(sys)
+            aviliable_sats[sys] = list(np.unique(curr_sys[:, 0]))
+        return aviliable_sats
+
+    def find_available_systems_in_eph_data(self):
+        """Returns a list of all avaiable systems in navigation file"""
+        satcodes = self.ephemerides[:, 0]
+        all_sys_codes = [code[0] for code in satcodes if code and code[0].isalpha()]
+        return list(set(all_sys_codes))
+
+    def get_closest_ephemerides_for_PRN_at_time(self, PRN:str, desired_tow:ndarray) -> ndarray:
+        """
+        Collects the ephemerides for the specified PRN number,
+        and retrieves those that are closest to the specified time.
+
+        INPUT:
+        -----
+        ephemerids : array containing ephemerides for all available systems and satellites
+        PRN        : str that specifies which satellite ehemerides shall be extracted for
+        desired_tow: array the time you are looking for a match for
+
+        OUTPUT:
+        ------
+        eph_closest_in_time: array containing ephemerides (closest in time) for the specific PRN
+
+        """
+        ephemerids_filtered = self.filter_array_on_PRN(PRN)
+        if ephemerids_filtered.size == 0:
+            return None
+        epochs_dates = ephemerids_filtered[:, 1:7].astype(int)
+        week_sat, tow_sat = date2gpstime_vectorized(epochs_dates)
+        diff = np.abs(tow_sat[:, np.newaxis] - desired_tow)
+        closest_indices  = np.argmin(diff, axis=0)
+        closest_indices_repeated = closest_indices[:, np.newaxis].repeat(desired_tow.shape[0], axis=1)[:,0]
+        modified_eph_array = ephemerids_filtered[closest_indices_repeated] # repeting array that contains all ephemerides needed in correct order wrt to observation epochs
+        return modified_eph_array
+
+
+
+    def compute_satellite_azimut_and_elevation_angle(self, drop_below_horizon=False):
+        """
+        Computes the satellites azimute and elevation angle based on satellitte and
+        receiver ECEF-coordinates. Utilizes vectorization (no for loops) for
+        better performance.
+
+        drop_below_horizon : Sets values to np.nan for satellites below horizon
+        """
+        # Check if satellite coordinates are computed
+        if not self.sat_coord_computed:
+            raise ValueError('Satellite coordinates are not computed. Please compute coordinates by calling the method "get_sat_ecef_coordinates" before performing this operation.')
+
+
+        ## -- WGS 84 ellipsoid:
+        a   =  6378137.0         # semi-major ax
+        b   =  6356752.314245    # semi minor ax
+
+        # Compute latitude and longitude for the receiver
+        lat,lon,h = ECEF2geodb(a, b, self.x_rec, self.y_rec, self.z_rec)
+
+        bar_format = '{desc}:{percentage:3.0f}%|{bar}|({n_fmt}/{total_fmt} satellites)'
+        desc = ', '.join(self.sys_names[:-1]) + (' and ' + self.sys_names[-1] if len(self.sys_names) > 1 else self.sys_names[0])
+        with tqdm(total = self.total_sats, desc =f"Satellite azimuth and elevation angles are being computed for {desc}" , position=0, leave=True, bar_format=bar_format) as pbar:
+            for sys in self.available_systems:
+                az_array = np.full((self.nepochs, self.max_sat_per_sys[sys]+1), np.nan)
+                el_array = np.full((self.nepochs, self.max_sat_per_sys[sys]+1), np.nan)
+                for sys_code in self.prn_overview[sys]:
+                    PRN = int(sys_code[1::])
+                    if self.sat_coord[sys]['position'][str(PRN)] is None:
+                        pbar.update(1) # Update the progress bar for each satellite processed
+                        continue
+                    X,Y,Z = self.sat_coord[sys]['position'][str(PRN)].T
+                    # Find coordinate difference between satellite and receiver
+                    dX = (X - self.x_rec)
+                    dY = (Y - self.y_rec)
+                    dZ = (Z - self.z_rec)
+
+                    # Convert from ECEF to ENU (east,north, up)
+                    east, north, up = np.vectorize(ECEF2enu)(lat,lon,dX,dY,dZ)
+
+                    # Calculate azimuth angle and correct for quadrants
+                    azimuth = np.rad2deg(np.arctan(east/north))
+                    azimuth = np.where((east > 0) & (north < 0) | ((east < 0) & (north < 0)), azimuth + 180, azimuth)
+                    azimuth = np.where((east < 0) & (north > 0), azimuth + 360, azimuth)
+
+                    # Calculate elevation angle
+                    elevation = np.rad2deg(np.arctan(up / np.sqrt(east**2 + north**2)))
+                    if drop_below_horizon:
+                        elevation = np.where((elevation <= 0) | (elevation >= 90), np.nan, elevation) # Set elevation angle to NaN if not in the range (0, 90)
+                        azimuth = np.where((elevation <= 0) | (elevation >= 90), np.nan, azimuth) # Set elevation angle to NaN if not in the range (0, 90)
+                    # Store azimuth and elevation in the numpy arrays
+                    az_array[:,PRN] = azimuth
+                    el_array[:,PRN] = elevation
+                    pbar.update(1) # Update the progress bar for each satellite processed
+                self.sat_coord[sys]['azimuth'] = az_array
+                self.sat_coord[sys]['elevation'] = el_array
+
+        return self.sat_coord
+
+
+
+    def get_sat_ecef_coordinates(self, desired_time:ndarray, time_fmt:Literal["TOW", "GREGORIAN"] = 'TOW', PRN: Optional[str] = None) -> ndarray:
+        """
+        Main method for converting from Kepler to ECEF and interpolate to the desired time. Input time can be both "Time of Week" (seconds) or
+        gregorian time. This is set by time_format. Set to "TOW" by default.". If no PRN is set, the method will compute the satellite coordiantes for
+        all available systems and satellites.
+
+        INPUT:
+        -----
+        desired_time: The desired time for interpolation. Given in "Time-of-Week" or gregorian time.
+
+                      Example on desired_time input fmt:
+                      ---------------------------------
+                      TOW: array([480134, 480135, 480136, 480137, 480138])
+
+                      GREOGORIAN: array([[2020,   10,   30,   13,   22,   14],
+                                         [2020,   10,   30,   13,   22,   15],
+                                         [2020,   10,   30,   13,   22,   16]])
+
+        time_fmt : Either Time-Of-Week (TOW) or gregorian calender (year, month, day, hh, min, sec)
+        PRN : Optional. The PRN number of the desired satellite included system code in front (ex: "G04" or "E19").
+
+        OUTPUT:
+        ------
+        sat_coord : dict containing X,Y,Z -coordinate for choosen system and satellites
+
+        """
+
+        if time_fmt == "GREGORIAN":
+            desired_time = np.atleast_2d(date2gpstime_vectorized(desired_time)).T
+        if PRN and 'R' not in PRN:
+            filtered_eph_data = self.get_closest_ephemerides_for_PRN_at_time(PRN, desired_time)
+            xs, ys, zs, dTrel = Kepler2ECEF(self.x_rec, self.y_rec, self.z_rec).kepler2ecef(filtered_eph_data,desired_time)
+            return np.array([xs, ys, zs]).T
+        elif PRN and 'R' in PRN:
+            filtered_eph_data = self.get_closest_ephemerides_for_PRN_at_time(PRN, desired_time)
+            current_sat_coord,_,_,_ = GLOStateVec2ECEF().interpolate_glonass_coord_runge_kutta(filtered_eph_data, desired_time)
+            return current_sat_coord
+        else:
+            bar_format = '{desc}:{percentage:3.0f}%|{bar}|({n_fmt}/{total_fmt} satellites)'
+            desc = ', '.join(self.sys_names[:-1]) + (' and ' + self.sys_names[-1] if len(self.sys_names) > 1 else self.sys_names[0])
+            with tqdm(total = self.total_sats, desc =f"Satellite coordinates are being computed for {desc}" , position=0, leave=True, bar_format=bar_format) as pbar:
+                for sys in self.available_systems:
+                    for sys_code in self.prn_overview[sys]:
+                        PRN = int(sys_code[1::])
+                        filtered_eph_data = self.get_closest_ephemerides_for_PRN_at_time(sys_code, desired_time)
+                        if filtered_eph_data is None:
+                            pbar.update(1)
+                            continue
+                        if 'R' in sys_code:
+                            current_sat_coord,_,_,_ = GLOStateVec2ECEF().interpolate_glonass_coord_runge_kutta(filtered_eph_data, desired_time)
+                        else:
+                            xs, ys, zs, dTrel = Kepler2ECEF(self.x_rec, self.y_rec, self.z_rec).kepler2ecef(filtered_eph_data, desired_time)
+                            current_sat_coord = np.array([xs, ys, zs]).T
+
+                        self.sat_coord[sys]['position'][str(PRN)] = current_sat_coord
+                        pbar.update(1) # Update the progress bar for each satellite processed
+        self.sat_coord_computed = True
+        self.nepochs = len(desired_time)
+        return self.sat_coord
+
+
+
+if __name__ == "__main__":
+    # nav1 = r"C:\Users\perhe\Desktop\TEST VEC\test_5\OPEC00NOR_S_20220010000_01D_GN.rnx"
+    # nav2 = r"C:\Users\perhe\Desktop\TEST VEC\test_5\OPEC00NOR_S_20220010000_01D_RN.rnx"
+    # nav3 = r"C:\Users\perhe\Desktop\TEST VEC\test_5\OPEC00NOR_S_20220010000_01D_CN.rnx"
+    # nav4 = r"C:\Users\perhe\Desktop\TEST VEC\test_5\OPEC00NOR_S_20220010000_01D_EN.rnx"
+
+    # rinnav = [nav1,nav2,nav3,nav4]
+    # x_rec = 3149785.9652
+    # y_rec = 598260.8822
+    # z_rec = 5495348.4927
+    # time_epochs = np.array([[2199, 45345]])
+    # time_epochs = np.array([[2020,   10,   30,   13,   22,   14],
+    #                         [2020,   10,   30,   13,   22,   15],
+    #                         [2020,   10,   30,   13,   22,   16]])
+    # time_epochs = np.load(r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\Admin\gammel kode ifm vectorisering\time_epochs.npy")
+    # CONVERTER = SatelliteEphemerisToECEF(rinnav, x_rec, y_rec, z_rec, data_rate=120)
+    # CONVERTER.get_sat_ecef_coordinates(time_epochs[:,1])
+    # sat_pos = CONVERTER.compute_satellite_azimut_and_elevation_angle()
+    # sat_pos = SatelliteEphemerisToECEF(rinnav, x_rec, y_rec, z_rec).get_sat_ecef_coordinates(time_epochs[:,1], PRN = 'R10')
+    # sat_pos = SatelliteEphemerisToECEF(rinnav, x_rec, y_rec, z_rec).get_sat_ecef_coordinates(time_epochs, time_fmt='GREGORIAN' ,PRN = 'R10')
+
+
+    pass
\ No newline at end of file
diff --git a/src/gnssmultipath/__init__.py b/src/gnssmultipath/__init__.py
index e3ae74f..629d5dd 100644
--- a/src/gnssmultipath/__init__.py
+++ b/src/gnssmultipath/__init__.py
@@ -2,4 +2,7 @@
 from .readRinexObs import readRinexObs
 from .RinexNav import Rinex_v3_Reader
 from .RinexNav import Rinex_v2_Reader
-from .PickleHandler import PickleHandler
\ No newline at end of file
+from .PickleHandler import PickleHandler
+from .read_SP3Nav import readSP3Nav
+from .Geodetic_functions import *
+from .SatelliteEphemerisToECEF import SatelliteEphemerisToECEF
\ No newline at end of file
diff --git a/src/gnssmultipath/computeDelayStats.py b/src/gnssmultipath/computeDelayStats.py
index 3674f71..8153c24 100644
--- a/src/gnssmultipath/computeDelayStats.py
+++ b/src/gnssmultipath/computeDelayStats.py
@@ -1,29 +1,38 @@
-from numpy import mean, sqrt, sin, nan,pi,isnan, sum, concatenate, zeros, where, intersect1d,nanmean,count_nonzero,union1d,arange
+"""
+This module is computing the stats.
+
+Made by: Per Helge Aarnes
+E-mail: per.helge.aarnes@gmail.com
+"""
+
+
+from numpy import mean, sqrt, sin, nan,pi,isnan, sum, concatenate, zeros, where, intersect1d,nanmean,count_nonzero,union1d,arange, nanstd
 import warnings
 warnings.filterwarnings("ignore")
 
 def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_angles, range1_slip_periods, ambiguity_slip_periods ,LLI_slip_periods, range1_observations, tInterval):
     """
     Function that computes statistical values on estimates of multipath delay, ionospheric delay and satellite elevation angles
-    
-    --------------------------------------------------------------------------------------------------------------------------
-     INPUTS 
-    
-    ion_delay_phase1:      matrix containing estimates of ionospheric delays 
+
+
+     INPUTS:
+    -------
+
+    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 
+
+     multipath_range1:     matrix containing estimates of multipath delays
                            on the first range signal for each PRN, at each epoch.
-    
+
                            multipath_range1(epoch, PRN)
-    
+
      current_sat_elevation_angles: Array contaning satellite elevation angles at each
-                                   epoch, for current GNSS system. 
-    
+                                   epoch, for current GNSS system.
+
                                    sat_elevation_angles(epoch, PRN)
-    
+
     range1_slip_periods:           dict, 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
@@ -33,174 +42,168 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
                                   every PRN. Each array contains epochs of detected cycle slips (ion & code phase)
                                   starts in first column and slip ends in second
                                   column.
-                                   
-                                  
-    
+
+
+
      LLI_slip_periods_per_sat:     dicct, each cell element contains one matrix for
                                    every PRN. Each matrix contains epochs of LLI
                                    slip starts in first column and LLI slip ends in second
                                    column. These may be often be empty depending on if
                                    RINEX observation files include LLI indicators or
                                    not.
-    
-    
+
+
      range1_observations:  matrix. Contains all range1 observations for all
-                           epochs and for all SVs. 
-    
+                           epochs and for all SVs.
+
                            range1_observations(epoch, PRN)
-    
+
      tInterval:            observation interval in seconds
-    --------------------------------------------------------------------------------------------------------------------------
-     OUTPUTS
-    
-     mean_multipath_range1:        array. contains mean values of estimates of 
+
+
+
+
+     OUTPUTS:
+     -------
+
+
+     mean_multipath_range1:        array. contains mean values of estimates of
                                    multipath delay of the first range signal,
                                    for each satellite
-                                   
-                                   mean_multipath_range1(PRN)                 
-                     
-    
-     overall_mean_multipath_range1: overall mean of all estimates of 
-                                    multipath delay of the first range signal        
-        
-     rms_multipath_range1:         rms values of estimates of multipath delay 
-                                   of the first range signal, for each 
+
+                                   mean_multipath_range1(PRN)
+
+
+     overall_mean_multipath_range1: overall mean of all estimates of
+                                    multipath delay of the first range signal
+
+     rms_multipath_range1:         rms values of estimates of multipath delay
+                                   of the first range signal, for each
                                    satellite.
-    
+
                                    rms_multipath_range1(PRN)
-    
-     average_rms_multipath_range1: average rms of estimates of multipath delay 
-                                   of the first range signal 
-    
+
+     average_rms_multipath_range1: average rms of estimates of multipath delay
+                                   of the first range signal
+
      elevation_weighted_rms_multipath_range1:  rms_multipath_range1, but
                                                weighted based on elevation
                                                angles
-    
+
      elevation_weighted_average_rms_multipath_range1:  average_rms_multipath_range1, but
                                                        weighted based on elevation
                                                        angles
-        
-     mean_ion_delay_phase1:        array. contains mean values of estimates of 
-                                   ionospheric delay on the first phase signal, 
+
+     mean_ion_delay_phase1:        array. contains mean values of estimates of
+                                   ionospheric delay on the first phase signal,
                                    for each satellite
-                                   
-                                   mean_ion_delay_phase1(PRN)  
-    
-     overall_mean_ion_delay_phase1: overall mean of all estimates of 
-                                    ionospheric delay on the first phase signal  
-    
+
+                                   mean_ion_delay_phase1(PRN)
+
+     overall_mean_ion_delay_phase1: overall mean of all estimates of
+                                    ionospheric delay on the first phase signal
+
      mean_sat_elevation_angles:    array. contains mean elevation angle, for
                                    each satellite.
-    
+
                                    mean_sat_elevation_angles(PRN)
-    
-     nEstimates:                   total amount of epochs with multipath estimates, double.  
-    
+
+     nEstimates:                   total amount of epochs with multipath estimates, double.
+
      nEstimates_per_sat:           amount of epochs with multipath estimates
                                    per sat.
-    
+
      nRange1Obs_Per_Sat:           array. each elements says how many range1
                                    observations for that SV
-    
+
                                    nRange1Obs_Per_Sat(PRN)
-    
+
      nRange1Obs:                   total number of range1 observations, all SVs
-    
-     range1_slip_distribution_per_sat:     cell. each element contains a structure for
+
+     range1_slip_distribution_per_sat:     dict. each element contains a structure for
                                    each satellite. each structure has
                                    information on number of range1 slips
                                    distributed over groups depending on the
                                    elevation angle of the satellite at the
                                    time of the slip.
-    
-     range1_slip_distribution:     structure. Stores information on number of 
-                                   range1 slips distributed over groups depending 
+
+     range1_slip_distribution:     dict. Stores information on number of
+                                   range1 slips distributed over groups depending
                                    on the elevation angle of the satellite at the
                                    time of the slip.
-    
-     LLI_slip_distribution_per_sat:     cell. each element contains a structure for
+
+     LLI_slip_distribution_per_sat:     dict. each element contains a structure for
                                    each satellite. each structure has
                                    information on number of LLI detected slips
                                    distributed over groups depending on the
                                    elevation angle of the satellite at the
                                    time of the slip.
-    
-     LLI_slip_distribution:        structure. Stores information on number of 
-                                   LLI detected slips distributed over groups depending 
+
+     LLI_slip_distribution:        dict. Stores information on number of
+                                   LLI detected slips distributed over groups depending
                                    on the elevation angle of the satellite at the
                                    time of the slip.
-    
-     combined_slip_distribution_per_sat:     cell. each element contains a structure for
+
+     combined_slip_distribution_per_sat:     dict. each element contains a structure for
                                    each satellite. each structure has
                                    information on number of slips detected
                                    both by LLI and this siftware analysis,
                                    distributed over groups depending on the
                                    elevation angle of the satellite at the
                                    time of the slip.
-    
-     combined_slip_distribution:   structure. Stores information on number of 
-                                   slips detected both ny LLI and this software 
-                                   analysis, distributed over groups depending 
+
+     combined_slip_distribution:   dict. Stores information on number of
+                                   slips detected both ny LLI and this software
+                                   analysis, distributed over groups depending
                                    on the elevation angle of the satellite at the
                                    time of the slip.
     --------------------------------------------------------------------------------------------------------------------------
     """
     nSat = len(range1_slip_periods)
 
-    ## set all 0 values to NaN so they are excluded from stats calculation
+    # set all 0 values to NaN so they are excluded from stats calculation
     ion_delay_phase1[ion_delay_phase1==0] = nan
     multipath_range1[multipath_range1==0] = nan
     current_sat_elevation_angles[current_sat_elevation_angles==0] = nan
-    
-    ## -- Compute mean
-    mean_multipath_range1 = nanmean(multipath_range1,axis=0)
-    
-    ## overall mean multipath, excluding NaN values
-    overall_mean_multipath_range1 = nanmean(mean_multipath_range1,axis=0)
-    
-    ## RMS multipath of each satellite, excluding NaN
-    rms_multipath_range1 = sqrt(nanmean(multipath_range1*multipath_range1,axis=0)) 
-    
-    ## Average RMS multipath, excluding NaN
-    average_rms_multipath_range1 = sqrt(nanmean(multipath_range1*multipath_range1)) 
-    
-
-    ## -- Weighted RMS multipath 
+
+    mean_multipath_range1 = nanmean(multipath_range1,axis=0)              # Compute mean
+    overall_mean_multipath_range1 = nanmean(mean_multipath_range1,axis=0) # overall mean multipath, excluding NaN values
+    rms_multipath_range1 = nanstd(multipath_range1,axis=0, ddof=0)        # RMS multipath of each satellite, excluding NaN
+    average_rms_multipath_range1 = nanstd(multipath_range1, ddof=0)       # Average RMS multipath, excluding NaN
+
+    #  Weighted RMS multipath
     weights     = current_sat_elevation_angles.copy()
-    crit_weight = 4*sin(30*pi/180)**2 
-    weights     = 4*sin(weights*pi/180)**2 
+    crit_weight = 4*sin(30*pi/180)**2
+    weights     = 4*sin(weights*pi/180)**2
     weights[weights > crit_weight] = 1
     elevation_weighted_multipath_range1 = multipath_range1*weights
-    
-    ## -- RMS multipath of each satellite, excluding NaN
+
+    # RMS multipath of each satellite, excluding NaN
     elevation_weighted_rms_multipath_range1 = sqrt(nanmean(elevation_weighted_multipath_range1*elevation_weighted_multipath_range1,axis=0))
-    
-    ## -- Average RMS multipath, excluding NaN
+    # Average RMS multipath, excluding NaN
     elevation_weighted_average_rms_multipath_range1 = sqrt(nanmean(elevation_weighted_multipath_range1*elevation_weighted_multipath_range1))
-    
-    ## --Ionosphere 
-    ## mean ionospheric delay for each satellite, excluding NaN
-    mean_ion_delay_phase1 = nanmean(ion_delay_phase1,axis=0)
-    
-    ## Overall mean ionospheric delay, excluding NaN
-    overall_mean_ion_delay_phase1 = nanmean(mean_ion_delay_phase1)
-    
+
+    # Ionosphere
+    mean_ion_delay_phase1 = nanmean(ion_delay_phase1,axis=0) # mean ionospheric delay for each satellite, excluding NaN
+    overall_mean_ion_delay_phase1 = nanmean(mean_ion_delay_phase1) # Overall mean ionospheric delay, excluding NaN
+
     ## Average elevation angle for each satellite, excluding NaN
-    dumm1 = (range1_observations!= 0)*1 # multiplying with 1 to get from True/False -> 1/0
+    dumm1 = (range1_observations!= 0)*1 # multiplying with 1 to get from True/False -> 1/0 # commented out 25.11.2023
     dumm2 = (~isnan(range1_observations))*1
     dummy = (dumm1 & dumm2)
     obs_elevations = current_sat_elevation_angles*dummy
     obs_elevations[obs_elevations==0] =nan
     mean_sat_elevation_angles = nanmean(obs_elevations,axis=0)
-    
+    mean_sat_elevation_angles = nanmean(current_sat_elevation_angles,axis=0)
+
     ## -- Amount of epochs with estimates
     nEstimates = sum(~isnan(multipath_range1))
     nEstimates_per_sat = sum(~isnan(multipath_range1),axis=0)
-        
+
     ## -- Cycle slip distribution based on code-phase difference ---
     range1_slip_distribution_per_sat = {}
     range1_slip_distribution         = {}
-    
+
     range1_slip_distribution['n_slips_0_10']   = 0
     range1_slip_distribution['n_slips_10_20']  = 0
     range1_slip_distribution['n_slips_20_30']  = 0
@@ -209,11 +212,11 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
     range1_slip_distribution['n_slips_over50'] = 0
     range1_slip_distribution['n_slips_NaN']    = 0
     range1_slip_distribution['n_slips_Tot']    = 0
-    
+
     ## -- Cycle slip distribution based on Loss of lock indicators---
     LLI_slip_distribution_per_sat = {}
     LLI_slip_distribution         = {}
-    
+
     LLI_slip_distribution['n_slips_0_10']   = 0
     LLI_slip_distribution['n_slips_10_20']  = 0
     LLI_slip_distribution['n_slips_20_30']  = 0
@@ -222,12 +225,12 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
     LLI_slip_distribution['n_slips_over50'] = 0
     LLI_slip_distribution['n_slips_NaN']    = 0
     LLI_slip_distribution['n_slips_Tot']    = 0
-    
-    
+
+
     ## -- Combining cycle slip based on code-phase difference and LLI ---
     combined_slip_distribution_per_sat = {}
     combined_slip_distribution         = {}
-    
+
     combined_slip_distribution['n_slips_0_10']   = 0
     combined_slip_distribution['n_slips_10_20']  = 0
     combined_slip_distribution['n_slips_20_30']  = 0
@@ -236,12 +239,12 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
     combined_slip_distribution['n_slips_over50'] = 0
     combined_slip_distribution['n_slips_NaN']    = 0
     combined_slip_distribution['n_slips_Tot']    = 0
-           
-    
+
+
     ## -- Detected Cycle slip distribution based on both code-phase difference and ionosphere residuals ---
     ambiguity_slip_distribution_per_sat = {}
     ambiguity_slip_distribution         = {}
-    
+
     ambiguity_slip_distribution['n_slips_0_10']   = 0
     ambiguity_slip_distribution['n_slips_10_20']  = 0
     ambiguity_slip_distribution['n_slips_20_30']  = 0
@@ -250,7 +253,7 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
     ambiguity_slip_distribution['n_slips_over50'] = 0
     ambiguity_slip_distribution['n_slips_NaN']    = 0
     ambiguity_slip_distribution['n_slips_Tot']    = 0
-   
+
     for i in arange(0,nSat):
         ## Create struct for current sat
         range1_slip_distribution_per_sat[i] = {}
@@ -261,12 +264,12 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
         for j in arange(0,nSlipPeriods):
             ## if slip period is shorter than 60 seconds
             slip_epochs.append(int(range1_slip_periods[i+1][j,0]))
-        
+
         ## Get elevation angles for every slip of current sat
-        slip_epoch_elevation_angles = current_sat_elevation_angles[slip_epochs, i+1] 
-        
-        ## -- Store number of slips into groups of their elevation angles. Groups are: 0-10, 10-20, 20-30, 30-40, 40-50, >50, and NaN     
-        range1_slip_distribution_per_sat[i]['n_slips_0_10']    = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=0)  & (slip_epoch_elevation_angles <10)]) # bufix
+        slip_epoch_elevation_angles = current_sat_elevation_angles[slip_epochs, i+1]
+
+        ## -- Store number of slips into groups of their elevation angles. Groups are: 0-10, 10-20, 20-30, 30-40, 40-50, >50, and NaN
+        range1_slip_distribution_per_sat[i]['n_slips_0_10']    = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=0)  & (slip_epoch_elevation_angles <10)])
         range1_slip_distribution_per_sat[i]['n_slips_10_20']   = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=10) & (slip_epoch_elevation_angles <20)])
         range1_slip_distribution_per_sat[i]['n_slips_20_30']   = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=20) & (slip_epoch_elevation_angles <30)])
         range1_slip_distribution_per_sat[i]['n_slips_30_40']   = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=30) & (slip_epoch_elevation_angles <40)])
@@ -274,7 +277,7 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
         range1_slip_distribution_per_sat[i]['n_slips_over50']  = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=50)])
         range1_slip_distribution_per_sat[i]['n_slips_NaN']     = len(slip_epoch_elevation_angles[isnan(slip_epoch_elevation_angles)])
         range1_slip_distribution_per_sat[i]['n_slips_Tot']     = len(slip_epoch_elevation_angles)
-        
+
         range1_slip_distribution['n_slips_0_10']   = range1_slip_distribution['n_slips_0_10']     + range1_slip_distribution_per_sat[i]['n_slips_0_10']
         range1_slip_distribution['n_slips_10_20']  = range1_slip_distribution['n_slips_10_20']    + range1_slip_distribution_per_sat[i]['n_slips_10_20']
         range1_slip_distribution['n_slips_20_30']  = range1_slip_distribution['n_slips_20_30']    + range1_slip_distribution_per_sat[i]['n_slips_20_30']
@@ -282,20 +285,20 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
         range1_slip_distribution['n_slips_40_50']  = range1_slip_distribution['n_slips_40_50']    + range1_slip_distribution_per_sat[i]['n_slips_40_50']
         range1_slip_distribution['n_slips_over50'] = range1_slip_distribution['n_slips_over50']   + range1_slip_distribution_per_sat[i]['n_slips_over50']
         range1_slip_distribution['n_slips_NaN']    = range1_slip_distribution['n_slips_NaN']      + range1_slip_distribution_per_sat[i]['n_slips_NaN']
-        range1_slip_distribution['n_slips_Tot']    = range1_slip_distribution['n_slips_Tot']      + range1_slip_distribution_per_sat[i]['n_slips_Tot']        
-        
-        ## Create a dictionary for current sat        
+        range1_slip_distribution['n_slips_Tot']    = range1_slip_distribution['n_slips_Tot']      + range1_slip_distribution_per_sat[i]['n_slips_Tot']
+
+        ## Create a dictionary for current sat
         LLI_slip_distribution_per_sat[i] = {}
         LLI_slip_epochs = []
         nSlipPeriods = len(LLI_slip_periods[i])
-        
+
         ## -- Get all LLI slip epochs of current sat.
         for j in arange(0,nSlipPeriods):
-            LLI_slip_epochs.append(int(LLI_slip_periods[i][j,0])) 
-            
+            LLI_slip_epochs.append(int(LLI_slip_periods[i][j,0]))
+
         ## Get elevation angles for every LLI slip of current sat
-        if not isnan(current_sat_elevation_angles[LLI_slip_epochs, i+1]):
-            LLI_slip_epoch_elevation_angles = current_sat_elevation_angles[LLI_slip_epochs, i+1].astype(int) # tror +1 int her git nan!!!! finn løsning
+        if not isnan(current_sat_elevation_angles[LLI_slip_epochs, i+1]).any():
+            LLI_slip_epoch_elevation_angles = current_sat_elevation_angles[LLI_slip_epochs, i+1].astype(int)
         else:
             LLI_slip_epoch_elevation_angles = current_sat_elevation_angles[LLI_slip_epochs, i+1]
 
@@ -308,7 +311,7 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
         LLI_slip_distribution_per_sat[i]['n_slips_over50']  = len(LLI_slip_epoch_elevation_angles[(LLI_slip_epoch_elevation_angles >=50)])
         LLI_slip_distribution_per_sat[i]['n_slips_NaN']     = len(LLI_slip_epoch_elevation_angles[isnan(LLI_slip_epoch_elevation_angles)])
         LLI_slip_distribution_per_sat[i]['n_slips_Tot']     = len(LLI_slip_epoch_elevation_angles)
-        
+
         LLI_slip_distribution['n_slips_0_10']   = LLI_slip_distribution['n_slips_0_10']     + LLI_slip_distribution_per_sat[i]['n_slips_0_10']
         LLI_slip_distribution['n_slips_10_20']  = LLI_slip_distribution['n_slips_10_20']    + LLI_slip_distribution_per_sat[i]['n_slips_10_20']
         LLI_slip_distribution['n_slips_20_30']  = LLI_slip_distribution['n_slips_20_30']    + LLI_slip_distribution_per_sat[i]['n_slips_20_30']
@@ -317,9 +320,9 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
         LLI_slip_distribution['n_slips_over50'] = LLI_slip_distribution['n_slips_over50']   + LLI_slip_distribution_per_sat[i]['n_slips_over50']
         LLI_slip_distribution['n_slips_NaN']    = LLI_slip_distribution['n_slips_NaN']      + LLI_slip_distribution_per_sat[i]['n_slips_NaN']
         LLI_slip_distribution['n_slips_Tot']    = LLI_slip_distribution['n_slips_Tot']      + LLI_slip_distribution_per_sat[i]['n_slips_Tot']
-        
-        
-        ## -- For ambiguity 
+
+
+        ## -- For ambiguity
         ambiguity_slip_distribution_per_sat[i] = {}
         ambiguity_slip_epochs = []
         ambiguity_nSlipPeriods = len(ambiguity_slip_periods[i+1])
@@ -328,12 +331,12 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
         for j in arange(0,ambiguity_nSlipPeriods):
             ## if slip period is shorter than 60 seconds
             ambiguity_slip_epochs.append(int(ambiguity_slip_periods[i+1][j,0]))
-        
+
         ## Get elevation angles for every slip of current sat
-        slip_epoch_elevation_angles = current_sat_elevation_angles[ambiguity_slip_epochs, i+1] 
-        
-        ## -- Store number of slips into groups of their elevation angles. Groups are: 0-10, 10-20, 20-30, 30-40, 40-50, >50, and NaN     
-        ambiguity_slip_distribution_per_sat[i]['n_slips_0_10']    = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=0)  & (slip_epoch_elevation_angles <10)]) 
+        slip_epoch_elevation_angles = current_sat_elevation_angles[ambiguity_slip_epochs, i+1]
+
+        ## -- Store number of slips into groups of their elevation angles. Groups are: 0-10, 10-20, 20-30, 30-40, 40-50, >50, and NaN
+        ambiguity_slip_distribution_per_sat[i]['n_slips_0_10']    = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=0)  & (slip_epoch_elevation_angles <10)])
         ambiguity_slip_distribution_per_sat[i]['n_slips_10_20']   = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=10) & (slip_epoch_elevation_angles <20)])
         ambiguity_slip_distribution_per_sat[i]['n_slips_20_30']   = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=20) & (slip_epoch_elevation_angles <30)])
         ambiguity_slip_distribution_per_sat[i]['n_slips_30_40']   = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=30) & (slip_epoch_elevation_angles <40)])
@@ -341,7 +344,7 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
         ambiguity_slip_distribution_per_sat[i]['n_slips_over50']  = len(slip_epoch_elevation_angles[(slip_epoch_elevation_angles >=50)])
         ambiguity_slip_distribution_per_sat[i]['n_slips_NaN']     = len(slip_epoch_elevation_angles[isnan(slip_epoch_elevation_angles)])
         ambiguity_slip_distribution_per_sat[i]['n_slips_Tot']     = len(slip_epoch_elevation_angles)
-        
+
         ambiguity_slip_distribution['n_slips_0_10']   = ambiguity_slip_distribution['n_slips_0_10']     + ambiguity_slip_distribution_per_sat[i]['n_slips_0_10']
         ambiguity_slip_distribution['n_slips_10_20']  = ambiguity_slip_distribution['n_slips_10_20']    + ambiguity_slip_distribution_per_sat[i]['n_slips_10_20']
         ambiguity_slip_distribution['n_slips_20_30']  = ambiguity_slip_distribution['n_slips_20_30']    + ambiguity_slip_distribution_per_sat[i]['n_slips_20_30']
@@ -349,24 +352,23 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
         ambiguity_slip_distribution['n_slips_40_50']  = ambiguity_slip_distribution['n_slips_40_50']    + ambiguity_slip_distribution_per_sat[i]['n_slips_40_50']
         ambiguity_slip_distribution['n_slips_over50'] = ambiguity_slip_distribution['n_slips_over50']   + ambiguity_slip_distribution_per_sat[i]['n_slips_over50']
         ambiguity_slip_distribution['n_slips_NaN']    = ambiguity_slip_distribution['n_slips_NaN']      + ambiguity_slip_distribution_per_sat[i]['n_slips_NaN']
-        ambiguity_slip_distribution['n_slips_Tot']    = ambiguity_slip_distribution['n_slips_Tot']      + ambiguity_slip_distribution_per_sat[i]['n_slips_Tot']        
-        
-        
-        
-        
+        ambiguity_slip_distribution['n_slips_Tot']    = ambiguity_slip_distribution['n_slips_Tot']      + ambiguity_slip_distribution_per_sat[i]['n_slips_Tot']
+
+
+
+
         ## create dict for current sat
         combined_slip_distribution_per_sat[i] = {}
-        
+
         ## -- Get all combined slips for current satellite
-        # combined_slip_epochs = intersect1d(slip_epochs, LLI_slip_epochs) ## intersect?? vil vel kobinere begge? union?
         combined_slip_epochs = union1d(slip_epochs, LLI_slip_epochs) ## må være union her??
 
         # get elevation angles for every combined slip of current sat
         if len(combined_slip_epochs) != 0:
             combined_slip_epochs = combined_slip_epochs.astype(int)
             combined_slip_epoch_elevation_angles = current_sat_elevation_angles[combined_slip_epochs, i+1].astype(int)
-            
-            # Store number of LLI slips into groups of their elevation angles. Groups are: 0-10, 10-20, 20-30, 30-40, 40-50, >50, and NaN           
+
+            # Store number of LLI slips into groups of their elevation angles. Groups are: 0-10, 10-20, 20-30, 30-40, 40-50, >50, and NaN
             combined_slip_distribution_per_sat[i]['n_slips_0_10']    = len(combined_slip_epoch_elevation_angles[(combined_slip_epoch_elevation_angles >=0)  & (combined_slip_epoch_elevation_angles <10)])
             combined_slip_distribution_per_sat[i]['n_slips_10_20']   = len(combined_slip_epoch_elevation_angles[(combined_slip_epoch_elevation_angles >=10) & (combined_slip_epoch_elevation_angles <20)])
             combined_slip_distribution_per_sat[i]['n_slips_20_30']   = len(combined_slip_epoch_elevation_angles[(combined_slip_epoch_elevation_angles >=20) & (combined_slip_epoch_elevation_angles <30)])
@@ -375,7 +377,7 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
             combined_slip_distribution_per_sat[i]['n_slips_over50']  = len(combined_slip_epoch_elevation_angles[(combined_slip_epoch_elevation_angles >=50)])
             combined_slip_distribution_per_sat[i]['n_slips_NaN']     = len(combined_slip_epoch_elevation_angles[isnan(combined_slip_epoch_elevation_angles)])
             combined_slip_distribution_per_sat[i]['n_slips_Tot']     = len(combined_slip_epoch_elevation_angles)
-            
+
             combined_slip_distribution['n_slips_0_10']   = combined_slip_distribution['n_slips_0_10']     + combined_slip_distribution_per_sat[i]['n_slips_0_10']
             combined_slip_distribution['n_slips_10_20']  = combined_slip_distribution['n_slips_10_20']    + combined_slip_distribution_per_sat[i]['n_slips_10_20']
             combined_slip_distribution['n_slips_20_30']  = combined_slip_distribution['n_slips_20_30']    + combined_slip_distribution_per_sat[i]['n_slips_20_30']
@@ -386,7 +388,7 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
             combined_slip_distribution['n_slips_Tot']    = combined_slip_distribution['n_slips_Tot']      + combined_slip_distribution_per_sat[i]['n_slips_Tot']
 
         else:
-    
+
             ## Set to zero
             combined_slip_distribution_per_sat[i]['n_slips_0_10']    = 0
             combined_slip_distribution_per_sat[i]['n_slips_10_20']   = 0
@@ -396,7 +398,7 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
             combined_slip_distribution_per_sat[i]['n_slips_over50']  = 0
             combined_slip_distribution_per_sat[i]['n_slips_NaN']     = 0
             combined_slip_distribution_per_sat[i]['n_slips_Tot']     = 0
-            
+
             combined_slip_distribution['n_slips_0_10']   = 0
             combined_slip_distribution['n_slips_10_20']  = 0
             combined_slip_distribution['n_slips_20_30']  = 0
@@ -405,14 +407,13 @@ def computeDelayStats(ion_delay_phase1, multipath_range1, current_sat_elevation_
             combined_slip_distribution['n_slips_over50'] = 0
             combined_slip_distribution['n_slips_NaN']    = 0
             combined_slip_distribution['n_slips_Tot']    = 0
-        
-    
-    ## --Amount of range1 observations 
-    nRange1Obs_Per_Sat = zeros([1, nSat+1])
-    for PRN in arange(0,nSat):
-        PRN = PRN + 1
-        nRange1Obs_Per_Sat[:,PRN] = count_nonzero(range1_observations[:,PRN] != 0)
-    nRange1Obs = sum(nRange1Obs_Per_Sat)
+
+
+    ## --Amount of range1 observations
+    nRange1Obs_Per_Sat = sum(~isnan(range1_observations),axis=0).reshape(1, -1) # reshape to get 2D array
+    nRange1Obs =  sum(~isnan(range1_observations))
+
+
     return mean_multipath_range1, overall_mean_multipath_range1, rms_multipath_range1, average_rms_multipath_range1,\
         mean_ion_delay_phase1, overall_mean_ion_delay_phase1, mean_sat_elevation_angles, nEstimates, nEstimates_per_sat,\
         nRange1Obs_Per_Sat, nRange1Obs, range1_slip_distribution_per_sat, range1_slip_distribution, ambiguity_slip_distribution_per_sat, ambiguity_slip_distribution,LLI_slip_distribution_per_sat, LLI_slip_distribution,\
diff --git a/src/gnssmultipath/computeSatElevAzimuth_fromNav.py b/src/gnssmultipath/computeSatElevAzimuth_fromNav.py
index ef15e50..3a397f6 100644
--- a/src/gnssmultipath/computeSatElevAzimuth_fromNav.py
+++ b/src/gnssmultipath/computeSatElevAzimuth_fromNav.py
@@ -1,236 +1,48 @@
-import numpy as np
-import pandas as pd, re, os
-from tqdm import tqdm
-from gnssmultipath.RinexNav import Rinex_v3_Reader
-from gnssmultipath.Geodetic_functions import *
+"""
+Module for satellite elevations and azimut angles based on broadcasted ephemerides.
 
-def computeSatElevAzimuth_fromNav(navigationFile,approxPosition,GNSS_SVs,GNSS_obs,time_epochs, nav_data_rate, tLim_GEC=None,tLim_R=None):
+Made by: Per Helge Aarnes
+E-mail: per.helge.aarnes@gmail.com
+"""
+
+from gnssmultipath.SatelliteEphemerisToECEF import SatelliteEphemerisToECEF
+
+def computeSatElevAzimuth_fromNav(rinex_nav_file, approxPosition, GNSS_SVs, time_epochs, nav_data_rate):
     """
     A function for computing satellite elevations and azimut angles based on
     broadcasted ephemerides. Support all global navigation systems (GPS,GLONASS,Galileo & BeiDou).
 
     Input:
     -----
-    navigationFile: list of navigation files
+    navigationFile: path to navigation file or list of paths to navigations files
+    approxPosition : Reciever position in ECEF from RINEX obseration file
+    GNSS_SVs : dictionary containing overview of observed satellites for each epoch and for each systems
+                {'G': array([[11., 30., 15., ...,  0.,  0.,  0.],
+                             [11., 30., 15., ...,  0.,  0.,  0.],
+                             [ 9., 31., 19., ...,  0.,  0.,  0.]]),
+                 'R': array([[ 8.,  8., 15., ...,  0.,  0.,  0.],
+                             [ 8.,  8., 15., ...,  0.,  0.,  0.],
+                             [ 8.,  8., 15., ...,  0.,  0.,  0.],
+                             [ 8.,  9.,  2., ...,  0.,  0.,  0.]])}
+    time_epochs : numpy array with GPS time for GNSS observations (week, time-of-week)  
+                  np.array([[  2190.        , 518399.99999988],
+                            [  2190.        , 518429.99999988],
+                            [  2190.        , 518459.99999988]])      
+    nav_data_rate : desired data rate for the broadcasted ephemerides (when reading the rinex nav file)
+
+    Return:
+    ------
+    sat_pos : Dictionary containing satellite posistion in ECEF, in addition to elavation and azimuth angels.
+    glo_fcn : Dictionary containing GLONASS frequency channel numbers extracted from RINEX navigation file
     """
 
-    ##--- Read rinex navigation file/files. If several files defined, the will
-    ## be added to the same data array for further analysis
-
-    nav_list = [i for i in navigationFile if i !=""] #remove "NONE" if exist in list
-    for idx, nav_file in enumerate(nav_list):
-        nav_data = Rinex_v3_Reader().read_rinex_nav(nav_file, data_rate=nav_data_rate)
-        data_ = nav_data['ephemerides']
-        glo_fcn = nav_data['glonass_fcn']
-        if idx == 0:
-            data = data_
-        else:
-            data = np.append(data,data_,axis=0)
-
-    ### ------- Compute satellite coordiantes in ECEF and compute azimut and elevation angels
-    ## -- Extracting approx postion from RINEX obs-file
-    x = float(approxPosition[0])
-    y = float(approxPosition[1])
-    z = float(approxPosition[2])
-
-
-    total_epochs = sum(len(GNSS_obs[sys]) for sys in GNSS_SVs)
-    pbar = tqdm(total=total_epochs, desc="Computing satellite coordinates, azimuth and elevation angles for desired systems", position=0, leave=True, bar_format='{desc}: {percentage:3.0f}%|{bar}| ({n_fmt}/{total_fmt})')
-
-    ## -- Running the computations for each system
-    GNSS_FullName = dict(list(zip(['G','R','E','C'],['GPS','GLONASS','Galileo','BeiDou']))) # dict for mapping from system code to full name
-    sat_pos = {}                           # Dict for storing all data
-    for sys in GNSS_SVs:
-        sat_pos[sys] = {}                  # new dict for each system
-        curr_pos = {}                      # Cells for storing data
-        nepochs = len(GNSS_obs[sys])       # total nr of epochs
-        X = np.zeros([nepochs,61])         # Array for storing X-coordinate
-        Y = np.zeros([nepochs,61])         # Array for storing X-coordinate
-        Z = np.zeros([nepochs,61])         # Array for storing X-coordinate
-        azimut    = np.zeros([nepochs,61]) # Array for storing satellites azimut angle
-        elevation = np.zeros([nepochs,61]) # Array for storing Satellites elevation angle
-        aktuelle_sat_list_all = []         # Dummy list for availible satellites
-        for ep in np.arange(0,nepochs):
-            aktuelle_sat = [PRN for PRN in list(GNSS_SVs[sys][ep][1::].astype(int)) if PRN !=0] ## added [1::] because first ele is nr of satellites
-            aktuelle_sat_list_all.append(aktuelle_sat)
-        aktuelle_sat_list_all = sorted(list({x for l in aktuelle_sat_list_all for x in l})) # finding uniqe value in list in list
+    sat_pos = {}
+    systems = GNSS_SVs.keys()
+    x_rec, y_rec, z_rec = approxPosition.flatten() # Extracting approx postion from RINEX obs-file
+    # Initialize a converter object to convert the broadcast ephemerides to ECEF and interpolate to "time_epochs"
+    CONVERTER = SatelliteEphemerisToECEF(rinex_nav_file, x_rec, y_rec, z_rec, desired_systems=systems, data_rate=nav_data_rate)
+    sat_pos = CONVERTER.get_sat_ecef_coordinates(time_epochs[:,1])
+    sat_pos = CONVERTER.compute_satellite_azimut_and_elevation_angle(drop_below_horizon = True)
+    glo_fcn = CONVERTER.glo_fcn
 
-        ## -- Run through all epochs
-        null_epoch = time_epochs[:,1][0]
-        for epoch in np.arange(0,nepochs):
-            aktuelle_sat_list = [PRN for PRN in list(GNSS_SVs[sys][epoch][1::].astype(int)) if PRN !=0]
-            ## -- Compute satellite coordinates for all availible satellites
-            curr_time = time_epochs[:,1][epoch] # Extracting time for RINEX obs-file
-            df_data = pd.DataFrame(data) # making dataframe of data
-            curr_data = df_data[df_data.iloc[:,0].str.contains(sys)].to_numpy() # extaction data for current system only
-            ## -- Extract satellites ephemerids
-            time_diff = curr_time - null_epoch
-            # Check to speed up the process. Ephemerides get extracted only if time_diff is greater than the pre defined limit tLim_GEC and tLim_R
-            if sys != 'R' and time_diff > tLim_GEC or sys == 'R' and time_diff > tLim_R or time_diff == 0:
-                eph_data = np.empty((36,1))        # Preallocation array for storing ephemerids
-                null_epoch = curr_time # resetting null_epoch
-                for PRN in aktuelle_sat_list_all:
-                    try: ## need try/except for satellites that in rinex obs, but not in nav file. consider adding functionality for his on the top script
-                        eph = extract_nav_message(curr_data,PRN,curr_time)
-                        eph[0] = eph[0][1::] # Removing system letter from number. Ex G10 -> 10
-                    except:
-                        eph = np.nan
-                        continue
-                    eph_data = np.column_stack((eph_data, eph))
-                eph_data = pd.DataFrame(eph_data)
-                eph_data.columns = eph_data.iloc[0].astype(int) # set first row as header
-            ## -- Updating aktuelle_sat_list based on availible sat from both obs file and nav file
-            aktuelle_sat_list = list(set(aktuelle_sat_list) & set(list(eph_data.columns))) # (prevent error message if in obs but not in nav)
-            for PRN in aktuelle_sat_list:
-                ephemerides = eph_data[PRN].to_numpy().reshape(36,1) # passing curr_data instead to get correct system
-                ephemerides = ephemerides.astype(float)
-                if sys != 'R':
-                    ## -- Computing satellite coordinates for GPS,Galileo and BeiDou if in Rinex file
-                    X[epoch,PRN], Y[epoch,PRN], Z[epoch,PRN],_ = Satkoord2(ephemerides, curr_time, x, y, z)
-                else:
-                    ## -- Computing satellite coordinates for GLONASS if in Rinex file
-                    curr_time = time_epochs[epoch]
-                    pos, _, _, _= compute_GLO_coord_from_nav(ephemerides, curr_time)
-                    X[epoch,PRN] = pos[0]
-                    Y[epoch,PRN] = pos[1]
-                    Z[epoch,PRN] = pos[2]
-                ## - Compute azimut and elevation angle
-                azimut[epoch,PRN],elevation[epoch,PRN] = compute_azimut_elev(X[epoch,PRN], Y[epoch,PRN], Z[epoch,PRN], x, y, z)
-                ## -- Assign the computed variable to temporarly dicts for storing results
-                curr_pos[str(PRN)] = np.array([X[:,PRN],Y[:,PRN],Z[:,PRN]]).T
-
-            ## -- Update dictionary with coordinates
-            sat_pos[sys]['Position']  = curr_pos
-            sat_pos[sys]['Azimut']    = azimut
-            sat_pos[sys]['Elevation'] = elevation
-            pbar.update(1)
-    pbar.close()
     return sat_pos, glo_fcn
-
-
-
-def shorten_navigation_file(navigationFile):
-    """
-    Function that is parsing a navigation file and removes system that
-    not GPS,GLONASS, GALILEO and BeiDou. In addition it removes epochs that
-    have same TOC (Time of clock) as the previous epoch. This function is only used
-    to shorted down the reading time of nav file. In addition it will shorten the time
-    it takes to compute satellite coordinates. (Shorter extraction time)
-
-    NO LONGER IN USE
-
-    Parameters:
-    ----------
-    tLim_GEC: Time limit for GPS,Galileo and BeiDou. If the time difference between
-              two epochs in the navigation file is less than this value, it will be removed
-              to speed up the processing time.
-
-    tLim_R  : Time limit for GLONASS. If the time difference between
-              two epochs in the navigation file is less than this value, it will be removed
-              to speed up the processing time. It recommended to have a smalle limit
-              for GLONASS since state vector are less suitet for interpolartion over
-              larger time windows.
-
-    """
-    fd = open(navigationFile,'r')
-    lines = fd.readlines()
-    del_start_indx = []
-    del_end_indx = []
-    pass_header = False
-
-    for idx in np.arange(0,len(lines)):
-        ## -- Checking if file contains 'S' code
-        if re.match(r'S\d+',lines[idx]):
-            del_start_indx.append(idx)
-            IDX = idx + 1
-            while 'S' not in lines[IDX]:
-                IDX = IDX + 1
-                if IDX == len(lines) or re.match(r'[J,S,G,R,I,E,C]\d+',lines[IDX]):
-                    break
-            del_end_indx.append(IDX)
-
-        ## -- Checking if file contains 'I' code
-        if re.match(r'I\d+',lines[idx]):
-            del_start_indx.append(idx)
-            IDX = idx + 1
-            while 'I' not in lines[IDX]:
-                IDX = IDX + 1
-                if IDX == len(lines) or re.match(r'[J,S,G,R,I,E,C]\d+',lines[IDX]):
-                    break
-            del_end_indx.append(IDX)
-        ## -- Checking if file contains 'J' code
-        if re.match(r'J\d+',lines[idx]):
-            del_start_indx.append(idx)
-            IDX = idx + 1
-            while 'J' not in lines[IDX]:
-                IDX = IDX + 1
-                if IDX == len(lines) or re.match(r'[J,S,G,R,I,E,C]\d+',lines[IDX]):
-                    break
-            del_end_indx.append(IDX)
-
-    ## -- Deleting the indencies detected above
-    del_start_indx.sort(reverse=True)
-    del_end_indx.sort(reverse=True)
-    for index,val in enumerate(del_start_indx):
-        idx_start = del_start_indx[index]
-        idx_end = del_end_indx[index]
-        del lines[idx_start:idx_end]
-
-    ## -- Some RINEX navigation files have several lines with same TOC
-    ## This makes the reading and extraction time of ephemerides in connention
-    ## with satelitte coordinates computation much slower. The script under
-    ## is removing epochs if it has the same TOC as the previous epoch.
-    del_start_indx = []
-    del_end_indx = []
-    for idx in np.arange(0,len(lines)):
-        ## -- Test if passed header
-        if 'END OF HEADER' in lines[idx]:
-            pass_header = True
-        ## -- Checking for GPS, Galileo and BeiDou
-        if pass_header==True and idx+8 < len(lines):
-            if re.match(r'G\d+',lines[idx]) or re.match(r'E\d+',lines[idx]) or re.match(r'C\d+',lines[idx]):
-                # line_curr = lines[idx][0:24]
-                line_curr = lines[idx][0:23] # change 29.01.2023 to remove "-" in end of seconds
-                line_curr = [el for el in line_curr.split(" ") if el != ""]
-                # line_next = lines[idx+8][0:24]
-                line_next = lines[idx+8][0:23] # change 29.01.2023 to remove "-" in end of seconds
-                line_next = [el for el in line_next.split(" ") if el != ""]
-                if line_curr == line_next:
-                    IDX_start = idx + 8
-                    IDX_end = IDX_start + 8
-                    del_start_indx.append(IDX_start)
-                    del_end_indx.append(IDX_end)
-        ## -- Checking for GLONASS (less lines in file)
-        if pass_header==True and idx+4 < len(lines):
-            if re.match(r'R\d+',lines[idx]):
-                line_curr = lines[idx][0:23] # change 29.01.2023  to remove "-" in end of seconds
-                line_curr = [el for el in line_curr.split(" ") if el != ""]
-                line_next = lines[idx+4][0:23] # change 29.01.2023  to remove "-" in end of seconds
-                line_next = [el for el in line_next.split(" ") if el != ""]
-                if line_curr == line_next:
-                    IDX_start = idx + 4
-                    IDX_end = IDX_start + 4
-                    del_start_indx.append(IDX_start)
-                    del_end_indx.append(IDX_end)
-
-    ## -- Deleting the indencies that have equal TOC
-    del_start_indx.sort(reverse=True)
-    del_end_indx.sort(reverse=True)
-    for index,val in enumerate(del_start_indx):
-        idx_start = del_start_indx[index]
-        idx_end = del_end_indx[index]
-        del lines[idx_start:idx_end]
-
-
-
-    ## -- Write to file
-    base_path = os.path.split(navigationFile)[0]
-    org_fileName = os.path.basename(navigationFile)
-    temp_fileName = org_fileName.split('.')[0] + '_' + 'temp.'+ org_fileName.split('.')[-1]
-    full_path = os.path.join(base_path,temp_fileName)
-    with open(full_path,'w') as fid:
-        for line in lines:
-            fid.write(line)
-    fid.close()
-    return
-
diff --git a/src/gnssmultipath/computeSatElevations.py b/src/gnssmultipath/computeSatElevations.py
index 734d385..a485aa5 100644
--- a/src/gnssmultipath/computeSatElevations.py
+++ b/src/gnssmultipath/computeSatElevations.py
@@ -19,55 +19,47 @@ def computeSatElevations(GNSS_SVs, GNSSsystems, approxPosition,\
     """
     Function that computes the satellite elevation angles of all satellites
     of all GNSS systems at each epoch of the observations period.
-    --------------------------------------------------------------------------------------------------------------------------
-     INPUTS
+
+     INPUTS:
+     -------
 
      GNSS_SVs:                 Dict containing number of satellites with
                                observations and PRN for each satellite,
                                for each GNSS system
-
                                GNSS_SVs[GNSSsystemIndex][epoch,j]
                                            j=0: (index null) number of observed GPS satellites
                                            j>1: PRN of observed satellites
-
      GNSSsystems:              Dict containing different GNSS
                                systems included in RINEX file. Elements are strings.
                                Must be either "G","R","E" or "C".
-
      approxPosition:           array containing approximate position from rinex
                                observation file header. [X, Y, Z]
-
      nepochs:                  number of epochs with observations in
                                rinex observation file.
-
      time_epochs:              matrix conatining gps-week and "time of week"
                                for each epoch
                                time_epochs(epoch,i),   i=1: week
                                                        i=2: time-of-week in seconds (tow)
-
      max_sat:                  array that stores max satellite PRN number for
                                each of the GNSS systems. Follows same order as GNSSsystems
-
      nav_filename:             string, path and filename of rinex3.xx navigation file.
-
      almanac_nav_filename:     string, path and filename of sen almanac
                                navigation filename for GLONASS
-    --------------------------------------------------------------------------------------------------------------------------
-     OUTPUTS
 
+
+     OUTPUTS:
+     --------
      sat_elevation_angles:     Dict contaning satellite elevation angles at each
                                epoch, for each GNSS system.
                                sat_elevation_angles[GNSSsystemIndex][epoch, PRN]
-
      sat_azimut_angles:        Dict contaning satellite azimuth angles at each
                                epoch, for each GNSS system.
                                sat_azimut_angles[GNSSsystemIndex][epoch, PRN]
-
-
      sat_coordinates:          Dict contaning satellite coordinates for each
                                epoch, for each GNSS system.
                                sat_coordinates[GNSSsystemIndex][epoch, PRN]
-    --------------------------------------------------------------------------------------------------------------------------
+
+
     """
 
     nGNSSsystems = len(GNSSsystems)
@@ -120,8 +112,6 @@ def computeSatElevations(GNSS_SVs, GNSSsystems, approxPosition,\
         nEpochs = nEpochs_1
         epochInterval = epochInterval_1
 
-    satMissingData = []
-
     bar_format = '{desc}: {percentage:3.0f}%|{bar}| ({n_fmt}/{total_fmt})'
     total_epochs = nGNSSsystems * nepochs
     pbar = tqdm(total=total_epochs, desc="Computing satellite coordinates, azimuth and elevation angles for desired systems", position=0, leave=True, bar_format=bar_format)
@@ -139,41 +129,27 @@ def computeSatElevations(GNSS_SVs, GNSSsystems, approxPosition,\
         if sys in navGNSSsystems:
             curr_pos = {}  # dict for storing data
             for epoch in np.arange(0,nepochs):
-                week = time_epochs[epoch,0] # GPS Week of current epoch
-                tow  = time_epochs[epoch,1] # Time of week of current epoch
-                ## -- Satellites in current epoch that should have elevation computed
-                SVs  = GNSS_SVs[sys][epoch,1::][GNSS_SVs[sys][epoch,1::] != 0].astype(int) #extract only nonzero and convert to integer
+                week, tow = time_epochs[epoch]  # "GPS Week" and "Time of week" of current epoch
+                SVs  = GNSS_SVs[sys][epoch,1::][GNSS_SVs[sys][epoch,1::] != 0].astype(int) # Satellites in current epoch that should have elevation computed
                 n_sat = len(SVs)
                 for sat in np.arange(0,n_sat):
                     # Get satellite elevation angle of current sat at current epoch
                     PRN = int(SVs[sat])
-                    elevation_angle, azimut_angle, missing_nav_data, X[epoch,PRN], Y[epoch,PRN], Z[epoch,PRN] = \
+                    elevation_angle, azimut_angle, X[epoch,PRN], Y[epoch,PRN], Z[epoch,PRN] = \
                         get_elevation_angle(GNSSsystems[k+1], SVs[sat], week, tow, sat_positions, nEpochs,\
                        epoch_dates, epochInterval, navGNSSsystems, approxPosition)
                     curr_pos[int(SVs[sat])] = np.array([X[:,PRN],Y[:,PRN],Z[:,PRN]]).T
                     sat_elevation_angles[k][epoch,SVs[sat]] = elevation_angle
                     sat_azimut_angles[k][epoch,SVs[sat]] = azimut_angle
-                    if missing_nav_data:
-                        ## Combine PRN number and GNSS system
-                        SVN = str(GNSSsystems[k]) + str(SVs(sat))
-                        if not any(SVN in satMissingData):
-                            satMissingData.append(SVN)
                 sat_coordinates[sys]  = curr_pos
                 pbar.update(1)
     pbar.close()
-
-    if satMissingData:
-        print('INFO(computeSatElevations): The following satellites had missing orbit data in SP3 file.',\
-        'Their elevation angles were set to 0 for these epochs')
-        print(satMissingData)
-
-
     print('INFO(computeSatElevations): Satellite elevation angles have been computed')
-
     return sat_elevation_angles, sat_azimut_angles, sat_coordinates
 
 
 
+
 def get_elevation_angle(sys, PRN, week, tow, sat_positions, nEpochs, epoch_dates, epochInterval, navGNSSsystems, x_e):
     """
     Calculates elevation angle of a satelite with specified PRN at specified
@@ -183,51 +159,31 @@ def get_elevation_angle(sys, PRN, week, tow, sat_positions, nEpochs, epoch_dates
     ------
 
     sys:              Satellite system, string. ex. "E" or "G"
-
     PRN:              Satellite identification number, integer
-
     week:             GPS-week number, float
-
     tow:              "time-of-week", float
-
-
     sat_positions:    dictionary containing satellite navigation ephemeris of all
                       satellites of observation period, of each GNSS system. The
                       structure is like this sat_positions[systemcode][epoch][PRN].
                       Then you get X,Y and Z coordinates. Ex: sat_positions['G'][100][24]
                       will extract GPS position at epoch 100 for PRN 24.
-
-
     nEpochs:          integer, The number of navigation ephemeris epochs for this satellite.
-
     epoch_dates:      The gregorian date for each epock in SP3 file
-
-
     epochInterval:    float. The epoch interval in seconds. The time difference between epochs.
-
-
-
     navGNSSsystems:   list, conatins codes of GNSS systems with navigation data.
                       ex: ['G', 'R']
-
     x_e:             Coordinates, in ECEF reference frame, of receiver station ex. [X,Y,Z]
 
-
-
     OUTPUTS:
     --------
-
     elevation_angle:  Elevation angle of specified satelite at specified
                       epoch, view from specified receiver station. Unit: Degrees
-
     missing_nav_data: Boolean, 1 if orbit data for current satellite is
                       missing from sp3 file, o otherwise
-
     Xs:               float. The computed X-coordinate
     Ys:               float. The computed Y-coordinate
     Zs:               float. The computed Z-coordinate
 
-
     """
 
     ## -- Define GRS80 ellipsoid parameters
@@ -235,12 +191,10 @@ def get_elevation_angle(sys, PRN, week, tow, sat_positions, nEpochs, epoch_dates
     f       = 1/298.257222100882711243
     b       = a*(1-f)
 
-    missing_nav_data = 0
     ## Get date in form of [year, month, week, day, min, sec] from GPS-week and tow
     date_ = gpstime2date(week, round(tow,1)) ## added round to prevent 59.99999 seconds
     Xs, Ys, Zs = preciseOrbits2ECEF(sys, PRN, date_, epoch_dates, epochInterval, nEpochs, sat_positions, navGNSSsystems)
     if all([Xs,Ys,Zs]) == 0:
-        missing_nav_data = 1
         elevation_angle = 0
         azimut_angle = 0
     else:
@@ -266,4 +220,4 @@ def get_elevation_angle(sys, PRN, week, tow, sat_positions, nEpochs, epoch_dates
             elevation_angle = np.nan
             azimut_angle = np.nan
 
-    return elevation_angle,azimut_angle, missing_nav_data, float(Xs), float(Ys), float(Zs)
+    return elevation_angle,azimut_angle, float(Xs), float(Ys), float(Zs)
diff --git a/src/gnssmultipath/detectCycleSlips.py b/src/gnssmultipath/detectCycleSlips.py
index ca854a8..1fed735 100644
--- a/src/gnssmultipath/detectCycleSlips.py
+++ b/src/gnssmultipath/detectCycleSlips.py
@@ -4,55 +4,75 @@ def detectCycleSlips(estimates, missing_obs_overview,epoch_first_obs, epoch_last
      Function that detects epochs with cycle slips, given test estimates
      and a critical rate of change
     --------------------------------------------------------------------------------------------------------------------------
-    
-     INPUTS:            
+
+     INPUTS:
     --------------------------------------------------------------------------------------------------------------------------
      estimates:            matrix containing estimates from a linear combination
                            for a specific PRN estimates(epoch, PRN)
-    
+
      missing_obs_overview: matrix of size nepochs x nPRN containing 1 or 0.
                            1 indicates that the satellite with this PRN has no
                            estimate at this epoch. There are no 1s before
-                           first observation epoch or after last observation epoch            
+                           first observation epoch or after last observation epoch
                            as lack of estimates are implied here.
-    
+
                            missing_obs_overview(epoch, PRN)
-    
-    
+
+
      epoch_first_obs:      epoch of first estimate for the current satellite
-    
+
      epoch_last_obs:       epoch of last estimate for the current satellite
-    
-     tInterval:            observations interval in seconds. 
-    
+
+     tInterval:            observations interval in seconds.
+
      crit_slip_rate:       critical rate of change of estimate to
-                           indicate an cycle slip. [m/seconds]. 
-    
+                           indicate an cycle slip. [m/seconds].
+
     --------------------------------------------------------------------------------------------------------------------------
-    
+
      OUTPUTS:
-    
+
      slip_epochs:          array, contains epochs with detected cycle slip
     --------------------------------------------------------------------------------------------------------------------------
     """
-    ## -- Preallocating
-    slips_from_missing_obs = []
+
     ## -- Calculate rate of change of estimates of ionospheric delay (time derivative)
-    estimates_rate_of_change = np.diff(estimates)/tInterval
+    estimates_rate_of_change = np.diff(estimates,axis=0)/tInterval
 
-    ## -- Detect slips for current sat. as epochs with either estimates_rate_of_change
+    ## -- Detect slips for epochs with either estimates_rate_of_change
     ## higher than critical value, or epochs with missing estimates
-    slips_from_crit_rate =  [idx for idx,val in enumerate(estimates_rate_of_change) if abs(val) > crit_slip_rate] # use listcomp instead
+    condition_met = np.abs(estimates_rate_of_change) > crit_slip_rate  # Create a boolean array where True indicates the condition is met
+    slips_from_crit_rate = np.where(condition_met)
 
-    if type(epoch_first_obs) !=np.ndarray and type(epoch_last_obs) != np.ndarray:
-        slips_from_missing_obs = (np.where(missing_obs_overview[epoch_first_obs:epoch_last_obs] == 1) + epoch_first_obs)[0].tolist() 
 
-    if len(slips_from_missing_obs) != 0:
-        slip_epochs = np.array(sorted(set(slips_from_crit_rate + slips_from_missing_obs)))
-    else:
-        slip_epochs = np.array(sorted(set(slips_from_crit_rate)))
+    slip_epochs = {str(key): [] for key in range(1, estimates.shape[1])}
+    if slips_from_crit_rate[0].size > 0:
+        for key, value in zip(slips_from_crit_rate[1].astype(str).tolist(), slips_from_crit_rate[0]):
+            slip_epochs[key].append(value)
+
 
-    return slip_epochs 
+    slips_from_missing_obs_curr_sat = []
+    if epoch_first_obs.size > 0:
+        for PRN in np.arange(len(epoch_first_obs)):
+            epoch_first_obs_temp = epoch_first_obs[PRN].astype(int) if not np.isnan(epoch_first_obs[PRN]) else None
+            epoch_last_obs_temp = epoch_last_obs[PRN].astype(int) if not np.isnan(epoch_last_obs[PRN]) else None
+
+            if any(item is None for item in [epoch_first_obs_temp, epoch_last_obs_temp]):
+                continue
+            else:
+                slips_from_missing_obs_curr_sat = (np.where(missing_obs_overview[epoch_first_obs_temp:epoch_last_obs_temp,PRN] == 1) + epoch_first_obs_temp)[0].tolist()
+
+            slips_from_crit_rate_curr = slip_epochs[str(PRN)]
+            if len(slips_from_missing_obs_curr_sat) != 0:
+                slip_epochs[str(PRN)].extend(list(sorted(set(slips_from_crit_rate_curr + slips_from_missing_obs_curr_sat))))
+            else:
+                slip_epochs[str(PRN)].extend(list(sorted(set(slips_from_missing_obs_curr_sat))))
+
+    # Ensure no duplicates
+    for key, value in slip_epochs.items():
+        slip_epochs[key] = list(sorted(set(value)))
+
+    return slip_epochs
 
 
 
@@ -62,29 +82,29 @@ def orgSlipEpochs(slip_epochs):
      Function that takes array of epochs with detected cycle slips and for
      one satellite and organizing them into slip periods
     --------------------------------------------------------------------------------------------------------------------------
-    
+
      INPUTS
-    
+
      slip_epochs:          array, contains epochs with detected cycle slip
     --------------------------------------------------------------------------------------------------------------------------
-    
+
      OUTPUTS
-    
+
      slip_periods:         Matrix that contains the start of periods with cycle slips in
                            the first column and the ends of the same periods in
                            the second column.
-    
+
                            cycle_slip_periods(ambiguity_slip_priod, j),
                            j = 1, ambiguity period start
                            j = 2, ambiguity period ends
-    
+
      n_slip_periods:       amount of slip periods for current satellite
-    
+
     --------------------------------------------------------------------------------------------------------------------------
     """
 
     ## -- If no slips occurs there are no slip periods for this sat.
-    if len(slip_epochs) != 0:    
+    if len(slip_epochs) != 0:
         # dummy is logical. It will be 1 at indices where the following
         # slip epoch is NOT the epoch following the current slip epoch.
         # These will therefor be the indices where slip periods end.
@@ -92,22 +112,22 @@ def orgSlipEpochs(slip_epochs):
         dummy = (np.diff(slip_epochs) != 1) * 1 # multiply with one to get from "False" and "True" to "0" and "1"
         dummy2 = np.where(dummy==1)
         slip_period_ends = np.append(slip_epochs[dummy2], np.array([slip_epochs[-1]]))
-        n_slip_periods = np.sum(dummy) + 1 
-        
-        ## -- Slip_periods = zeros(n_slip_periods,2); 
+        n_slip_periods = np.sum(dummy) + 1
+
+        ## -- Slip_periods = zeros(n_slip_periods,2);
         slip_periods = np.zeros([n_slip_periods,2])
         ## -- Store slip ends
         slip_periods[:,1] = slip_period_ends
         ## -- Store first slip start manually
         slip_periods[0,0] = slip_epochs[0]
-       
+
         ## Insert remaining slip period starts
         for k in range(1,n_slip_periods):
             indx = [x+1 for x, val in enumerate(slip_epochs) if val == slip_periods[k-1, 1]]
             if len(indx) != 0:
                 indx = indx[0]
                 slip_periods[k, 0] = slip_epochs[indx]
-             
+
     else:
         slip_periods = []
         n_slip_periods = 0
diff --git a/src/gnssmultipath/estimateSignalDelays.py b/src/gnssmultipath/estimateSignalDelays.py
index bd8ee14..51024d2 100644
--- a/src/gnssmultipath/estimateSignalDelays.py
+++ b/src/gnssmultipath/estimateSignalDelays.py
@@ -96,35 +96,35 @@ def estimateSignalDelays(range1_Code, range2_Code,phase1_Code, phase2_Code, carr
     --------------------------------------------------------------------------------------------------------------------------
      OUTPUTS
 
-     ion_delay_phase1:     matrix containing estimates of ionospheric delays
+     ion_delay_phase1:     2D array containing estimates of ionospheric delays
                            on the first phase signal for each PRN, at each epoch.
 
-                           ion_delay_phase1(epoch, PRN)
+                           ion_delay_phase1[epoch, PRN]
 
-     multipath_range1:     matrix containing estimates of multipath delays
+     multipath_range1:     2D array 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
+     multipath_range2:     2D array containing estimates of multipath delays
                            on the second range signal for each PRN, at each epoch.
 
-                           multipath_range2(epoch, PRN)
+                           multipath_range2[epoch, PRN]
 
-     range1_slip_periods:  cell, each cell element contains one matrix for
+     range1_slip_periods:  array, 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
+     range1_observations:  2D array. Contains all range1 observations for all
                            epochs and for all SVs.
 
-                           range1_observations(epoch, PRN)
+                           range1_observations[epoch, PRN]
 
-     phase1_observations:  matrix. Contains all phase1 observations for all
+     phase1_observations:  2D array. Contains all phase1 observations for all
                            epochs and for all SVs.
 
-                           range1_observations(epoch, PRN)
+                           phase1_observations[epoch, PRN]
 
      success:              boolean, 1 if no error occurs, 0 otherwise
 
@@ -155,108 +155,75 @@ def estimateSignalDelays(range1_Code, range2_Code,phase1_Code, phase2_Code, carr
                                     # the range1/phase1 signal
 
     ## -- Initialize data matrices
-    ion_delay_phase1        = np.zeros([nepochs, max_sat+1]) # Addin 1 to since python is nullindexed
+    ion_delay_phase1        = np.zeros([nepochs, max_sat+1])
     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): Some observations are missing. Check for missing data in RINEX observation file!')
-                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
+    ambiguity_slip_periods = {int(key): [] for key in range(1,  max_sat+1)} # Initialize cell for storing phase slip periods
+    range1_slip_periods = {int(key): [] for key in range(1,  max_sat+1)} # Initialize cell for storing slip periods for only range1/phase1
+
+    if FDMA_used: # denne er inødnendig etter vektoriseringen
+        carrier_freq1 = carrier_freq1_list
+        carrier_freq2 = carrier_freq2_list
+        alpha = carrier_freq1**2/carrier_freq2**2 # amplfication factor
+
+
+    #  Get observations
+    range1 = create_array_for_current_obscode(GNSS_obs, ismember(obsCodes[currentGNSSsystem],range1_Code))
+    range2 = create_array_for_current_obscode(GNSS_obs, ismember(obsCodes[currentGNSSsystem],range2_Code))
+
+    phase1 = create_array_for_current_obscode(GNSS_obs, ismember(obsCodes[currentGNSSsystem],phase1_Code))*c/carrier_freq1
+    phase2 = create_array_for_current_obscode(GNSS_obs, ismember(obsCodes[currentGNSSsystem],phase2_Code))*c/carrier_freq2
+
+    if all(v is None for v in [range1, range2, phase1, phase2]):
+        print('ERROR(estimateSignalDelays): Some observations are missing. Check for missing data in RINEX observation file!')
+        success = 0
+        return success
+
+
+    ## -- Calculate estimate for Ionospheric delay on phase 1 signal
+    ion_delay_phase1 = 1/(alpha-1)*(phase1-phase2)
+    ## -- Calculate multipath on both code signals
+    multipath_range1 = 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 = phase1 - range1
+
+    # If any of the four observations are missing, ie np.nan, estimate remains 0 for that epoch and satellite
+    missing_obs_overview = find_missing_observation(range1, range2, array3=phase1, array4=phase2)
+    missing_range1_overview = find_missing_observation(range1,phase1)
+
+    ## -- Get first and last epoch with observations for current PRN
+    epoch_first_obs = find_epoch_of_first_obs(ion_delay_phase1)
+    epoch_last_obs = find_epoch_of_last_obs(ion_delay_phase1)
+
+    ## -- Get first and last epoch with range 1 observations for current PRN
+    epoch_first_range1obs = find_epoch_of_first_obs(N1_pseudo_estimate)
+    epoch_last_range1obs =  find_epoch_of_last_obs(N1_pseudo_estimate)
+
+    range1_slip_epochs = detectCycleSlips(N1_pseudo_estimate, missing_range1_overview, epoch_first_range1obs, epoch_last_range1obs, tInterval, phaseCodeLimit) # detect cycle slips for current epoch for range1/phase1 only
+    ionosphere_slip_epochs = detectCycleSlips(ion_delay_phase1, missing_obs_overview, epoch_first_obs, epoch_last_obs, tInterval, ionLimit) # detect cycle slips for current epoch for either range1/phase1 signal, range2/phase2 signal, or both
 
 
     ## -- 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]
-
+        if np.isnan(epoch_first_obs[PRN]):
+            continue  # Skip the current iteration if it's np.nan
+        epoch_first_obs_prn = int(epoch_first_obs[PRN])
+        epoch_last_obs_prn = int(epoch_last_obs[PRN])
 
-        ## -- 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)
-        ## -- 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
+        ambiguity_slip_epochs = np.union1d(np.array(range1_slip_epochs[str(PRN)]),np.array(ionosphere_slip_epochs[str(PRN)])) # Make combined array of slip epochs from both lin. combinations used to detects slips
+        range1_slip_periods[int(PRN)],_ = orgSlipEpochs(np.array(range1_slip_epochs[str(PRN)])) # Organize slips detected on range1/phase1 signal only
+        ambiguity_slip_periods[int(PRN)], n_slip_periods = orgSlipEpochs(ambiguity_slip_epochs) # Orginize combined slips detected on range1/phase1 signal only
 
         ## 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.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]))
+        if len(ambiguity_slip_periods[int(PRN)]) == 0:
+            ion_delay_phase1[epoch_first_obs_prn::, PRN] = ion_delay_phase1[epoch_first_obs_prn::, PRN] - ion_delay_phase1[epoch_first_obs_prn, PRN]
+            multipath_range1[epoch_first_obs_prn::, PRN] = multipath_range1[epoch_first_obs_prn::, PRN] - np.nanmean(multipath_range1[epoch_first_obs_prn::, 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):
@@ -267,22 +234,22 @@ def estimateSignalDelays(range1_Code, range2_Code,phase1_Code, phase2_Code, carr
                     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)
+                    ion_delay_phase1[slip_start:slip_end+1, PRN] = np.nan
+                    multipath_range1[slip_start:slip_end+1, PRN] = np.nan
                     # 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
+            for ambiguity_period in np.arange(0,n_slip_periods+1):
                 if ambiguity_period == 0:
-                    ambiguity_period_start  = epoch_first_obs
-                    ambiguity_period_end    = int(ambiguity_slip_periods[PRN][0,0])   #INK -1 igjen???
+                    ambiguity_period_start  = epoch_first_obs_prn
+                    ambiguity_period_end    = int(ambiguity_slip_periods[PRN][0,0])
                 ## -- 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
+                    ambiguity_period_end         = epoch_last_obs_prn +1
 
                     ## If last epoch with observation is a slip, then there is no  last ambiguity period
-                    if ambiguity_period_start > epoch_last_obs:
+                    if ambiguity_period_start > epoch_last_obs_prn:
                         ambiguity_period_start = []
                         ambiguity_period_end = []
 
@@ -313,17 +280,13 @@ def estimateSignalDelays(range1_Code, range2_Code,phase1_Code, phase2_Code, carr
                     # 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)]
-
+    # Create array of all code and phase observation for the current obscodes
+    range1_observations = create_array_for_current_obscode(GNSS_obs, ismember(obsCodes[currentGNSSsystem],range1_Code))
+    phase1_observations = create_array_for_current_obscode(GNSS_obs, ismember(obsCodes[currentGNSSsystem],phase1_Code))
 
     # 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
+    return ion_delay_phase1, multipath_range1, range1_slip_periods, ambiguity_slip_periods, range1_observations, phase1_observations, success #removed multipath_range2
+
 
 
 def ismember(list_,code):
@@ -335,3 +298,70 @@ def ismember(list_,code):
         indx = indx[0]
     return indx
 
+
+def create_array_for_current_obscode(GNSS_obs, obscode_idx):
+    """
+    Takes inn a dict of numpy arrays and returns a
+    numpy array where the keys are rows in the array.
+    Replaces all zeros with np.nan
+
+    """
+    try:
+        code_array = np.stack(list(GNSS_obs.values()))[:, :, obscode_idx]
+        code_array = np.squeeze(code_array)
+        code_array[code_array == 0] = np.nan
+    except:
+        code_array = None
+    return code_array
+
+
+
+def find_epoch_of_first_obs(ion_delay_phase1):
+    """
+    Finds the first epoch with observation and
+    stores the epoch/row index for each satellite as a single-column array.
+    A column with no observation at all, is stored as  np.nan
+    """
+    not_nan_indices = np.isfinite(ion_delay_phase1)
+    first_epochs_array = np.argmax(not_nan_indices, axis=0).astype(float)
+    no_true_columns = np.all(~not_nan_indices, axis=0) # Find columns with no True values and set their indices to np.nan
+    first_epochs_array[no_true_columns] = np.nan
+    return first_epochs_array
+
+
+def find_epoch_of_last_obs(ion_delay_phase1):
+    """
+    Finds the last epoch with observations and
+    stores the epoch/row index for each satellite as a single-column array.
+    A column with no observation at all, is stored as  np.nan
+    """
+    not_nan_indices = np.isfinite(ion_delay_phase1)
+    reversed_array = np.flip(not_nan_indices, axis=0) # Reverse the boolean array along the rows
+    last_epochs_array = not_nan_indices.shape[0] - 1 - np.argmax(reversed_array, axis=0).astype(float) # Find the indices of the first True element in each column
+    last_epochs_array[~np.any(not_nan_indices, axis=0)] =  np.nan # Mask columns with no True values with np.nan
+    return last_epochs_array
+
+
+
+def find_missing_observation(array1,array2,array3=None,array4=None):
+    """
+    Creates a array with overview of epoch and sats
+    that are missing observations.
+    """
+    if array3 is not None:
+        mask1 = np.isnan(array1).astype(int)
+        mask2 = np.isnan(array2).astype(int)
+        mask3 = np.isnan(array3).astype(int)
+        mask4 = np.isnan(array4).astype(int)
+        # Combine the masks to create the final result
+        missing_obs_overview = np.maximum.reduce([mask1, mask2, mask3, mask4])
+    else:
+        mask1 = np.isnan(array1).astype(int)
+        mask2 = np.isnan(array2).astype(int)
+        # Combine the masks to create the final result
+        missing_obs_overview = np.maximum.reduce([mask1, mask2])
+    return missing_obs_overview
+
+
+
+
diff --git a/src/gnssmultipath/make_polarplot.py b/src/gnssmultipath/make_polarplot.py
index 0db83b5..8a3cb12 100644
--- a/src/gnssmultipath/make_polarplot.py
+++ b/src/gnssmultipath/make_polarplot.py
@@ -25,8 +25,8 @@ def make_polarplot(analysisResults, graph_dir):
         curr_sys = GNSS_Name2Code[system]
         bands_curr_sys = analysisResults[system]['Bands']
         try:
-            sat_elevation = analysisResults['Sat_position'][curr_sys]['Elevation']
-            sat_azimut = analysisResults['Sat_position'][curr_sys]['Azimut']
+            sat_elevation = analysisResults['Sat_position'][curr_sys]['elevation']
+            sat_azimut = analysisResults['Sat_position'][curr_sys]['azimuth']
         except:
             logger.warning("INFO(GNSS_MultipathAnalysis): Polarplot of multipath is not possible for %s. Satellite azimuth and elevation angles are missing.", system)
             continue
@@ -166,8 +166,8 @@ def make_polarplot_SNR(analysisResults, GNSS_obs,GNSSsystems, obsCodes, graphDir
     for system in analysisResults['GNSSsystems']:
         curr_sys = GNSS_Name2Code[system]
         try:
-            sat_elevation = analysisResults['Sat_position'][curr_sys]['Elevation']
-            sat_azimut = analysisResults['Sat_position'][curr_sys]['Azimut']
+            sat_elevation = analysisResults['Sat_position'][curr_sys]['elevation']
+            sat_azimut = analysisResults['Sat_position'][curr_sys]['azimuth']
         except:
             logger.warning(f"INFO(GNSS_MultipathAnalysis): Polarplot of SNR is not possible for {system}. Satellite azimuth and elevation angles are missing." )
             continue
@@ -260,7 +260,7 @@ def plot_SNR_wrt_elev(analysisResults,GNSS_obs, GNSSsystems, obsCodes, graphDir,
     for system in analysisResults['GNSSsystems']:
         curr_sys = GNSS_Name2Code[system]
         try:
-            sat_elevation = analysisResults['Sat_position'][curr_sys]['Elevation']
+            sat_elevation = analysisResults['Sat_position'][curr_sys]['elevation']
             sat_elevation[(sat_elevation < 0) | (sat_elevation > 90)] = np.nan # removes elevation angels when the sat not visable (below the horizon)
         except:
             logger.warning(f"INFO(GNSS_MultipathAnalysis): Plot of SNR wrt elevation angle is not possible for {system}. Satellite azimuth and elevation angles are missing." )
@@ -300,7 +300,7 @@ def plot_SNR_wrt_elev(analysisResults,GNSS_obs, GNSSsystems, obsCodes, graphDir,
             ## -- Subplot 2
             ax[1].grid(True,linewidth=0.3)
             ax[1].set_xlim([0,90])
-            ax[1].set_title("Signal to noise ratio (SNR) as funtion elevation angle for \n Signal: %s (%s)" % (range1_code, system), va='bottom',fontsize=22)
+            ax[1].set_title("Signal to noise ratio (SNR) as funtion of elevation angle for \n Signal: %s (%s)" % (range1_code, system), va='bottom',fontsize=22)
             ax[1].set_ylim(0,np.nanmax(SNR)+10)
             ax[1].set_xlabel('Elevation angle $[^{\circ}]$',fontsize=18,labelpad=10)
             ax[1].set_ylabel('[dB-Hz]',fontsize=18,labelpad=10)
diff --git a/src/gnssmultipath/make_polarplot_dont_use_TEX.py b/src/gnssmultipath/make_polarplot_dont_use_TEX.py
index 0a326a8..3a21ece 100644
--- a/src/gnssmultipath/make_polarplot_dont_use_TEX.py
+++ b/src/gnssmultipath/make_polarplot_dont_use_TEX.py
@@ -25,8 +25,8 @@ def make_polarplot_dont_use_TEX(analysisResults, graph_dir):
         curr_sys = GNSS_Name2Code[system]
         bands_curr_sys = analysisResults[system]['Bands']
         try:
-            sat_elevation = analysisResults['Sat_position'][curr_sys]['Elevation']
-            sat_azimut = analysisResults['Sat_position'][curr_sys]['Azimut']
+            sat_elevation = analysisResults['Sat_position'][curr_sys]['elevation']
+            sat_azimut = analysisResults['Sat_position'][curr_sys]['azimuth']
         except:
             logger.warning("INFO(GNSS_MultipathAnalysis): Polarplot of multipath is not possible for %s. Satellite azimuth and elevation angles are missing.", system)
             continue
@@ -164,8 +164,8 @@ def make_polarplot_SNR_dont_use_TEX(analysisResults, GNSS_obs,GNSSsystems, obsCo
     for system in analysisResults['GNSSsystems']:
         curr_sys = GNSS_Name2Code[system]
         try:
-            sat_elevation = analysisResults['Sat_position'][curr_sys]['Elevation']
-            sat_azimut = analysisResults['Sat_position'][curr_sys]['Azimut']
+            sat_elevation = analysisResults['Sat_position'][curr_sys]['elevation']
+            sat_azimut = analysisResults['Sat_position'][curr_sys]['azimuth']
         except:
             logger.warning(f"INFO(GNSS_MultipathAnalysis): Polarplot of SNR is not possible for {system}. Satellite azimuth and elevation angles are missing." )
             continue
@@ -258,7 +258,7 @@ def plot_SNR_wrt_elev_dont_use_TEX(analysisResults,GNSS_obs, GNSSsystems, obsCod
     for system in analysisResults['GNSSsystems']:
         curr_sys = GNSS_Name2Code[system]
         try:
-            sat_elevation = analysisResults['Sat_position'][curr_sys]['Elevation']
+            sat_elevation = analysisResults['Sat_position'][curr_sys]['elevation']
             sat_elevation[(sat_elevation < 0) | (sat_elevation > 90)] = np.nan # removes elevation angels when the sat not visable (below the horizon)
         except:
             logger.warning("INFO(GNSS_MultipathAnalysis): Plot of SNR wrt elevation angle is not possible for %s. Satellite azimuth and elevation angles are missing.", system)
@@ -298,7 +298,7 @@ def plot_SNR_wrt_elev_dont_use_TEX(analysisResults,GNSS_obs, GNSSsystems, obsCod
             ## -- Subplot 2
             ax[1].grid(True,linewidth=0.3)
             ax[1].set_xlim([0,90])
-            ax[1].set_title("Signal to noise ratio (SNR) as funtion elevation angle for \n Signal: %s (%s)" % (range1_code, system), va='bottom',fontsize=22)
+            ax[1].set_title("Signal to noise ratio (SNR) as funtion of elevation angle for \n Signal: %s (%s)" % (range1_code, system), va='bottom',fontsize=22)
             ax[1].set_ylim(0,np.nanmax(SNR)+10)
             ax[1].set_xlabel('Elevation angle [degree]',fontsize=18,labelpad=10)
             ax[1].set_ylabel('[dB-Hz]',fontsize=18,labelpad=10)
diff --git a/src/gnssmultipath/writeOutputFile.py b/src/gnssmultipath/writeOutputFile.py
index 5251bc0..5bb2368 100644
--- a/src/gnssmultipath/writeOutputFile.py
+++ b/src/gnssmultipath/writeOutputFile.py
@@ -107,8 +107,8 @@ def writeOutputFile(outputFilename, outputDir, analysisResults, includeResultSum
     fid = open(outputFilename, 'w+')
 
     fid.write('GNSS_MultipathAnalysis\n')
-    fid.write('Software version: 1.2.0\n')
-    fid.write('Last software version release: 12/03/2023\n\n')
+    fid.write('Software version: 1.4.3\n')
+    fid.write('Last software version release: 26/11/2023\n\n')
     fid.write('Software developed by Per Helge Aarnes (per.helge.aarnes@gmail.com) \n\n')
     fid.write('RINEX observation filename:\t\t %s\n' % (rinex_obs_filename))
     if sp3_filename is not None:
@@ -373,7 +373,7 @@ def writeOutputFile(outputFilename, outputDir, analysisResults, includeResultSum
                     nSlips30_40msg              = nSlips30_40msg + '%8d|' % (current_code_struct['cycle_slip_distribution']['n_slips_30_40'])
                     nSlips40_50msg              = nSlips40_50msg + '%8d|' % (current_code_struct['cycle_slip_distribution']['n_slips_40_50'])
                     nSlipsOver50msg             = nSlipsOver50msg + '%8d|' % (current_code_struct['cycle_slip_distribution']['n_slips_over50'])
-                    nSlipsNaNmsg                = nSlipsNaNmsg + '%8d|' % (current_code_struct['cycle_slip_distribution']['n_slips_NaN'])
+                    # nSlipsNaNmsg                = nSlipsNaNmsg + '%8d|' % (current_code_struct['cycle_slip_distribution']['n_slips_NaN'])
 
             fid.write(topline + '\n')
             fid.write(headermsg + '\n')
@@ -392,8 +392,8 @@ def writeOutputFile(outputFilename, outputDir, analysisResults, includeResultSum
             fid.write(bottomline + '\n')
             fid.write(nSlipsOver50msg + '\n')
             fid.write(bottomline + '\n')
-            fid.write(nSlipsNaNmsg + '\n')
-            fid.write(bottomline + '\n')
+            # fid.write(nSlipsNaNmsg + '\n')
+            # fid.write(bottomline + '\n')
             fid.write(slipRatiomsg + '\n')
             fid.write(bottomline + '\n')
 


Test patch:
diff --git a/.github/workflows/pypi-package-upload-testpypi.yml b/.github/workflows/pypi-package-upload-testpypi.yml
new file mode 100644
index 0000000..8e6e7b5
--- /dev/null
+++ b/.github/workflows/pypi-package-upload-testpypi.yml
@@ -0,0 +1,41 @@
+name: Upload Python Package to TEST PYPI
+
+on:
+  pull_request:
+    branches:
+      - master
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v2
+
+    - name: Set up Python
+      uses: actions/setup-python@v2
+      with:
+        python-version: 3.8
+
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install --upgrade setuptools wheel build twine
+
+    - name: Get current version from pyproject.toml
+      id: get_version
+      run: echo "::set-output name=version::$(grep -oP '(?<=version = ")[^"]*' pyproject.toml)"
+
+    - name: Increment version
+      id: increment_version
+      run: echo "::set-output name=version::$(python -c "version='${{ steps.get_version.outputs.version }}'.split('.'); version[-1] = str(int(version[-1]) + 1); print('.'.join(version))")"
+
+    - name: Update version in pyproject.toml
+      run: sed -i "s/version = \".*\"/version = \"${{ steps.increment_version.outputs.version }}\"/" pyproject.toml
+
+    - name: Build and upload to Test PyPI (master branch)
+      if: github.ref == 'refs/heads/master'
+      run: |
+        python -m build
+        twine upload --repository-url https://test.pypi.org/legacy/ --username __token__ --password ${{ secrets.TEST_PYPI_TOKEN }} dist/*
diff --git a/tests/analysisResults_NMBUS.pkl b/tests/analysisResults_NMBUS.pkl
index de92237..f6dd00f 100644
Binary files a/tests/analysisResults_NMBUS.pkl and b/tests/analysisResults_NMBUS.pkl differ
diff --git a/tests/analysisResults_OPEC00NOR_S_20220010000_01D_30S_MO.pkl b/tests/analysisResults_OPEC00NOR_S_20220010000_01D_30S_MO.pkl
deleted file mode 100644
index 3f58c21..0000000
Binary files a/tests/analysisResults_OPEC00NOR_S_20220010000_01D_30S_MO.pkl and /dev/null differ
diff --git a/tests/analysisResults_OPEC00NOR_S_20220010000_01D_30S_MO_3.04_croped.pkl b/tests/analysisResults_OPEC00NOR_S_20220010000_01D_30S_MO_3.04_croped.pkl
index 529628c..2540498 100644
Binary files a/tests/analysisResults_OPEC00NOR_S_20220010000_01D_30S_MO_3.04_croped.pkl and b/tests/analysisResults_OPEC00NOR_S_20220010000_01D_30S_MO_3.04_croped.pkl differ
diff --git a/tests/test_GNSS_MultipathAnalysis.py b/tests/test_GNSS_MultipathAnalysis.py
index d9140f0..615a453 100644
--- a/tests/test_GNSS_MultipathAnalysis.py
+++ b/tests/test_GNSS_MultipathAnalysis.py
@@ -1,12 +1,24 @@
 """
-This module is using pytest to run tests on the software. Each pull request to the master branch will
-trigger a workflow that run these tests.
+This module is using pytest to run tests on the software.
+Each pull request to the master branch will trigger
+a workflow that run these tests.
 
 Made by: Per Helge Aarnes
 E-mail: per.helge.aarnes@gmail.com
+
+Example on how to run the tests:
+
+To run all tests:
+    pytest test_GNSS_MultipathAnalysis.py  -vv
+
+To run a specific test:
+    pytest test_GNSS_MultipathAnalysis.py::test_GNSS_MultipathAnalysis_sp3_file  -vv
+
+To run a specific test in debug mode to see where it fails:
+    pytest test_GNSS_MultipathAnalysis.py::test_GNSS_MultipathAnalysis_sp3_file  -vv --pdb  (type "quit" to exit the debug mode)
+
 """
 
-import pickle
 import sys
 import os
 import pytest
@@ -17,6 +29,7 @@
 
 from gnssmultipath import GNSS_MultipathAnalysis
 from gnssmultipath.Geodetic_functions import ECEF2geodb, ECEF2enu, compute_azimut_elev
+from gnssmultipath.PickleHandler import PickleHandler
 
 
 test_data_ECEF2geodb = [
@@ -51,13 +64,11 @@ def test_compute_azimut_elev(X, Y, Z, xm, ym, zm, expected_output):
     result = compute_azimut_elev(X, Y, Z, xm, ym, zm)
     assert result == pytest.approx(expected_output, abs=1e-4)
 
-def read_pickle(file_path):
-    with open(file_path, 'rb') as file:
-        data = pickle.load(file)
-    return data
-
 
 os.chdir(os.path.join(project_path,'src'))
+
+@pytest.mark.filterwarnings("ignore::RuntimeWarning")
+@pytest.mark.filterwarnings("ignore::DeprecationWarning")
 def test_GNSS_MultipathAnalysis_sp3_file():
     """
     Test the results from OPEC2022 and sp3 files
@@ -70,18 +81,15 @@ def test_GNSS_MultipathAnalysis_sp3_file():
                                     plotEstimates=False,
                                     plot_polarplot=False)
 
-    expected_result = read_pickle(expected_res)
+    expected_result = PickleHandler.read_zstd_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']['C']['Azimut'],  result['Sat_position']['C']['Azimut'],decimal=4)
-
+    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']['azimuth'],  result['Sat_position']['G']['azimuth'],decimal=4)
+    assert_almost_equal(expected_result['Sat_position']['E']['azimuth'],  result['Sat_position']['E']['azimuth'],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']['C']['azimuth'],  result['Sat_position']['C']['azimuth'],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)
 
@@ -89,6 +97,12 @@ def test_GNSS_MultipathAnalysis_sp3_file():
 
 
 
+
+
+
+
+@pytest.mark.filterwarnings("ignore::DeprecationWarning")
+@pytest.mark.filterwarnings("ignore::RuntimeWarning")
 def test_GNSS_MultipathAnalysis_broadcast_navfile():
     """
     Test the results from OPEC2022 and use of broadcast ephemerids
@@ -101,32 +115,32 @@ def test_GNSS_MultipathAnalysis_broadcast_navfile():
     result = GNSS_MultipathAnalysis(rinObsFilename=rinObs_file, broadcastNav1=broadNav_file,
                                     plotEstimates=False,
                                     plot_polarplot=False,
-                                    nav_data_rate=0)
+                                    nav_data_rate=120)
 
-    expected_result = read_pickle(expected_res)
+    expected_result = PickleHandler.read_zstd_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['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']['azimuth'],  result['Sat_position']['G']['azimuth'],decimal=4)
+    assert_almost_equal(expected_result['Sat_position']['R']['azimuth'],  result['Sat_position']['R']['azimuth'],decimal=4)
+    assert_almost_equal(expected_result['Sat_position']['E']['azimuth'],  result['Sat_position']['E']['azimuth'],decimal=4)
+    assert_almost_equal(expected_result['Sat_position']['C']['azimuth'],  result['Sat_position']['C']['azimuth'],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)