Quintec writeups and other thoughts

Mr. Radar

27 Jun 2021

We are given the location of a radar and several radar pulse returns of a satellite, and want to determine its orbital parameters. We first trace out the path of the satellite, using Skyfield to do the appropriate calculations. The functions latlon and from_altaz were particularly helpful to be able to convert the data given directly into ICRF coordinates.

from skyfield.api import *
from skyfield.units import *
from skyfield.positionlib import *

ts = load.timescale(builtin=True)

radar = wgs84.latlon(8.7256, 167.715)
radar.elevation = Distance(m=35)

with open('radar_data.txt', 'r') as f:
	dat = f.readlines()[1:]

path = []

xs = []
ys = []
zs = []

for d in dat:
	time, az, el, rng = d.split()
	t2 = time.split(':')
	mins = int(t2[1])
	secs = int(t2[2].split('.')[0])
	stamp = ts.utc(2021, 6, 27, 0, mins, secs)
	rp = radar.at(stamp)
	pos = rp.from_altaz(alt_degrees=float(el), az_degrees=float(az), distance=Distance(km=float(rng)))
	sat = pos.position.km + rp.position.km
	path.append(sat)
	xs.append(sat[0])
	ys.append(sat[1])
	zs.append(sat[2])

As a sanity check, we visualize the orbit with matplotlib.

from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlim(4500, 5500)
ax.set_ylim(6700, 7700)
ax.set_zlim(-2000, -3000)
ax.scatter3D(xs, ys, zs)
plt.show()

We then use a least squares regression to fit an elliptical orbit to our calculated coordinates, using Skyfield’s Satrec to quickly trace out an orbit given the orbital parameters. We use this calculator to estimate initial parameters for the regression, using the first two data points to generate a position + velocity vector to convert into orbital parameters.

def cost(params, debug=False):
	M = 5.9722e24
	G = 6.6743015e-11
	a, e, i, Omega, omega, M0 = params
	P = np.sqrt(4*np.pi**2 / (G*M) * (1000*a)**3) / 60

	satrec = Satrec()
	satrec.sgp4init(
		WGS72,
		'i',
		1337,
		t0.tdb - 2433281.5,
		0.0, # drag
		0.0,
		0.0,
		e,
		np.radians(omega),
		np.radians(i),
		np.radians(M0),
		2*np.pi / P,
		np.radians(Omega),
	)
	sat = EarthSatellite.from_satrec(satrec, ts)
	if debug:
		return sat

	err = []
	for i in range(100):
		t = ts.utc(2021, 6, 27, 0, 8, 12+i)
		err += list(path[i] - sat.at(t).position.km)

	return err

# print(path[0])
# print((path[1] - path[0]))

#initial parameters
a = 20804.737560848476
e = 0.6610543599443122
i = 33.99839311086964
Omega = 78.00087926523055
omega = 271.65368202323253
M0 = 10.704246168858505
sol = least_squares(cost, np.array([a, e, i, Omega, omega, M0]))
print(sol)

params = [2.29939092e+04, 7.00821959e-01, 3.38773730e+01, 7.82120294e+01, 2.70127740e+02, 9.77239867e+00] #from solution
print(cost(params))
sat = cost(params, debug=True)
t = ts.utc(2021, 6, 27, 0, 9, 52)
#sanity checks
print(sat.at(t).position.km)
print(sat.at(t).velocity.km_per_s)

The regression gives us the orbital parameters that best fit our data, and submitting those to the server gives us the flag.

Related posts