66_python_gps_91485
eceive
#!/us
in/python3.4
#Receiver Program
import sys
import math
import numpy
data_file = './data.dat'
class Sat_trans:
def __init__(self, num):
self.num = num
##########Helper Functions##########
#Converts geo coordinates to cartesian at time tv
def geo_to_cart(tv, lat_d, lat_m, lat_s, NS, lon_d, lon_m, lon_s, EW, h):
psi = NS*deg_to_rad(lat_d, lat_m, lat_s)
lam = EW*deg_to_rad(lon_d, lon_m, lon_s) + 2*pi*tv/s
r = R+h
x = r * math.cos(psi) * math.cos(lam)
y = r * math.cos(psi) * math.sin(lam)
z = r * math.sin(psi)
return [x,y,z]
#Converts cartesian coordinates to geo at time tv
def cart_to_geo(tv, xv):
x = xv[0]
y = xv[1]
z = xv[2]
h = numpy.linalg.norm(xv) - R
if (x**2 + y**2) is not 0:
psi = math.atan( z/((x**2 + y**2)**(.5)) )
elif z > 0:
psi = pi/2
else:
psi = -pi/2
if x == 0:
if y > 0:
lam = pi/2
else:
lam = -pi/2
elif x > 0 and y > 0:
lam = math.atan(y/x)
elif x < 0:
lam = pi + math.atan(y/x)
else:
lam = 2*pi + math.atan(y/x)
lam = lam - 2*pi*tv/s
while lam < -1*pi:
lam = lam + 2*pi
while lam > pi:
lam = lam - 2*pi
EW = 1
if lam < 0:
EW = -1
NS = 1
if psi < 0:
NS = -1
(lat_d, lat_m, lat_s) = rad_to_deg(math.fabs(psi))
(lon_d, lon_m, lon_s) = rad_to_deg(math.fabs(lam))
return (lat_d, lat_m, lat_s, NS, lon_d, lon_m, lon_s, EW, h)
#Converts degrees to radians
def deg_to_rad(d, m, s):
tmp = d + m/60 + s/3600
return tmp*math.pi/180
#Converts radians to degrees
def rad_to_deg(r):
tmp = r*180/math.pi
d = math.floor(tmp)
m = math.floor((tmp - d)*60)
s = round(((tmp - d)*60 - m)*60, 4)
#NOTE: Edge case: rounding may
ing the values up to their max
if s == 60.0:
s = 0.0
m += 1
if m == 60:
m = 0
d += 1
return (d, m, s)
#Calculates the vehicle time given veh and sat postion and sat time
def calc_veh_time(xs, xv, sat_time):
veh_sat_dis = numpy.linalg.norm(xs-xv)
return veh_sat_dis/c + sat_time
#Calulates the position of the vehicle given the sat data
def calc_veh_pos(Sat_data):
#initialize the step norm to -1
s_norm = -1
#initialize xv to geo location of b12 lamp post at last ts
xv = geo_to_cart(prev_time, 40, 45, 55, 1, 111, 50, 58.0, -1, 1372.00)
while s_norm == -1 or s_norm > 0.01:
J = build_jacob(Sat_data, xv)
F = build_resid(Sat_data, xv)
#Calculate the step using numpy's least squares function
#J(xv)*s = -F(xv)
s = numpy.linalg.lstsq(J, -F)[0]
#xv = xv + s
xv = numpy.add(xv, s)
s_norm = numpy.linalg.norm(s)
tv = calc_veh_time(Sat_data[0].position, xv, Sat_data[0].time)
(lat_d, lat_m, lat_s, NS, lon_d, lon_m, lon_s, EW, h) = cart_to_geo(tv, xv)
lat_s = round(lat_s, 2)
lon_s = round(lon_s, 2)
h = round(h, 2)
return str(tv) + ' ' + str(lat_d) + ' ' + str(lat_m) + ' ' + str(lat_s) + ' ' + str(NS) + ' ' + str(lon_d) + ' ' + str(lon_m) + ' ' + str(lon_s) + ' ' + str(EW) + ' ' + str(h)
#Builds the Jacobian matrix
def build_jacob(Sat_data, xv):
jacobian = []
for sat in range(0, len(Sat_data) - 1):
Xs1 = numpy.a
ay(Sat_data[sat].position)
Xs2 = numpy.a
ay(Sat_data[sat+1].position)
Xs1_n = numpy.linalg.norm(Xs1 - xv)
Xs2_n = numpy.linalg.norm(Xs2 - xv)
J1 = (Xs1[0] - xv[0])/Xs1_n
J1 = J1 - ((Xs2[0] - xv[0])/Xs2_n)
J2 = (Xs1[1] - xv[1])/Xs1_n
J2 = J2 - ((Xs2[1] - xv[1])/Xs2_n)
J3 = (Xs1[2] - xv[2])/Xs1_n
J3 = J3 - ((Xs2[2] - xv[2])/Xs2_n)
jacobian.append( [J1,J2,J3] )
return numpy.a
ay(jacobian)
#Builds the residual vecto
def build_resid(Sat_data, xv):
residual = []
for sat in range(0, len(Sat_data) - 1):
Xs1 = numpy.a
ay(Sat_data[sat].position)
Xs2 = numpy.a
ay(Sat_data[sat+1].position)
Xv = numpy.a
ay(xv)
ts1 = Sat_data[sat].time
ts2 = Sat_data[sat+1].time
D1 = numpy.linalg.norm(Xs2 - Xv)
D2 = numpy.linalg.norm(Xs1 - Xv)
row = D1 - D2 - c*(ts1 - ts2)
residual.append(row)
return numpy.a
ay(residual)
####################################
#Parse the data file and set the variables
data_file = open(data_file, 'r')
for line in data_file:
line_splt = line.replace(',', '').split()
component = line_splt[2]
value = float(line_splt[0])
if component == 'pi':
pi = float(value)
elif component == 'c':
c = float(value)
elif component == 'R':
R = float(value)
elif component == 's':
s = float(value)
#Parse the piped input from satellite
prev_time = 0
Sat_data = []
for line in sys.stdin:
line_split = line.split()
ts = float(line_split[1])
#NOTE: A time diff of 1.0 seconds was given in the instructions,
# but .999 is used in np.dat, so we're setting .95 to be safe.
if prev_time is not 0 and abs(prev_time - ts) > .95:
#A new set of satellite transmissions has begun, process the cu
ent set before reading more input
print( calc_veh_pos(Sat_data) )
Sat_data = []
prev_time = ts
sat_num = int(line_split[0])
xs = float(line_split[2])
ys = float(line_split[3])
zs = float(line_split[4])
Xs = [xs,ys,zs]
new_sat_trans = Sat_trans(sat_num)
new_sat_trans.position = Xs
new_sat_trans.time = ts
new_sat_trans.sat_num = sat_num
Sat_data.append(new_sat_trans)
#Last set of transmissions
print( calc_veh_pos(Sat_data) )
66_python_gps_91485/satellite
#!/us
in/python3.4
#Satellite Program
import sys
import math
data_file = './data.dat'
class Satellite:
def __init__(self, name):
self.name = name
def position(self, t):
h = self.altitude
p = self.periodicity
phase = self.phase
u1 = self.u1
u2 = self.u2
u3 = self.u3
v1 = self.v1
v2 = self.v2
v3 = self.v3
x = (R+h)*(u1*math.cos((2*pi*t/p)+phase) + v1*math.sin((2*pi*t/p)+phase))
y = (R+h)*(u2*math.cos((2*pi*t/p)+phase) + v2*math.sin((2*pi*t/p)+phase))
z = (R+h)*(u3*math.cos((2*pi*t/p)+phase) + v3*math.sin((2*pi*t/p)+phase))
return (x,y,z)
def is_above_horizon(self, xv,yv,zv, t):
(xs,ys,zs) = self.position(t)
xs_1 = xs-xv
ys_1 = ys-yv
zs_1 = zs-zv
return (xv*xs_1)+(yv*ys_1)+(zv*zs_1) > 0
##########Helper Functions##########
def geo_to_cart(tv, lat_d, lat_m, lat_s, NS, lon_d, lon_m, lon_s, EW, h):
psi = NS*deg_to_rad(lat_d, lat_m, lat_s)
lam = EW*deg_to_rad(lon_d, lon_m, lon_s) + 2*pi*tv/s
r = (R+h)
x = r * math.cos(psi) * math.cos(lam)
y = r * math.cos(psi) * math.sin(lam)
z = r * math.sin(psi)
return (x,y,z)
#Converts degrees to radians
def deg_to_rad(d, m, s):
tmp = d + m/60 + s/3600
return tmp*math.pi/180
####################################
satellites = {}
data_file = open(data_file, 'r')
for line in data_file:
line_splt = line.replace(',', '').split()
component = line_splt[2]
value = float(line_splt[0])
if component == 'pi':
pi = float(value)
elif component == 'c':
c = float(value)
elif component == 'R':
R = float(value)
elif component == 's':
s = float(value)
else:
sat = int(line_splt[5])
if sat not in satellites:
new_sat = Satellite(sat)
satellites[sat] = new_sat
if component == 'v1':
satellites[sat].v1 = value;
elif component == 'v2':
satellites[sat].v2 = value;
elif component == 'v3':
satellites[sat].v3 = value;
elif component == 'u1':
satellites[sat].u1 = value;
elif component == 'u2':
satellites[sat].u2 = value;
elif component == 'u3':
satellites[sat].u3 = value;
elif component == 'periodicity':
satellites[sat].periodicity = value
elif component == 'altitude':
satellites[sat].altitude = value
elif component == 'phase':
satellites[sat].phase = value
#Parse the piped input from vehicle
for line in sys.stdin:
line_split = line.split()
#Vehical location at time tv
tv = float(line_split[0])
lat_d = float(line_split[1])
lat_m = float(line_split[2])
lat_s = float(line_split[3])
NS = int(line_split[4])
lon_d = float(line_split[5])
lon_m = float(line_split[6])
lon_s = float(line_split[7])
EW = int(line_split[8])
h = float(line_split[9])
#Compute the position of the vehicle in cartesian coordinates at time tv
(xv,yv,zv) = geo_to_cart(tv, lat_d, lat_m, lat_s, NS, lon_d, lon_m, lon_s, EW, h)
#Compute the satellites above the horizon at time tv
for sat_index in satellites:
sat = satellites[sat_index]
if not sat.is_above_horizon(xv,yv,zv, tv):
continue
#Compute ts and the satellites position at ts
tk = tv
tk_prev = -1
while tk - tk_prev:
(xs,ys,zs) = sat.position(tk)
distance = ( (xs-xv)**2 + (ys-yv)**2 + (zs-zv)**2 )**(.5)
tk_prev = tk
tk = tv - distance/c
ret = str(sat_index) + ' ' + str(tk) + ' ' + str(xs) + ' ' + str(ys) + ' ' + str(zs)
print(ret)
66_python_gps_91485/angles.class
public synchronized class angles {
double radians;
double seconds;
int degrees;
int minutes;
boolean plus;
static final double pi;
static final double twopi;
public void angles(double);
public void angles(int, int, int, boolean);
static double rad(int, int, double, int);
static double rad(int, int, double, boolean);
void report(String);
public static void main(String[]);
void check();
static void
();
}
66_python_gps_91485
12.dat
12123.0 12123.0 0 40 45 55 1 111 50 58 -1 1372
66_python_gps_91485
12t0.dat
0.0 0.0 0 40 45 55 1 111 50 58 -1 1372
66_python_gps_91485
m.dat
102123.0 112000.0 90 40 48 44 1 111 48 40 -1 2377
66_python_gps_91485/data.dat
3.141592653589793116E+00 /= pi
2.997924580000000000E+08 /= c, speed of light, [m/s]
6.367444500000000000E+06 /= R, radius of earth, [m]
8.616408999999999651E+04 /= s, length of a sidereal day, [s]
1.000000000000000000E+00 /= u1 of Sat. 0
0.000000000000000000E+00 /= u2 of Sat. 0
0.000000000000000000E+00 /= u3 of Sat. 0
0.000000000000000000E+00 /= v1 of Sat. 0
5.735764363510461594E-01 /= v2 of Sat. 0
8.191520442889917986E-01 /= v3 of Sat. 0
4.308204499999999825E+04 /= periodicity of Sat. 0 [s]
2.020000000000000000E+07 /= altitude of Sat. 0 [m]
0.000000000000000000E+00 /= phase of Sat. 0 [rad]
1.000000000000000000E+00 /= u1 of Sat. 1
0.000000000000000000E+00 /= u2 of Sat. 1
0.000000000000000000E+00 /= u3 of Sat. 1
0.000000000000000000E+00 /= v1 of Sat. 1
5.735764363510461594E-01 /= v2 of Sat. 1
8.191520442889917986E-01 /= v3 of Sat. 1
4.308204499999999825E+04 /= periodicity of Sat. 1 [s]
2.020000000000000000E+07 /= altitude of Sat. 1 [m]
1.570796326794896558E+00 /= phase of Sat. 1 [rad]
1.000000000000000000E+00 /= u1 of Sat. 2
0.000000000000000000E+00 /= u2 of Sat. 2
0.000000000000000000E+00 /= u3 of Sat. 2
0.000000000000000000E+00 /= v1 of Sat. 2
5.735764363510461594E-01 /= v2 of Sat. 2
8.191520442889917986E-01 /= v3 of Sat. 2
4.308204499999999825E+04 /= periodicity of Sat. 2 [s]
2.020000000000000000E+07 /= altitude of Sat. 2 [m]
3.141592653589793116E+00 /= phase of Sat. 2 [rad]
1.000000000000000000E+00 /= u1 of Sat. 3
0.000000000000000000E+00 /= u2 of Sat. 3
0.000000000000000000E+00 /= u3 of Sat. 3
0.000000000000000000E+00 /= v1 of Sat. 3
5.735764363510461594E-01 /= v2 of Sat. 3
8.191520442889917986E-01 /= v3 of Sat. 3
4.308204499999999825E+04 /= periodicity of Sat. 3 [s]
2.020000000000000000E+07 /= altitude of Sat. 3 [m]
4.712388980384689674E+00 /= phase of Sat. 3 [rad]
5.000000000000000000E-01 /= u1 of Sat. 4
8.660254037844385966E-01 /= u2 of Sat. 4
0.000000000000000000E+00 /= u3 of Sat. 4
-4.967317648921540929E-01 /= v1 of Sat. 4
2.867882181755230797E-01 /= v2 of Sat. 4
8.191520442889917986E-01 /= v3 of Sat. 4
4.308204499999999825E+04 /= periodicity of Sat. 4 [s]
2.020000000000000000E+07 /= altitude of Sat. 4 [m]
1.000000000000000000E+00 /= phase of Sat. 4 [rad]
5.000000000000000000E-01 /= u1 of Sat. 5
8.660254037844385966E-01 /= u2 of Sat. 5
0.000000000000000000E+00 /= u3 of Sat. 5
-4.967317648921540929E-01 /= v1 of Sat. 5
2.867882181755230797E-01 /= v2 of Sat. 5
8.191520442889917986E-01 /= v3 of Sat. 5
4.308204499999999825E+04 /= periodicity of Sat. 5 [s]
2.020000000000000000E+07 /= altitude of Sat. 5 [m]
2.570796326794896558E+00 /= phase of Sat. 5 [rad]
5.000000000000000000E-01 /= u1 of Sat. 6
8.660254037844385966E-01 /= u2 of Sat. 6
0.000000000000000000E+00 /= u3 of Sat. 6
-4.967317648921540929E-01 /= v1 of Sat. 6
2.867882181755230797E-01 /= v2 of Sat. 6
8.191520442889917986E-01 /= v3 of Sat. 6
4.308204499999999825E+04 /= periodicity of Sat. 6 [s]
2.020000000000000000E+07 /= altitude of Sat. 6 [m]
4.141592653589793116E+00 /= phase of Sat. 6 [rad]
5.000000000000000000E-01 /= u1 of Sat. 7
8.660254037844385966E-01 /= u2 of Sat. 7
0.000000000000000000E+00 /= u3 of Sat. 7
-4.967317648921540929E-01 /= v1 of Sat. 7
2.867882181755230797E-01 /= v2 of Sat. 7
8.191520442889917986E-01 /= v3 of Sat. 7
4.308204499999999825E+04 /= periodicity of Sat. 7 [s]
2.020000000000000000E+07 /= altitude of Sat. 7 [m]
5.712388980384689674E+00 /= phase of Sat. 7 [rad]
-4.999999999999998890E-01 /= u1 of Sat. 8
8.660254037844385966E-01 /= u2 of Sat. 8
0.000000000000000000E+00 /= u3 of Sat. 8
-4.967317648921540929E-01 /= v1 of Sat. 8
-2.867882181755230242E-01 /= v2 of Sat. 8
8.191520442889917986E-01 /= v3 of Sat. 8
4.308204499999999825E+04 /= periodicity of Sat. 8 [s]
2.020000000000000000E+07 /= altitude of Sat. 8 [m]
2.000000000000000000E+00 /= phase of Sat. 8 [rad]
-4.999999999999998890E-01 /= u1 of Sat. 9
8.660254037844385966E-01 /= u2 of Sat. 9
0.000000000000000000E+00 /= u3 of Sat. 9
-4.967317648921540929E-01 /= v1 of Sat. 9
-2.867882181755230242E-01 /= v2 of Sat. 9
8.191520442889917986E-01 /= v3 of Sat. 9
4.308204499999999825E+04 /= periodicity of Sat. 9 [s]
2.020000000000000000E+07 /= altitude of Sat. 9 [m]
3.570796326794896558E+00 /= phase of Sat. 9 [rad]
-4.999999999999998890E-01 /= u1 of Sat. 10
8.660254037844385966E-01 /= u2 of Sat. 10
0.000000000000000000E+00 /= u3 of Sat. 10
-4.967317648921540929E-01 /= v1 of Sat. 10
-2.867882181755230242E-01 /= v2 of Sat. 10
8.191520442889917986E-01 /= v3 of Sat. 10
4.308204499999999825E+04 /= periodicity of Sat. 10 [s]
2.020000000000000000E+07 /= altitude of Sat. 10 [m]
5.141592653589793116E+00 /= phase of Sat. 10 [rad]
-4.999999999999998890E-01 /= u1 of Sat. 11
8.660254037844385966E-01 /= u2 of Sat. 11
0.000000000000000000E+00 /= u3 of Sat. 11
-4.967317648921540929E-01 /= v1 of Sat. 11
-2.867882181755230242E-01 /= v2 of Sat. 11
8.191520442889917986E-01 /= v3 of Sat. 11
4.308204499999999825E+04 /= periodicity of Sat. 11 [s]
2.020000000000000000E+07 /= altitude of Sat. 11 [m]
6.712388980384689674E+00 /= phase of Sat. 11 [rad]
-9.999999999999997780E-01 /= u1 of Sat. 12
1.110223024625156540E-16 /= u2 of Sat. 12
0.000000000000000000E+00 /= u3 of Sat. 12
-5.551115123125782702E-17 /= v1 of Sat. 12
-5.735764363510460484E-01 /= v2 of Sat. 12
8.191520442889917986E-01 /= v3 of Sat. 12
4.308204499999999825E+04 /= periodicity of Sat. 12 [s]
2.020000000000000000E+07 /= altitude of Sat. 12 [m]
3.000000000000000000E+00 /= phase of Sat. 12 [rad]
-9.999999999999997780E-01 /= u1 of Sat. 13
1.110223024625156540E-16 /= u2 of Sat. 13
0.000000000000000000E+00 /= u3 of Sat. 13
-5.551115123125782702E-17 /= v1 of Sat. 13
-5.735764363510460484E-01 /= v2 of Sat. 13
8.191520442889917986E-01 /= v3 of Sat. 13
4.308204499999999825E+04 /= periodicity of Sat. 13 [s]
2.020000000000000000E+07 /= altitude of Sat. 13 [m]
4.570796326794896558E+00 /= phase of Sat. 13 [rad]
-9.999999999999997780E-01 /= u1 of Sat. 14
1.110223024625156540E-16 /= u2 of Sat. 14
0.000000000000000000E+00 /= u3 of Sat. 14
-5.551115123125782702E-17 /= v1 of Sat. 14
-5.735764363510460484E-01 /= v2 of Sat. 14
8.191520442889917986E-01 /= v3 of Sat. 14
4.308204499999999825E+04 /= periodicity of Sat. 14 [s]
2.020000000000000000E+07 /= altitude of Sat. 14 [m]
6.141592653589793116E+00 /= phase of Sat. 14 [rad]
-9.999999999999997780E-01 /= u1 of Sat. 15
1.110223024625156540E-16 /= u2 of Sat. 15
0.000000000000000000E+00 /= u3 of Sat. 15
-5.551115123125782702E-17 /= v1 of Sat. 15
-5.735764363510460484E-01 /= v2 of Sat. 15
8.191520442889917986E-01 /= v3 of Sat. 15
4.308204499999999825E+04 /= periodicity of Sat. 15 [s]
2.020000000000000000E+07 /= altitude of Sat. 15 [m]
7.712388980384689674E+00 /= phase of Sat. 15 [rad]
-5.000000000000000000E-01 /= u1 of Sat. 16
-8.660254037844383745E-01 /= u2 of Sat. 16
0.000000000000000000E+00 /= u3 of Sat. 16
4.967317648921539819E-01 /= v1 of Sat. 16
-2.867882181755230797E-01 /= v2 of Sat. 16
8.191520442889917986E-01 /= v3 of Sat. 16
4.308204499999999825E+04 /= periodicity of Sat. 16 [s]
2.020000000000000000E+07 /= altitude of Sat. 16 [m]
4.000000000000000000E+00 /= phase of Sat. 16 [rad]
-5.000000000000000000E-01 /= u1 of Sat. 17
-8.660254037844383745E-01 /= u2 of Sat. 17
0.000000000000000000E+00 /= u3 of Sat. 17
4.967317648921539819E-01 /= v1 of Sat. 17
-2.867882181755230797E-01 /= v2 of Sat. 17
8.191520442889917986E-01 /= v3 of Sat. 17
4.308204499999999825E+04 /= periodicity of Sat. 17 [s]
2.020000000000000000E+07 /= altitude of Sat. 17 [m]
5.570796326794896558E+00 /= phase of Sat. 17 [rad]
-5.000000000000000000E-01 /= u1 of Sat. 18
-8.660254037844383745E-01 /= u2 of Sat. 18
0.000000000000000000E+00 /= u3 of Sat. 18
4.967317648921539819E-01 /= v1 of Sat. 18
-2.867882181755230797E-01 /= v2 of Sat. 18
8.191520442889917986E-01 /= v3 of Sat. 18
4.308204499999999825E+04 /= periodicity of Sat. 18 [s]
2.020000000000000000E+07 /= altitude of Sat. 18 [m]
7.141592653589793116E+00 /= phase of Sat. 18 [rad]
-5.000000000000000000E-01 /= u1 of Sat. 19
-8.660254037844383745E-01 /= u2 of Sat. 19
0.000000000000000000E+00 /= u3 of Sat. 19
4.967317648921539819E-01 /= v1 of Sat. 19
-2.867882181755230797E-01 /= v2 of Sat. 19
8.191520442889917986E-01 /= v3 of Sat. 19
4.308204499999999825E+04 /= periodicity of Sat. 19 [s]
2.020000000000000000E+07 /= altitude of Sat. 19 [m]
8.712388980384689674E+00 /= phase of Sat. 19 [rad]
4.999999999999996669E-01 /= u1 of Sat. 20
-8.660254037844384856E-01 /= u2 of Sat. 20
0.000000000000000000E+00 /= u3 of Sat. 20
4.967317648921540374E-01 /= v1 of Sat. 20
2.867882181755229132E-01 /= v2 of Sat. 20
8.191520442889917986E-01 /= v3 of Sat. 20
4.308204499999999825E+04 /= periodicity of Sat. 20 [s]
2.020000000000000000E+07 /= altitude of Sat. 20 [m]
5.000000000000000000E+00 /= phase of Sat. 20 [rad]
4.999999999999996669E-01 /= u1 of Sat. 21
-8.660254037844384856E-01 /= u2 of Sat. 21
0.000000000000000000E+00 /= u3 of Sat. 21
4.967317648921540374E-01 /= v1 of Sat. 21
2.867882181755229132E-01 /= v2 of Sat. 21
8.191520442889917986E-01 /= v3 of Sat. 21
4.308204499999999825E+04 /= periodicity of Sat. 21 [s]
2.020000000000000000E+07 /= altitude of Sat. 21 [m]
6.570796326794896558E+00 /= phase of Sat. 21 [rad]
4.999999999999996669E-01 /= u1 of Sat. 22
-8.660254037844384856E-01 /= u2 of Sat. 22
0.000000000000000000E+00 /= u3 of Sat. 22
4.967317648921540374E-01 /= v1 of Sat. 22
2.867882181755229132E-01 /= v2 of Sat. 22
8.191520442889917986E-01 /= v3 of Sat. 22
4.308204499999999825E+04 /= periodicity of Sat. 22 [s]
2.020000000000000000E+07 /= altitude of Sat. 22 [m]
8.141592653589793116E+00 /= phase of Sat. 22 [rad]
4.999999999999996669E-01 /= u1 of Sat. 23
-8.660254037844384856E-01 /= u2 of Sat. 23
0.000000000000000000E+00 /= u3 of Sat. 23
4.967317648921540374E-01 /= v1 of Sat. 23
2.867882181755229132E-01 /= v2 of Sat. 23
8.191520442889917986E-01 /= v3 of Sat. 23
4.308204499999999825E+04 /= periodicity of Sat. 23 [s]
2.020000000000000000E+07 /= altitude of Sat. 23 [m]
9.712388980384689674E+00 /= phase of Sat. 23 [rad]
66_python_gps_91485/data.f
implicit double precision (a - h, o - z)
dimension u(3), v(3)
open (unit = 20, file = 'data.dat')
open (unit = 21, file = 'data.sht')
pi = 2.0d0*dacos(0.0d0)
write (20, 1000) pi
write (21, 1000) pi
c = 2.99792458d08
write (20, 1010) c
write (21, 1010) c
r = 6.3674445d6
write (20, 1020)
write (21, 1020)
s = 8.616409d04
write (20, 1030) s
write (21, 1030) s
v(1) = 0.0d0
v(3) = dsin(55.0d0/360.0d0*2.0d0*pi)
v(2) = dsqrt(1.0d0 - v(3)**2)
u(1) = 1.0d0
u(2) = 0.0d0
u(3) = 0.0d0
h = 2.0200d7
p = s/2.0d0
*
do 40 i = 1, 6
do 30 j = 1, 4
isat = (i - 1)*4 + j-1
do 10 k = 1, 3
write (20, 1040) u(k), k, isat
if (isat .eq. 0) write (21, 1040) u(k), k, isat
10 continue
do 20 k = 1, 3
write (20, 1050) v(k), k, isat
if (isat .eq. 0) write (21, 1050) v(k), k, isat
20 continue
write (20, 1070) p, isat
write (20, 1080) h, isat
phase = dfloat(i - 1) + dfloat(j - 1)*pi/2.0d0
write (20, 1090) phase, isat
if (isat .eq. 1) then
write (21, 1070) p, isat
write (21, 1080) h, isat
write (21, 1090) phase, isat
end if
30 continue
if (i .lt. 6) then
u1 = u(1)/2.0d0 - dsqrt(0.75d0)*u(2)
u2 = dsqrt(0.75d0)*u(1) + u(2)/2.0d0
u(1) = u1
u(2) = u2
v1 = v(1)/2.0d0 - dsqrt(0.75d0)*v(2)
v2 = dsqrt(0.75d0)*v(1) + v(2)/2.0d0
v(1) = v1
v(2) = v2
end if
40 continue
close (20)
close (21)
stop
1000 format (1pd26.18, 10x, ' /= pi ')
1010 format (1pd26.18, 10x, ' /= c, speed of light, [m/s] ')
1020 format (1pd26.18, 10x, ' /= R, radius of earth, [m] ')
1030 format (1pd26.18, 10x, ' /= s, length of a sidereal day, [s]')
1040 format (1pd26.18, 10x, ' /= u', i1, ' of Sat. ', i2)
1050 format (1pd26.18, 10x, ' /= v', i1, ' of Sat. ', i2)
1060 format (' check: plane ', i2, ' satellite ', i4, 2(' 1 = ', f4.2),
+ ' 0 /= ', d12.4)
1070 format (1pd26.18, 10x, ' /= periodicity of Sat. ', i2, ' [s]')
1080 format (1pd26.18, 10x, ' /= altitude of Sat. ', i2, ' [m]')
1090 format (1pd26.18, 10x, ' /= phase of Sat. ', i2, ' [rad]')
end
66_python_gps_91485/np.dat
152.3 152.3 15 90 0 0 1 150 0 0 1 0
66_python_gps_91485/o.dat
102123.2123 112000 9 0 0 0 1 0 0 0 1 0
66_python_gps_91485/README
This directory contains the following files (in addition to others):
README: this file
12.dat: trip to B12 in one step *
12t0.dat: trip to b12 in one step, at time 0 *
m.dat: trip to Black Mountain, in 90 steps *
data.dat: the data file
data.f: a program that writes the data file
np.dat: a trip to the north pole in 15 steps *
sp.dat: a trip to the South Pole in 15 steps *
o.dat: a trip to the point O in 9 steps *
eceiver: executable of my own receive
satellite: executable of my own satellite
seed.dat: used by my satellite for random number generation
vehicle: executable of my vehicle
Note: The executables run on our Sparc stations. You can use these
machines over the network. Let me know if you need an account. Use
the command
ssh xserver.math.utah.edu.
If for some reason you are unable to use them let me know. I may be
able to generate executables for your prefe
ed machine.
You can give commands like
cat bm.dat | vehicle | satellite | receiver (**)
This will create a file {\tt vehicle.dat} which should be identical to
the output generated by the receiver (except for the time
discrepancy). It should also generate files "satellite.log" and
"receiver.log". You may want to look at some of these.
The file 'bm.dat' in the command (**) can be replaced by any of the
above data files marked with an asterisk *. You can also create you
own data files. They have the following contents:
- beginning time
- ending time
- number of steps of trip
- degrees lattitude
- minutes lattidude
- seconds lattidude
- hemisphere (1 if Northern, -1 if Southern)
- degrees longitude
- minutes longitude
- seconds longitude
- hemisphere (1 if Eastern, -1 if Western)
This will generate a trip to the given position in the given number of
steps, starting at the lamp post B12.
Any of the three programs, vehicle, satellite, receiver, could be
eplaced by your own version.
To use these files first copy them into your own directory.
66_python_gps_91485
eceiver.class
public synchronized class receiver {
public void receiver();
public static void main(String[]);
static double cos(double);
static double sin(double);
static double sqrt(double);
static double sqrt(int);
static double atan(double);
static double abs(double);
static void write(String);
static void debug(String);
static double round(double);
}
66_python_gps_91485/satellite.class
public synchronized class satellite {
public void satellite();
public static void main(String[]);
static double cos(double);
static double sin(double);
static double sqrt(double);
static double abs(double);
static void write(String);
static void debug(String);
}
66_python_gps_91485/seed.dat
34
66_python_gps_91485/sp.dat
152.3 152.3 150 90 0 0 -1 150 0 0 1 0
66_python_gps_91485/sp.dat~
152.3 152.3 15 90 0 0 -1 150 0 0 1 0
66_python_gps_91485/tp.pdf
The future of GPS is really fascinating.
Hofmann-Wellenhof et al (2nd ed.), p. 288
Mathematics 5600—Summer 2015 −1−
List of Contents
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Sources of Information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
The Space Segment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
Coordinate Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
The Programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
vehicle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
satellite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
eceiver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Mathematical Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
What To Do . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Procrastination is Hazardous . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Simplifying Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
Appendix: The Data File, data.dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
Introduction
The primary numerical task in this semester project is the solution of certain systems
of nonlinear equations. However, beyond that the project illustrates a recent technological
development which I hope you will find interesting and motivating. To get everything to
work in this project we will have to put many pieces together, and they have to work just
ight, as in most real life problems. The results of your work, a couple of programs, will
need to work with the programs of your class mates and my own. When we reach that stage
we will be able to say with assurance that we truly understand the problem and its solution.
The original US Global Positioning System−2− (GPS) consists of 24 satellites that o
it
the earth and constantly transmit their position and the precise time of that transmission,
in addition to other data. Commercially available receivers, (Global Positioning Devices)
process that information and, using the differences in the runtimes of the satellite signals,
compute their position and the precise time. The position can be expressed, e.g., as longi-
tude, latitude, and altitude. Its accuracy depends on the type and quality of the device and
its software, and other circumstances. It ranges from less than one centimeter to about 100
meters−3−.
In this project, we will write two programs that simulate a much simplified version of
this system.
GPS was developed originally for the US armed forces, but there are of course also
civilian and scientific applications, many of which are still being discovered. As far as volume
and resources involved, civilian applications now vastly outstrip military applications. The
following is an incomplete list of applications I have read or heard about (or, in the case of
the first three applications, actually practiced):
−1− by Peter Alfeld, [email protected], 801-581-6842, JWB 127, TEX processing date:
May 20, 2015. Please let me know of any e
ors or inconsistencies you find in this
document.
−2− Full Disclosure: GPS is a multibillion dollar industry. I have no commercial connections
or interests in GPS. This project is strictly for scientific purposes, and specific GPS devices
or publications are mentioned or shown only for illustration and information. I’m not
endorsing or recommending any specific GPS devices.
−3− A meter is a little more than 3 feet. Most of the literature on GPS uses the metric
system, and we follow that convention.
Mathematics 5600 Summer 2015 Semester Project page 1
• Navigating in a car. This may be the largest civilian application at present. You can
uy devices that will tell you where you are and what’s near you, and give you oral
turn by turn instructions on how to get to your destination. Many high end models
now come with built-in GPS based navigation systems.
• Finding your way in the wilderness.
• Marine navigation.
• Geocaching is a high-tech treasure hunting game played throughout the world by ad-
venture seekers equipped with GPS devices. This is a highly popular past time where
people search for a cache established by others, given it’s coordinates, remove a token
item, and replace it with another. The basic idea is to locate hidden containers, called
geocaches, outdoors and then share your experiences online. The official web site is
http:
www.geocaching.com/ .
• Electronically marking a spot (a good fishing place, a buried treasure or mine, man o
item ove
oard) in a featureless area such as the ocean, a lake, a desert, or the Salt
Flats.
• Aerial navigation, including precision approach and landing.
• Measuring continental drift and expansion.
• Automatic grading and paving in road construction.
• Telling blind people (via Braille or Audio) where they are.
• Steering a tractor planting vine plants along a suitable trajectory.
• Surveying with high accuracy.
• Controlling a fleet of vehicles (e.g., Police, Taxi, Truck or Bus companies).
• Transfer of extremely accurate time information.
• Installation of several GPS devices in a vehicle, like a ship or plane, and determination
not just of position but also of attitude, roll, pitch, and yaw.
• By coupling a cell phone with a GPS device, in an emergency it can tell your location
to the dispatcher even if you don’t know it yourself. If the phone is installed in you
car, and the car is stolen, you might be able to call your car and ask where it is.
GPS is available anywhere on or near earth, 24 hours a day, in fog, space, or darkness.
By making several measurements over time one can compute and display such quantities as
speed, bearing, or estimated time of a
ival. It’s eerie to have your car give you detailed
instructions on how to get where you are going, as you are going. Or, for example, you can
transform a simple rental boat into a sophisticated yacht with compass and speedomete
simply by placing a GPS device next to the wheel.
Sources of Information
• You can download my software to test yours. For details see
http:
www.math.utah.edu/~pa/5600/tp
If you are interested in GPS and like to learn more, here are some places to get started:
• A standard text on GPS is B. Hofmann-Wellenhof, H. Lichtenegger, and J. Collins:
GPS, Theory and Practice, Springer-Verlag Wien New York, 5th edition, 2001, ISBN
3211835342. This excellent but very technical book gives lots of details and references.
• Another standard detailed and technical text is Elliott D. Kaplan (ed.), Understand-
ing GPS Principles and Applications, Artech House, Boston, London, 1996, ISBN
0890067937.
Mathematics 5600 Summer 2015 Semester Project page 2
• There’s a great deal of information on the web. The following Table shows the numbe
of results found by a Google search for GPS.
May 2004 17, 000, 000
May 2005 58, 600, 000
August 2007 257, 000, 000
May 2008 407, 000, 000
August 2009 244, 000, 000
May 2012 1, 330, 000, 000
May 2013 720, 000, 000
May 2014 298, 000, 000
May 2015 620, 000, 000
Table: Number of Google results for GPS
• There is excellent information in a concise and readable form about the shape of the
earth and its gravitational field in the Encyclopedia Britannica, particularly in the
article The Earth in the macropedia.
• A useful text giving all kinds of mathematical formulas and information is the VNR
Concise Encyclopedia of Mathematics. It’s published by Van Nostrand Reinhold Com-
pany. My copy was published in 1975. The ISBN is 0442226462.
Note on notation
Notation can get quite cumbersome in this project. I use boldface letters to denote
vectors in cartesian coordinates. The coordinates of these vectors themselves are lower case
Roman letters that may be subscripted. For example, the vector x is usually denoted by
x
y
z
(1)
ut sometimes, for example in the discussion of Newton’s method, by
x1
x2
x3
. (2)
Subscripts are often used with an S to denote co
espondence to a satellite. for example,
xSi denotes co
espondence to the i-th of several satellites.
The Space Segment
The number and configuration of the o
iting satellites is evolving. We will consider a
set of 24 satellites o
iting in six planes at an inclination of 55◦ to the equator at an altitude
of about 20,200 km.
To help you get going on the project, sprinkled throughout this assignment are a total
of 17 exercises. You will have to solve these exercises or closely related problems during
the work on the project, so they collectively form our first home work in this class.
The Model
The model will consist of three modules, i.e., programs, the vehicle, the satellite,
and the receiver. The vehicle generates positional data, the satellite converts those
into data of the kind processed by GPS receivers, and the receiver converts these back
into positional data. The three modules are piped together in the standard Unix fashion:
vehicle | satellite | receiver (3)
Mathematics 5600 Summer 2015 Semester Project page 3
If all goes well the standard output of receiver will equal (almost) the standard output
of vehicle. Your job is to write satellite and receiver (and perhaps a rudimentary
vehicle for testing purposes).
The model depends on a set of data provided in a file called
data.dat
which is read by satellite. The beginning part of that file is also read by receiver.
Figure 1 summarizes the model. A
ows indicate data flow.
data.dat
↓ ց
vehicle −→ satellite −→ receive
Figure 1: The Model
Following are some of the assumptions we’ll make to construct our model. Numerical
values are given here for illustration. They will be actually read from data.dat and the
precise numerical values are subject to change (with appropriate notification). You or I may
also change those parameters for the purpose of testing our programs.
1. The Earth is perfectly spherical and has a radius R of
R = 6, 367, 444.50 m. (4)
2. The Earth turns at a constant rate and completes one revolution in one sidereal day−4−
of
s = 86, 164.09 seconds. (5)
3. The satellites move at a constant speed in perfectly circular o
its at an altitude of
h = 20, 200, 000 m (6)
and with an o
ital period
p =
s
2
(7)
of exactly half a sidereal day.
4. The satellites are evenly spaced in groups of four in six planes, each of which is inclined
55◦ to the equator. According to Hofmann-Wellenhof et al this provides global coverage
of four to eight satellites being more than 15◦ above the horizon at all times and in all
places on earth.
−4− A sidereal day is the time required for a complete rotation of the earth in reference to
a fixed star. It differs from a standard solar day of 24 hours because in the course of a day
the earth not only rotates, but also changes position in its o
it around the Sun. Hence the
sun changes its position in the sky and the earth needs to rotate an additional 4 minutes o
so to catch up with the sun before it is noon again.
Mathematics 5600 Summer 2015 Semester Project page 4
Coordinate Systems
We use two coordinate systems: longitude, latitude and altitude to express positions
on and near earth, and a cartesian coordinate system to describe the space in which the
satellites o
it and the earth rotates. Like commercial global positioning devices, we express
angles in degrees, minutes of degree, and seconds of degree. The two coordinate systems
are linked by the following conventions:
The North Pole is located at (0, 0, R), the South Pole at (0, 0,−R), (assuming both
have altitude 0, which is definitely wrong for the South Pole). At time 0 the point O of zero
longitude, latitude, and altitude, is at (R,0,0), and the Earth rotates (west to east) once in
a sidereal day.
Exercise 1: Find a formula that describes the trajectory of the point O in cartesian coordinates
as a function of time.
The Programs
As stated above, our model will consist of three programs:
1. The Vehicle. This program, which you can download from our web page, produces a
stream of data sets (to standard output) each of which has the form:
tV ψd ψm ψs NS λd λm λs EW h (8)
where tV denotes Universal time
−5− in seconds, ψd, ψm, ψs is latitude in degrees,
minutes, and seconds, and, similarly, λd, λm, λs is longitude in degrees, minutes, and
seconds. More specifically,
tV is a real number given to an accuracy of 10
−2 seconds ranging from 0 to 106. This
is the time at which the vehicle is at the specified position.
ψd is an integer ranging from 0
◦ (i.e., the Equator) to +90◦ (i.e., the North or South
Pole).
ψm is an integer ranging from 0 to 59 minutes of degree.
ψs is a real number ranging from 0 to 59.9999 seconds of degree. It should be given
to an accuracy of 10−2 (which co
esponds to an accuracy of about a foot).
NS is an integer that is +1 North of the equator and −1 South of the equator.
λd is an integer ranging from 0
◦ (i.e., the meridian of Greenwich) to 180 (i.e., 180
degrees east, or west, the date line).
λm is an integer ranging from 0 to 59 minutes of degree.
λs is a real number ranging from 0 to 59.9999 seconds of degree, given to the same
accuracy as ψs.
EW is an integer that is +1 east of Greenwich and −1 west of Greenwich.
h is a real number giving the altitude in meters, to an accuracy of 1cm.
−5− The mean solar time of the Greenwich meridian (0◦ longitude). Universal time replaced
the designation Greenwich mean time in 1928; it is now used to denote the solar time when
an accuracy of about 1 second suffices. In 1955 the International Astronomical Union defined
several categories of Universal Time of successively increasing accuracy. UT0 represents the
initial values of Universal Time obtained by optical observations of star transits at various
astronomical observatories. These values differ slightly from each other because of the effects
of polar motion. UT1, which gives the precise angular coordinate of the Earth about its spin
axis, is obtained by co
ecting UT0 for the effects of polar motion. Finally, an empirical
co
ection to take account of annual changes in the Earth’s speed of rotation is added to
UT1 to convert it into UT2. Coordinated Universal Time, the international basis of civil
and scientific time, is obtained from an atomic clock that is adjusted so as to remain close
to UT1; in this way, the solar time that is indicated by Universal time is kept in close
coordination with atomic time. [Encyclopedia Britannica, 1992]. We will assume that all
the times described here are identical and call them Universal Time.
Mathematics 5600 Summer 2015 Semester Project page 5
For example, according to my Magellan Trailblazer the street light labeled B12−6− in
front of the South Window of my office (at time t) is located at:
t 40 45 55.0 1 111 50 58.0 − 1 1372.00 (9)
i.e., at latitude 40◦ 45’ 55” North, longitude 111◦ 50’ 58” West, and an altitude of 1372
m.
Exercise 2: Write a program that converts angles from degrees, minutes, and seconds to radi-
ans, and vice versa. Make sure your program does what it’s supposed to do.
For the following four exercises assume that tV equals true Universal time and denote
it by t.
Exercise 3: Find a formula that converts position as given in (8) at time t = 0 into cartesian
coordinates.
Exercise 4: Find a formula that converts position and general time t as given in (8) into
cartesian coordinates.
Exercise 5: Find a formula that converts a position given in cartesian coordinates at time t = 0
into a position of the form (8).
Exercise 6: Find a formula that converts general time t and a position given in cartesian
coordinates into a position of the form (8).
Exercise 7: Find a formula that describes the trajectory of lamp post B12 in cartesian coordi-
nates as a function of time.
The vehicle should not produce impossible positions (like a latitude of 90◦ 59’ 59”).
For testing purposes we will use vehicles that produce a stream of data co
esponding,
for example, to a walk from JWB to the Ma
iott Li
ary, a hike up to Mount Olympus,
or a flight from Salt Lake City to the North Pole. The data sets should be spaced at
least 1 second (of time) apart.
In addition vehicle writes a copy of the standard input and the standard output into
the file
vehicle.log.
2. The Satellite. This program reads data from data.dat and then reads the data
generated by a vehicle and processes those data as follows:
A. It computes the position xV of the vehicle in cartesian coordinates (measured in
meters) at the time tV .
B. For each satellite S that is above the horizon−7− at the vehicle’s position it com-
putes the time tS at which a signal needs to be sent to reach the vehicle at time
tV and position xV , assuming the signal moves at the speed c of light in vacuum
where
c = 2.99792458× 108m/s. (10)
The satellite needs to be above the horizon at time tS , The satellite also computes
the position xS of the satellite at the time tS It then writes the following data to
standard output:
iS tS xS (11)
where iS is the reference number of the satellite (ranging from 0 to 23). The time
tS should be given to an accuracy of 10
−11 seconds, and the position xS should be
given as a triple of cartesian coordinates with an accuracy of 1cm.
−6− If you look around on campus you’ll see that all street lights have been labeled with a
letter and a number. On some lights, the labels have disappeared, or, as in the case of B12,
obstructed by an attachment.
−7− My own receiver will check this condition. I figured that assuming a transparent earth
would be too much of a simplification. However, calculations in the actual GPS system are
usually based on satellites at least 15◦ above the horizon.
Mathematics 5600 Summer 2015 Semester Project page 6
C. In addition, satellite writes a copy of the standard input and the standard
output into the file
satellite.log.
It should terminate gracefully when the data stream from the satellite program is
terminated.
Exercise 8: Given a point x on earth and a point s in space, both in cartesian coordinates,
find a condition that tells you whether s as viewed from x is above the...