This Jupyter notebook simulates and compares the DOP depending on the user latitude for the Galileo and GPS constellations. Both the reference constellations (with and without spares/additional satellites) and the real constellations taken from the almanac at 2020-01-01 are used.
In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from astropy.time import Time
import astropy.units as u
import astropy.constants as const
from astropy.coordinates import Angle, EarthLocation, ITRS
import cartopy.crs as ccrs
from multiprocessing import Pool
from functools import partial
In [2]:
def era(t):
# formula from https://en.wikipedia.org/wiki/Sidereal_time#Earth_Rotation_Angle
# direct computation in ITRS using LAAN has been checked to yield
# similar results to computation in GCRS and conversion to ITRS with Astropy
tU = t - Time(2451545, format = 'jd')
return Angle(2*np.pi*(0.7790572732640+1.00273781191135448*tU/u.day)*u.rad).wrap_at('180d')
class Constellation:
def __init__(self, a, e, i, raan0, omega0, M0, t0, draan = 0*u.deg/u.day, domega = 0*u.deg/u.day, dM = None):
self.a = a
self.e = e
self.i = i
self.raan0 = raan0
self.omega0 = omega0
self.M0 = M0
self.t0 = t0
self.draan = draan
self.domega = domega
self.dM = dM
def xyz(self, t):
t = t.reshape((-1,1))
# use mean anomaly rate if provided, else derive it from Keplerian orbit
n = self.dM if self.dM is not None else np.sqrt(const.GM_earth/self.a**3) * u.rad
tk = (t - self.t0)
M = (self.M0 + n*tk).wrap_at('180d')
# series expansion in e for nu
nu = M + (self.e * np.sin(M) + self.e**2*5/4*np.sin(2*M) - self.e**3*(1/4*np.sin(M) \
- 13/12*np.sin(3*M)))* u.rad if self.e else M
omega = self.omega0 + self.domega * tk
uu = nu + omega
r = self.a * (1 - self.e**2) / (1 + self.e*np.cos(nu)) if self.e else self.a
x_ = r * np.cos(uu)
y_ = r * np.sin(uu)
raan = self.raan0 + self.draan * tk
laan = raan - era(t)
x = x_ * np.cos(laan) - y_ * np.cos(self.i) * np.sin(laan)
y = x_ * np.sin(laan) + y_ * np.cos(self.i) * np.cos(laan)
z = y_ * np.sin(self.i)
return ITRS(x = x, y = y, z = z, obstime = t, representation_type = 'cartesian')
@property
def size(self):
return self.M0.size
Galileo reference constellation information taken from GSC Orbital and technical parameters. Note that there is a typo and slots B2 and B7 are swapped. This constellation is a Walker 24/3/1 constellation.
In [3]:
# information taken from https://www.gsc-europa.eu/system-service-status/orbital-and-technical-parameters
galileo_data = """GSAT0101 E11 B05 29599.8 0.0 56.0 77.632 0.0 15.153
GSAT0102 E12 B06 29599.8 0.0 56.0 77.632 0.0 60.153
GSAT0103 E19 C04 29599.8 0.0 56.0 197.632 0.0 345.153
GSAT0104 E20 C05 29599.8 0.0 56.0 197.632 0.0 30.153
GSAT0203 E26 B08 29599.8 0.0 56.0 77.632 0.0 150.153
GSAT0204 E22 B03 29599.8 0.0 56.0 77.632 0.0 285.153
GSAT0205 E24 A08 29599.8 0.0 56.0 317.632 0.0 135.153
GSAT0206 E30 A05 29599.8 0.0 56.0 317.632 0.0 0.153
GSAT0208 E08 C07 29599.8 0.0 56.0 197.632 0.0 120.153
GSAT0209 E09 C02 29599.8 0.0 56.0 197.632 0.0 255.153
GSAT0210 E01 A02 29599.8 0.0 56.0 317.632 0.0 225.153
GSAT0211 E02 A06 29599.8 0.0 56.0 317.632 0.0 45.153
GSAT0207 E07 C06 29599.8 0.0 56.0 197.632 0.0 75.153
GSAT0212 E03 C08 29599.8 0.0 56.0 197.632 0.0 165.153
GSAT0213 E04 C03 29599.8 0.0 56.0 197.632 0.0 300.153
GSAT0214 E05 C01 29599.8 0.0 56.0 197.632 0.0 210.153
GSAT0215 E21 A03 29599.8 0.0 56.0 317.632 0.0 270.153
GSAT0216 E25 A07 29599.8 0.0 56.0 317.632 0.0 90.153
GSAT0217 E27 A04 29599.8 0.0 56.0 317.632 0.0 315.153
GSAT0218 E31 A01 29599.8 0.0 56.0 317.632 0.0 180.153
GSAT0219 E36 B04 29599.8 0.0 56.0 77.632 0.0 330.153
GSAT0220 E13 B01 29599.8 0.0 56.0 77.632 0.0 195.153
GSAT0221 E15 B02 29599.8 0.0 56.0 77.632 0.0 105.153
GSAT0222 E33 B07 29599.8 0.0 56.0 77.632 0.0 240.153"""
The positions of the 6 future spares satellites haven't been defined anywhere. We use the example given by Teunissen etl. to position the spares.
In [4]:
raan = [float(l.split()[-3]) for l in galileo_data.split('\n')]
galileo_raan = Angle(np.array([raan])*u.deg)
M = [float(l.split()[-1]) for l in galileo_data.split('\n')]
galileo_M = Angle(np.array([M])*u.deg)
galileo_names = [l.split()[0] + f' ({l.split()[1]})' for l in galileo_data.split('\n')]
galileo_slots = [l.split()[2] for l in galileo_data.split('\n')]
# constellation with 6 spares. the location of the spares is indicative, as it hasn't been published
# spares between A2,A3, A6,A7, B1,B2, B5,B6, C8,C1, C4,C5
# note that B2 and B7 are swapped in galileo_data listing above
raan_spare = np.repeat([317.632, 77.632, 197.632], 2)
M_spare = [(225.153+270.153)/2, (45.153+90.153)/2, (195.153+240.153)/2, (15.153+60.153)/2,\
(165.153+210.153)/2, (345.153-360+30.153)/2]
galileo_30_raan = Angle(np.concatenate((raan, raan_spare)).reshape((1,-1))*u.deg)
galileo_30_M = Angle(np.concatenate((M, M_spare)).reshape((1,-1))*u.deg)
galileo_a = 29599.8 * u.km
galileo_i = Angle(56.0 * u.deg)
galileo_t0 = Time('2016-01-01T00:00:00')
galileo_draan = -0.02764398 * u.deg / u.day
galileo_dM = 613.72253566 * u.deg / u.day
We load up the almanac for the first week of 2020. The almanac contains the LAAN at the start of the Galileo week (2019-12-29T00:00:00 GST = 2019-12-28T23:59:42 UTC), which must be converted to RAAN using the Earth rotation angle. The t0a of the almanac is 463200 which gives 2020-01-03T08:40:00 GST = 2020-01-03T08:39:42 UTC.
The almanac contains 22 satellites, which are all the active satellites.
In [5]:
def get_xml_field(lines, field):
return [l.strip().lstrip(f'<{field}>').rstrip(f'</{field}>') for l in lines if l.strip().startswith(f'<{field}>')]
with open('2020-01-03.xml') as f:
galileo_alm = f.readlines()
galileo_alm_prns = [int(n) for n in get_xml_field(galileo_alm, 'SVID')]
galileo_alm_laan = Angle(np.array([[float(f) for f in get_xml_field(galileo_alm, 'omega0')]])*180*u.deg)
galileo_alm_raan = galileo_alm_laan + era(Time('2019-12-28T23:59:42'))
galileo_alm_omega = Angle(np.array([[float(f) for f in get_xml_field(galileo_alm, 'w')]])*180*u.deg)
galileo_alm_M = Angle(np.array([[float(f) for f in get_xml_field(galileo_alm, 'm0')]])*180*u.deg)
galileo_alm_t0 = Time('2020-01-03T08:39:42')
We plot and compare the almanac constellation with the reference constellation. To do so, we move the almanac constellation to the timestamp of the reference constellation, updating the RAAN and mean anomalies at their nominal rates.
In [6]:
alpha = (galileo_alm_t0 - galileo_t0) * galileo_draan
beta = (galileo_alm_t0 - galileo_t0) * galileo_dM
plt.figure(figsize = (10,6))
plt.plot(galileo_30_raan.ravel(), galileo_30_M.ravel(), 'o', color = 'C2', label = 'Reference spares')
plt.plot((galileo_alm_raan - alpha).wrap_at('360d').ravel(),\
(galileo_alm_omega + galileo_alm_M - beta).wrap_at('360d').ravel(), 'o', color = 'C0',\
label = 'Almanac 2020-01-01');
plt.plot(galileo_raan.ravel(), galileo_M.ravel(), 'o', color = 'C1', label = 'Reference constellation')
for j in range(galileo_alm_raan.size):
plt.text((galileo_alm_raan - alpha).wrap_at('360d')[0,j].deg,\
(galileo_alm_omega + galileo_alm_M - beta).wrap_at('360d')[0,j].deg,\
f'E{galileo_alm_prns[j]:02d}')
for j in range(galileo_raan.size):
plt.text(galileo_raan[0,j].deg,\
galileo_M[0,j].deg, galileo_slots[j] + ' ' + galileo_names[j])
plt.xlim((0,360))
plt.ylim((0,360))
plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.12),
fancybox=True, shadow=True, ncol=5)
plt.title(f'Galileo constellation at {galileo_t0.datetime}')
plt.xlabel('RAAN (deg)')
plt.ylabel('Mean anomaly (deg)');
In [7]:
# Galileo constellation
galileo_constellation = Constellation(a = galileo_a,\
e = 0,\
i = galileo_i,\
raan0 = galileo_raan,\
omega0 = Angle(0.0 * u.deg),\
M0 = galileo_M,\
t0 = galileo_t0,\
draan = galileo_draan,\
dM = galileo_dM)
galileo_30_constellation = Constellation(a = galileo_a,\
e = 0,\
i = galileo_i,\
raan0 = galileo_30_raan,\
omega0 = Angle(0.0 * u.deg),\
M0 = galileo_30_M,\
t0 = galileo_t0,\
draan = galileo_draan,\
dM = galileo_dM)
galileo_alm_constellation = Constellation(a = galileo_a,\
e = 0,\
i = galileo_i,\
raan0 = galileo_alm_raan,\
omega0 = galileo_alm_omega,\
M0 = galileo_alm_M,\
t0 = galileo_alm_t0,\
draan = galileo_draan,\
dM = galileo_dM)
The reference constellation for GPS is taken from the 2008 SPS performance standard. Two constellations are defined: a 24 satellite constellation and an expandable constellation where 3 of the slots in the 24 satellite constellation can be split in two to accomodate up to 27 satellites. We consider the 24 satellite constellation and the fully split 27 satellite constellation only.
Update: In April 2020 the 2020 SPS performance standard defined a 30 satellite expanded constellation, here 6 slots of the 24 satellite base constellation are split into two. Moreover, the RAANs of the constellation have precessed 16.003 degrees. The new document also gives updated values for the semi-major axis and the drift of the RAAN. This new version of the notebook uses the constellations as defined in the 2020 SPS performance standard.
In [8]:
# values taken from https://www.gps.gov/technical/ps/2020-SPS-performance-standard.pdf
gps_raan = Angle(np.repeat([288.85, 348.85, 48.85, 108.85, 168.85, 228.85], 4).reshape((1,-1)) * u.deg)
gps_M = Angle(np.array([[239.54, 133.20, 343.09, 13.22,
52.37, 144.75, 281.39, 175.79,
83.29, 343.21, 311.08, 212.97,
106.64, 236.86, 6.57, 138.77,
168.46, 274.01, 37.48, 305.10,
210.30, 316.64, 76.62, 106.76]]) * u.deg)
gps_raan_exp30 = Angle(np.concatenate((
np.repeat([288.85, 348.85, 48.85, 108.85, 168.85, 228.85], 4),
[348.85, 108.85, 228.85, 288.85, 48.85, 168.85])).reshape((1,-1)) * u.deg)
gps_M_exp30 = Angle(np.array([[239.54, 144.40, 343.09, 13.22,
66.33, 144.75, 281.39, 175.79,
83.29, 343.21, 311.08, 220.43,
106.64, 254.09, 6.57, 138.77,
168.46, 274.01, 52.08, 305.10,
210.30, 331.87, 76.62, 106.76,
37.77, 229.39, 305.43,
117.96, 195.73, 23.52]]) * u.deg)
gps_a = 26559.8 * u.km
gps_i = Angle(55.0 * u.deg)
gps_t0 = Time('2016-12-31T23:59:43')
gps_draan = -0.0402 * u.deg / u.day
We load up the almanac for the first week of 2020. Note that the parameter "Right Ascen" is actually the LAAN at the start of the reference week (2019-12-29T00:00:00 GPST = 2019-12-28T23:59:42 UTC). The t0a is 503808, which gives 2020-01-03T19:56:48 GPST = 2020-01-03T19:56:30 UTC.
The almanac contains 31 satellites, including G04 (the new block III satellite), which was unhealthy at the beginning of the year but became usable after a few days, so it is also considered here.
In [9]:
with open('001.ALM') as f:
gps_alm = f.readlines()
gps_alm_prns = [int(l.split(':')[1]) for l in gps_alm if l.startswith('ID')]
gps_alm_laan = Angle((np.array([[float(l.split(':')[1]) for l in gps_alm if l.startswith('Right Ascen')]])*u.rad).to(u.deg)).wrap_at('360d')
gps_alm_raan = gps_alm_laan + era(Time('2019-12-28T23:59:42'))
gps_alm_omega = Angle((np.array([[float(l.split(':')[1]) for l in gps_alm if l.startswith('Argument of Perigee')]])*u.rad).to(u.deg)).wrap_at('360d')
gps_alm_M = Angle((np.array([[float(l.split(':')[1]) for l in gps_alm if l.startswith('Mean Anom')]])*u.rad).to(u.deg)).wrap_at('360d')
gps_alm_t0 = Time('2020-01-03T19:56:30')
gps_alm_slots = ['D2', 'D1', 'E1', 'A6', 'E3', 'D4', 'A4', 'C3', 'F3', 'E2', 'D5', 'B4', 'F6', 'F5', 'F2', 'B1', 'C4', 'C5', 'E4', 'D3', 'E6', 'F4', 'A1', 'B2', 'B5', 'C2', 'B3', 'C1', 'A3', 'A2', 'F1']
gps_slots = [plane + str(n) for plane in 'ABCDEF' for n in range(1,5)]
We plot and compare the almanac constellation with the reference constellation. To do so, we move the almanac constellation to the timestamp of the reference constellation, updating the RAAN and mean anomalies at their nominal rates. The nominal mean anomaly rate (mean motion) is taken to be 2 revolutions per sidereal day.
Update: The 2020 SPS gives a slightly different orbital radius that causes a slightly different dM. This is used below to propagate the almanac constellation to the reference constellation epoch, obtain good agreement. The new. As seen below, the old orbital radius of 26559.7km made a full revolution take 5.02 seconds less than a sidereal day, while the new orbital radiuse of 26559.8km makes a full revolution take 4.77 seconds less than a sidereal day. The difference is small, but noticeable when propagated over several years.
In [10]:
gps_dM = (np.sqrt(const.GM_earth/(26559.7*u.km)**3) * u.rad).to(u.rad/u.s)
((360*u.deg/gps_dM) - 0.5*u.sday).to(u.s)
Out[10]:
In [11]:
gps_dM = (np.sqrt(const.GM_earth/gps_a**3) * u.rad).to(u.rad/u.s)
((360*u.deg/gps_dM) - 0.5*u.sday).to(u.s)
Out[11]:
In [12]:
((1/((360*u.deg/(0.5*u.sday)/u.rad)**2/const.GM_earth)).to(u.km**3))**(1/3)
Out[12]:
In [13]:
alpha = (gps_alm_t0 - gps_t0) * gps_draan
# This was the old beta of 2 revolutions / sday
# beta = (gps_alm_t0 - gps_t0) * 2*360*u.deg/u.sday
# This new beta is not exactly 2 revolutions per sday
beta = (gps_alm_t0 - gps_t0) * gps_dM
plt.figure(figsize = (10,6))
plt.plot((gps_alm_raan - alpha).wrap_at('360d').ravel(),\
(gps_alm_M + gps_alm_omega - beta).wrap_at('360d').ravel(), 'o', color = 'C0', label = 'Almanac 2020-01-01')
plt.plot(gps_raan.ravel(), gps_M.ravel(), 'o', color = 'C2', label = 'Reference constellation')
plt.plot(gps_raan_exp30.ravel(), gps_M_exp30.ravel(), 'o', color = 'C1', label = 'Expandable constellation')
for j in range(gps_raan.size):
plt.text(gps_raan[0,j].deg, gps_M[0,j].deg, gps_slots[j])
for j in range(gps_alm_raan.size):
plt.annotate(f'{gps_alm_slots[j]} G{gps_alm_prns[j]:02d}',
xy = ((gps_alm_raan - alpha).wrap_at('360d')[0,j].deg,
(gps_alm_M + gps_alm_omega - beta).wrap_at('360d')[0,j].deg),
xytext = ((gps_alm_raan - alpha).wrap_at('360d')[0,j].deg - 30,
(gps_alm_M + gps_alm_omega - beta).wrap_at('360d')[0,j].deg))
plt.xlim((0,360))
plt.ylim((0,360))
# Put a legend below current axis
plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.12),
fancybox=True, shadow=True, ncol=5)
plt.title(f'GPS constellation at {gps_t0.datetime}')
plt.xlabel('RAAN (deg)')
plt.ylabel('Mean anomaly (deg)');
In [14]:
gps_constellation = Constellation(a = gps_a,\
e = 0.0,\
i = gps_i,\
raan0 = gps_raan,\
omega0 = Angle(0.0 * u.deg),\
M0 = gps_M,\
t0 = gps_t0,\
draan = gps_draan)
gps_exp30_constellation = Constellation(a = gps_a,\
e = 0.0,\
i = gps_i,\
raan0 = gps_raan_exp30,\
omega0 = Angle(0.0 * u.deg),\
M0 = gps_M_exp30,\
t0 = gps_t0,\
draan = gps_draan)
gps_alm_constellation = Constellation(a = 26559.7 * u.km,\
e = 0.0,\
i = gps_i,\
raan0 = gps_alm_raan,\
omega0 = gps_alm_omega,\
M0 = gps_alm_M,\
t0 = gps_alm_t0,\
draan = gps_draan)
Below we plot the position of the reference and almanac constellations for the current time. The almanac constellations can be compared to real time status in galmon.eu to validate the orbit propagation and constellation parameters.
In [15]:
def plot_constellation(constellation, names, t):
plt.figure(figsize = (14,8))
ax = plt.axes(projection=ccrs.PlateCarree())
xyz = constellation.xyz(t)
plt.plot(xyz.earth_location.geodetic.lon, xyz.earth_location.geodetic.lat, 'o', color = 'C0')
for j,name in enumerate(names):
plt.text(xyz[0,j].earth_location.geodetic.lon.deg, xyz[0,j].earth_location.geodetic.lat.deg, name)
ax.coastlines()
ax.set_global()
In [16]:
plot_constellation(galileo_constellation, galileo_names, Time.now())
In [17]:
plot_constellation(galileo_alm_constellation, [f'E{n:02d}' for n in galileo_alm_prns], Time.now())
In [18]:
plot_constellation(gps_constellation, gps_slots, Time.now())
In [19]:
plot_constellation(gps_alm_constellation, [f'G{n:02d} {slot}' for n,slot in zip(gps_alm_prns, gps_alm_slots)],\
Time.now())
We only consider satellites above an elevation mask of 5 degrees
In [20]:
elevation_mask = 5 * u.deg
In [21]:
# This is a reference implementation, used to validate the optimized implementation given below
def compute_dops(xyz, observer):
xyz = np.einsum('jki->ijk', xyz)
uvw = xyz - np.array(observer.get_itrs().cartesian.xyz.to(u.km))[:,np.newaxis,np.newaxis]
# rotation matrix taken from https://scivision.github.io/pymap3d/ecef.html uvw2enu()
uvw2enu = np.array([[-np.sin(observer.lon), np.cos(observer.lon), 0],\
[-np.sin(observer.lat)*np.cos(observer.lon), -np.sin(observer.lat)*np.sin(observer.lon), np.cos(observer.lat)],\
[np.cos(observer.lat)*np.cos(observer.lon), np.cos(observer.lat)*np.sin(observer.lon), np.sin(observer.lat)]])
enu = np.einsum('ijk,li->jkl', uvw, uvw2enu)
ej = enu/np.sqrt(np.sum(enu**2, axis = 2))[...,np.newaxis]
visible = ej[:,:,2] >= np.sin(elevation_mask)
J0 = np.empty((enu.shape[0], enu.shape[1], 4))
J0[:,:,:3] = -ej
J0[:,:,3] = 1
J0[~visible, :] = 0
DOPs = np.linalg.inv(np.einsum('ijl,ijk->ilk', J0, J0))[:,np.arange(4),np.arange(4)]
return DOPs
def alloc_J0(time_steps, nsats):
return np.empty((time_steps, nsats, 4))
# optimized implementation using Cholesky, an external J0 array and some code rewrite
# J0 is passed to prevent realocating all the time
def compute_dops_cholesky(xyz, observer, J0):
uvw = xyz - np.array(observer.get_itrs().cartesian.xyz.to(u.km))[np.newaxis,np.newaxis,:]
# rotation matrix taken from https://scivision.github.io/pymap3d/ecef.html uvw2enu()
#return None
lon = observer.lon
lat = observer.lat
slon = np.sin(lon)
slat = np.sin(lat)
clon = np.cos(lon)
clat = np.cos(lat)
uvw2enu = np.array([[-slon, clon, 0],\
[-slat*clon, -slat*slon, clat],\
[clat*clon, clat*slon, slat]])
enu = np.einsum('ijk,lk->ijl', uvw, -uvw2enu) # really -enu
J0[:,:,:3] = enu/np.sqrt(np.einsum('ijk,ijk->ij', enu, enu))[...,np.newaxis]
visible = J0[:,:,2] < -np.sin(elevation_mask) # note we have -ej in J0
J0[:,:,3] = 1
J0[~visible, :] = 0
# Cholesky decomposition L*L = J0*J0
L = np.linalg.cholesky(np.einsum('ijl,ijk->ilk', J0, J0))
DOPs = np.empty((J0.shape[0],4))
for j in range(4):
en = np.zeros((1,4))
en[0,j] = 1
v = np.linalg.solve(L, en)
DOPs[:, j] = np.einsum('ij,ij->i', v.conj(), v)
return DOPs
The simulation time is the constellation repetition rate (10 sidereal days for Galileo, 1 sidereal day for GPS), in steps of approximately 1 minute.
In [22]:
repetition_gal = 10 * u.sday
t_start = Time('2020-01-01T00:00:00')
t_gal = t_start + np.linspace(0, repetition_gal, int(np.ceil(repetition_gal.to(u.min)/u.min)))
repetition_gps = 1 * u.sday
t_gps = t_start + np.linspace(0, repetition_gps, int(np.ceil(repetition_gps.to(u.min)/u.min)))
In [23]:
dop_types = ['GDOP', 'PDOP', 'HDOP', 'VDOP']
def dop(dops, dop_type):
if dop_type == 'GDOP':
return np.sqrt(np.sum(dops, axis = 1))
if dop_type == 'PDOP':
return np.sqrt(np.sum(dops[:,:3], axis = 1))
if dop_type == 'HDOP':
return np.sqrt(np.sum(dops[:, :2], axis = 1))
if dop_type == 'VDOP':
return np.sqrt(dops[:, 2])
raise ValueError('Invalid DOP type')
def process_dops(dops):
return [np.average(dop(dops, dop_type)) for dop_type in dop_types] + \
[np.max(dop(dops, dop_type)) for dop_type in dop_types]
def dops_vs_lat(xyz, J0, lats, lon):
return np.array([process_dops(compute_dops_cholesky(xyz, EarthLocation.from_geodetic(lon = lon, lat = lat), J0))\
for lat in lats])
def dops_vs_geo(constellation, t, lats, lons):
xyz = np.empty((t.size, constellation.size, 3))
xyz[:] = np.einsum('ijk->jki', constellation.xyz(t).cartesian.xyz.to(u.km))
J0 = alloc_J0(t.size, constellation.size)
with Pool() as p:
return np.array(p.map(partial(dops_vs_lat, xyz, J0, lats), lons))
We simulate, computing DOP across the globe, with a 1º step in both latitude and longitude.
In [24]:
lats = np.arange(-90,91,1)
lons = np.arange(-180,180,1)
Uncomment to re-run the simulation. This takes 4 hours to run in my laptop.
In [25]:
%%time
#print('Simulating Galileo reference (24)...')
#galileo_dops = dops_vs_geo(galileo_constellation, t_gal, lats, lons)
#np.save('galileo_dops.npy', galileo_dops)
#print('Simulating Galileo reference + spares (30)...')
#galileo_30_dops = dops_vs_geo(galileo_30_constellation, t_gal, lats, lons)
#np.save('galileo_30_dops.npy', galileo_30_dops)
#print('Simulating Galileo almanac (22)...')
#galileo_alm_dops = dops_vs_geo(galileo_alm_constellation, t_gal, lats, lons)
#np.save('galileo_alm_dops.npy', galileo_alm_dops)
#print('Simulating GPS reference (24)...')
#gps_dops = dops_vs_geo(gps_constellation, t_gps, lats, lons)
#np.save('gps_dops.npy', gps_dops)
#print('Simulating GPS expandable (30)...')
#gps_exp30_dops = dops_vs_geo(gps_exp30_constellation, t_gps, lats, lons)
#np.save('gps_exp30_dops.npy', gps_exp30_dops)
#print('Simulating GPS almanac (31)...')
#gps_alm_dops = dops_vs_geo(gps_alm_constellation, t_gps, lats, lons)
#np.save('gps_alm_dops.npy', gps_alm_dops)
In [26]:
galileo_dops = np.load('galileo_dops.npy')
galileo_30_dops = np.load('galileo_30_dops.npy')
galileo_alm_dops = np.load('galileo_alm_dops.npy')
gps_dops = np.load('gps_dops.npy')
gps_exp30_dops = np.load('gps_exp30_dops.npy')
gps_alm_dops = np.load('gps_alm_dops.npy')
In [27]:
dop_list = [galileo_dops, galileo_30_dops, galileo_alm_dops, gps_dops, gps_exp30_dops, gps_alm_dops]
names = ['Galileo (24)', 'Galileo + spares (30)', 'Galileo almanac (22)', 'GPS (24)', 'GPS expanded (30)', 'GPS almanac (31)']
layout = [1,3,5,2,4,6]
In [28]:
vmins = [1.6, 1.45, 0.7, 1.2, 1.95, 1.8, 0.8, 1.6]
vmaxs = [2.5, 2.2, 1.2, 2, 7, 6, 2.5, 5.8]
In [29]:
for i in range(8):
fig = plt.figure(figsize = (14,12), facecolor = 'w')
for j, constellation in enumerate(names):
ax = plt.subplot(3, 2, layout[j], projection = ccrs.PlateCarree())
gridlon, gridlat = np.meshgrid(lons, lats)
img = ax.pcolormesh(gridlon, gridlat, dop_list[j][:,:,i].T, vmin = vmins[i], vmax = vmaxs[i])
ax.coastlines()
ax.set_title(constellation)
cax = plt.axes([0.15, 0.11, 0.72, 0.02])
fig.colorbar(img, cax = cax, orientation = 'horizontal', ax = ax)
plt.subplots_adjust(wspace = 0.02, hspace = 0.03)
plot_type = 'Average' if i < 4 else 'Worst'
fig.suptitle(f'{plot_type} {dop_types[i%4]}', y = 0.9, fontsize = 14)
In [30]:
colors = ['C0', 'C0', 'C0', 'C1', 'C1', 'C1']
styles = [':', '--', '-', ':', '--', '-']
In [31]:
for i in range(8):
plt.figure(figsize = (10,6), facecolor = 'w')
for j, constellation in enumerate(names):
combine = np.average if i < 4 else np.max
plt.plot(lats[lats >= 0], combine(dop_list[j][:,lats>=0, i], axis = 0), styles[j], color = colors[j])
plt.legend(names)
plt.grid('both')
plt.xticks(np.arange(0,90+10,10))
plt.xlabel('Latitude (deg)')
plt.ylabel('DOP')
plot_type = 'Average' if i < 4 else 'Worst'
plt.title(f'{plot_type} {dop_types[i%4]}')
if i >= 4:
plt.ylim(top = vmaxs[i])