Skip to content

Commit

Permalink
LEOGPS v1.0 (Stable)
Browse files Browse the repository at this point in the history 
Official Version 1.0 Release!
  • Loading branch information
@sammmlow
sammmlow committed Mar 29, 2021
1 parent 2800f56 commit 0c58ba3
Show file tree
Hide file tree
Showing 113 changed files with 27,570 additions and 133 deletions.
30 changes: 11 additions & 19 deletions README.md
View file
Original file line number Diff line number Diff line change
@@ -1,41 +1,33 @@
![LEOGPS - GPS-Aided Relative Satellite Navigation in Python](https://raw.githubusercontent.com/sammmlow/LEOGPS/master/gui/logo.png)

### What is LEOGPS?
### LEOGPS
LEOGPS is an open-source Python package for absolute and relative satellite navigation, with specific intent for relative navigation between formation flying satellites. It currently supports only observations from the GPS constellation, up to L1/L2 frequency, with input files as RINEX v2.XX. It uses ephemeris and clock files from the University of Bern, CODE.

### How do I use LEOGPS?

First, clone this repository. Then, that's it! The user first runs the application by running 'leogps.py', the main python file. The app interfaces with the inputs configuration file, allowing the user to save or load parameters in the 'config.txt' file, through the GUI.
First, for the full documentation, please refer to it here:

Then, the user simply provides the program with two RINEX (v2.xx) observation files, one for each LEO satellite, by pasting it in the 'inputs' folder. By default, LEOGPS should be provided with two example RINEX files comprising GPS pseudorange and carrier phase measurements from two satellites in formation - GRACE A and B. Ground truths are also provided in the 'outputs' directory.
First, clone this repository. Then, that's it! The user first runs the application by running 'leogps.py', the main python file. The app interfaces with your inputs, and saves it in a configuration file, allowing the user to save or load parameters through the GUI.

Next, the user simply has to run LEOGPS via the GUI, once the user saves the desired configuration options. Keep an internet connection online, as LEOGPS will attempt to download precise ephemerides and clock bias data from the University of Bern's COD FTP (unless you already have those files). Note that in the previous versions, LEOGPS used International GNSS Service (IGS) products. We have now switched to Bern's COD in v0.3.
Then, the user simply provides the program with two RINEX (v2.xx) observation files, one for each LEO satellite, by pasting it in the 'inputs' folder. By default, LEOGPS should be provided with two example RINEX files comprising GPS pseudorange and carrier phase measurements from two satellites in formation - GRACE A and B. Ground truths are also provided in the 'outputs' directory.

The app should look something like the image below.
Next, the user can run LEOGPS. The UI of the app should pop out and look something like below:

![LEOGPS - Graphical User Interface](https://raw.githubusercontent.com/sammmlow/LEOGPS/master/gui/gui.jpg)

Keep an internet connection online, as LEOGPS will attempt to download precise ephemerides and clock bias data from the University of Bern's COD FTP (unless you already have those files). Note that in the previous versions, LEOGPS used International GNSS Service (IGS) products. We have now switched to Bern's COD in v0.3 onwards.

LEOGPS will process the raw GPS measurements to produce a report comprising:

- The absolute positions and absolute velocities of both LEOs.
- Precise (centimeter-level) baseline estimation (relative position vector).
- The single point positions and velocities of both LEOs.
- Precise baseline vectors between the two LEOs.
- Dilution of precision values.

LEOGPS will also (optionally, depending on your choice in the GUI) output plots or reports on the interpolated GPS satellite ephemeris and clock biases. Currently ephemeris errors are still present and we are working to resolve them.

### What is the method of navigation used behind LEOGPS?

Absolute positioning and velocity estimation is performed by trilaterating GPS pseudorange measurements, and Doppler (pseudorange rate) measurements respectively. Single-frequency measurements employ the GRAPHIC linear combination (T. P. Yunck, 1993), while dual frequency measurements employ the ionosphere-free linear combination (Misra & Enge, 2006). If Doppler values are missing in the RINEX observation file, LEOGPS will also attempt to estimate them using a gradient estimation via first order Taylor expansion approximation of the carrier phase values (in v0.2 onwards). However, the velocity estimation results so far are not promising.

The relative navigation between LEO satellites are performed using a double-differencing of carrier phase values, and using the float ambiguities directly. An optional module exists ('ambfix.py') in the codes folder for Dr Peter Teunissen's LAMBDA method for estimating integer ambiguities given the float ambiguities and covariance matrices.

### Why did you create LEOGPS?

LEOGPS is a personal pet project of mine, useful as a modular software used to study relative satellite navigation techniques for distributed space systems with a GPS receiver, such as formation missions like GRACE A and B. LEOGPS' sole purpose is as an adaptable platform for the rapid prototyping of navigation algorithms in Python, or for testing out integer ambiguity resolution techniques. A Pythonic translation of Professor Peter Teunissen's LAMBDA method has also been adapted, in the 'file ambfix.py'

### Requirements for LEOGPS:

Requires Python version 3.6.5 (Anaconda with Spyder recommended).
Tested on Python version 3.6.5 (Anaconda with Spyder).

Core libraries necessary: NumPy (v1.14 and above) and matplotlib

Expand All @@ -45,7 +37,7 @@ Libraries for GUI: PIL, tkinter

The above are all standard Python libraries that should come with a typical Anaconda installation.

This README was last updated on: 12-Jan-2021.
This README was last updated on: 29-March-2021.

### Contact:

Expand Down
47 changes: 44 additions & 3 deletions codes/ambest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
## | | | __ / \ / __| _ | __| ##
## | |__| __ ( ) | (_ | _|__ \ ##
## |____|___ \___/ \___|_| \___/ ##
## v 0.3 (Alpha) ##
## v 1.0 (Stable) ##
## ##
## FILE DESCRIPTION: ##
## ##
Expand Down Expand Up @@ -134,7 +134,7 @@ def ambest( epoch, gps, rx1, rx2, pos1, pos2, inps ):

''' We now proceed to solve the least squares problem. '''

A = getgeom( gps, glist, grefr, grefi, pos1, pos2 ) # Geometry matrix
A = getgeom( gps, glist, grefr, grefi, pos1, pos2, inps ) # Geometry matrix
B = wavelength * (carrDD - zfix) # Carrier phase minus integer ambiguities

# Solve the equation now!
Expand Down Expand Up @@ -200,10 +200,11 @@ def getref( gps, rx1, rx2, pos1, pos2 ):

''' Define a function that finds the DD geometry matrix '''

def getgeom( gps, glist, grefr, grefi, pos1, pos2 ):
def getgeom( gps, glist, grefr, grefi, pos1, pos2, inps ):

DD_geom_mat = np.array([]) # Matrix for DD geometry vector.
leopos = (0.5 * (pos1 + pos2))[:3] # Formation center
erp = inps['earthrotation'] # Offset Earth rotation?

for k in glist:

Expand All @@ -218,6 +219,46 @@ def getgeom( gps, glist, grefr, grefi, pos1, pos2 ):
gps[grefr]['py'],
gps[grefr]['pz']]) * 1000

# Get the GPS satellite velocities
gpskv = np.array([gps[k]['vx'],
gps[k]['vy'],
gps[k]['vz']]) * 1000
gpsrv = np.array([gps[grefr]['vx'],
gps[grefr]['vy'],
gps[grefr]['vz']]) * 1000

# Now, we need to compensate for the GPS satellite positions due
# to Earth's rotation and the movement of the GPS satellites
# during the signal time-of-flight.
if erp == 'True':

# Get the signal time of flights for the k-th GPS SV
sigTOF_k = np.linalg.norm( gpsk - leopos ) / consts.C

# Get the signal time of flights for the k-th GPS SV
sigTOF_r = np.linalg.norm( gpsr - leopos ) / consts.C

# Rotation rate * TOF
Rk = consts.ERR * sigTOF_k
Rr = consts.ERR * sigTOF_r

# Get the Earth rotation direction cosine matrices,
# using small angle approximations for simplicity.
dcm_k = np.array([[ 1.0, Rk, 0.0 ],
[ -1*Rk, 1.0, 0.0 ],
[ 0.0, 0.0, 1.0 ]])
dcm_r = np.array([[ 1.0, Rr, 0.0 ],
[ -1*Rr, 1.0, 0.0 ],
[ 0.0, 0.0, 1.0 ]])

# Rotate the coordinate frame for GPS satellite positions
gpsk = np.matmul(dcm_k, gpsk)
gpsr = np.matmul(dcm_r, gpsr)

# Compensate for GPS satellite motion during signal TOF
gpsk = gpsk - (gpskv * sigTOF_k)
gpsr = gpsr - (gpsrv * sigTOF_r)

# First, obtain the undifferenced relative vector of nth GPS.
ZD_geom1 = gpsk - leopos
ZD_geom1 = ZD_geom1 / np.linalg.norm(ZD_geom1)
Expand Down
5 changes: 3 additions & 2 deletions codes/ambfix.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
## | | | __ / \ / __| _ | __| ##
## | |__| __ ( ) | (_ | _|__ \ ##
## |____|___ \___/ \___|_| \___/ ##
## v 0.3 (Alpha) ##
## v 1.0 (Stable) ##
## ##
## FILE DESCRIPTION: ##
## ##
Expand Down Expand Up @@ -46,7 +46,8 @@
## - Equivalently in Python: for i in range(0,5) => {0,1,2,3,4} ##
## - Indices are thus updated accordingly. ##
## ##
## ORIGINAL AUTHOR: Sandra Verhagen and Bofeng Li ##
## DEVELOPER: Professor Peter Teunissen (TU Delft) ##
## ORIGINAL AUTHOR: Sandra Verhagen and Bofeng Li (TU Delft) ##
## AUTHOR MODIFIED: 26-07-2019, by Samuel Y.W. Low, with permissions. ##
## ##
###############################################################################
Expand Down
2 changes: 1 addition & 1 deletion codes/azimel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
## | | | __ / \ / __| _ | __| ##
## | |__| __ ( ) | (_ | _|__ \ ##
## |____|___ \___/ \___|_| \___/ ##
## v 0.3 (Alpha) ##
## v 1.0 (Stable) ##
## ##
## FILE DESCRIPTION: ##
## ##
Expand Down
2 changes: 1 addition & 1 deletion codes/consts.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
## | | | __ / \ / __| _ | __| ##
## | |__| __ ( ) | (_ | _|__ \ ##
## |____|___ \___/ \___|_| \___/ ##
## v 0.3 (Alpha) ##
## v 1.0 (Stable) ##
## ##
## FILE DESCRIPTION: ##
## ##
Expand Down
4 changes: 2 additions & 2 deletions codes/dopest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
## | | | __ / \ / __| _ | __| ##
## | |__| __ ( ) | (_ | _|__ \ ##
## |____|___ \___/ \___|_| \___/ ##
## v 0.3 (Alpha) ##
## v 1.0 (Stable) ##
## ##
## FILE DESCRIPTION: ##
## ##
Expand Down Expand Up @@ -108,7 +108,7 @@ def dopest(rnxdata, goodsats, tstart, tstop, rnxstep, inps):

# Next, we get an array for an infinitesimal step of
# both time and phase (t + delta-t, and L + delta-L).
delta = 0.00004 # Delta timestep
delta = 0.01 # Delta timestep
N_delta = N + delta # N + delta-N
L_delta = np.polyval(coeff, N_delta) # L + delta-L

Expand Down
14 changes: 8 additions & 6 deletions codes/einstn.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
## | | | __ / \ / __| _ | __| ##
## | |__| __ ( ) | (_ | _|__ \ ##
## |____|___ \___/ \___|_| \___/ ##
## v 0.3 (Alpha) ##
## v 1.0 (Stable) ##
## ##
## FILE DESCRIPTION: ##
## ##
Expand Down Expand Up @@ -51,16 +51,17 @@
# IMPORT LOCAL LIBRARIES
from codes import consts

gm = consts.GM # GRAVITY CONSTANT*EARTH MASS
c = consts.C # VELOCITY OF LIGHT
gm = consts.GM # GRAVITY CONSTANT * EARTH MASS
c = consts.C # VELOCITY OF LIGHT (M/S)
err = consts.ERR # # EARTH INERTIAL ROTATION RATE (RAD/S)

''' CLOCK ADVANCE EFFECT '''

def clockadv(gpspos,gpsvel):

# This effect is more dominant than the Shapiro Delay effect.

return -2*np.dot(gpspos,gpsvel)/c # Units in meters
return -2*np.dot(gpspos,gpsvel)/c # Units in meters, thus c is not squared

''' SHAPIRO DELAY EFFECT '''

Expand All @@ -74,6 +75,7 @@ def shapiro(leopos,gpspos):
gps = np.linalg.norm(gpspos)
rel = np.linalg.norm(gpspos-leopos)
coeff = 2*gm/(c**2)
ratio = ( gps + leo + rel / gps + leo - rel )
ratio = ( gps + leo + rel ) / ( gps + leo - rel )
shapiro_path_correction = coeff * np.log( ratio )

return coeff*math.log(ratio) # Units in meters
return shapiro_path_correction # Units in meters
44 changes: 38 additions & 6 deletions codes/gpsxtr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
## | | | __ / \ / __| _ | __| ##
## | |__| __ ( ) | (_ | _|__ \ ##
## |____|___ \___/ \___|_| \___/ ##
## v 0.3 (Alpha) ##
## v 1.0 (Stable) ##
## ##
## FILE DESCRIPTION: ##
## ##
Expand Down Expand Up @@ -51,6 +51,18 @@
import numpy as np
import urllib.request





# DELETE AFTER
import matplotlib.pyplot as plt






# IMPORT LOCAL LIBRARIES
from codes import pubplt

Expand Down Expand Up @@ -246,7 +258,7 @@ def gpsxtr(inps, tstart, tstop, tstep):

# Sort the list of good satellites found.
goodsats.sort() # Sort in ascending order

###########################################################################
###########################################################################

Expand Down Expand Up @@ -408,8 +420,8 @@ def gpsxtr(inps, tstart, tstop, tstep):
# period for interpolation. See research paper "Polynomial interpolation
# of GPS satellite coordinates" by Milan Horemuz (Feb 2006).
# Thus, we use a sliding window interpolation, with fit interval of 4h,
# and a validity window in the middle of 15 minutes only, leaving the rest
# of the 4h as interpolation buffer to prevent Runge's phenomenon.
# and a validity window in the middle of 15 minutes only, leaving the
# rest of the 4h as interpolation buffer to prevent Runge's phenomenon.

# Our first task is to generate the IGS ephemeris time axis.
gpsephm_epoch_dt = [] # List of datetime objects.
Expand All @@ -422,7 +434,7 @@ def gpsxtr(inps, tstart, tstop, tstep):
gpsephm_epoch_dt.append(gpsephm_dt)
gpsephm_epoch_ss.append(gpsephm_ss)
gpsephm_dt = gpsephm_dt + gps_tstep
gpsephm_ss = gpsephm_ss + 900
gpsephm_ss = gpsephm_ss + gps_tstep.seconds

# We also initialise the starting time in GPS ephemeris interpolation.
# This block of code is meant to create a time axis in integer seconds.
Expand Down Expand Up @@ -505,7 +517,27 @@ def gpsxtr(inps, tstart, tstop, tstep):
break # End the while loop if we are.
else:
windex += 1 # Update the sliding window filter by 2 hours


# # For debugging
# xt = []
# xx= []
# yp = []
# ypi = []
# count = t_offset_eph
# for t in gpsdata.keys():
# xt.append(count)
# ypi.append(gpsdata[t][30]['px'])
# count += tstep.seconds
# count = 0
# for t1 in gpsephm_epoch_dt:
# xx.append(count)
# yp.append(gpsephm[30][t1]['px'])
# count += gps_tstep.seconds
# plt.figure(1)
# plt.plot(xt,ypi)
# plt.scatter(xx,yp)
# plt.show()

###########################################################################
###########################################################################

Expand Down
2 changes: 1 addition & 1 deletion codes/inpxtr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
## | | | __ / \ / __| _ | __| ##
## | |__| __ ( ) | (_ | _|__ \ ##
## |____|___ \___/ \___|_| \___/ ##
## v 0.3 (Alpha) ##
## v 1.0 (Stable) ##
## ##
## FILE DESCRIPTION: ##
## ##
Expand Down
12 changes: 9 additions & 3 deletions codes/leogui.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
## | | | __ / \ / __| _ | __| ##
## | |__| __ ( ) | (_ | _|__ \ ##
## |____|___ \___/ \___|_| \___/ ##
## v 0.3 (Alpha) ##
## v 1.0 (Stable) ##
## ##
## FILE DESCRIPTION: ##
## ##
Expand Down Expand Up @@ -35,7 +35,7 @@ class run_gui:
def __init__(self, master):

# Create the main frame and window.
master.title('LEOGPS v0.3 Alpha')
master.title('LEOGPS v1.0')
master.geometry('1024x896')

# Initialise the basic text labels (found in the configuration file):
Expand Down Expand Up @@ -646,6 +646,12 @@ def cfg_W(self):

def run(self):

leorun.run()
try:
self.cfg_W()
leorun.run()
except Exception as excpt:
print('Error in running!')
print(excpt)
pass

return None
Loading

0 comments on commit 0c58ba3

Please sign in to comment.