### Quintecwriteups and other thoughts

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 *

path = []

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

for d in dat:
time, az, el, rng = d.split()
t2 = time.split(':')
mins = int(t2)
secs = int(t2.split('.'))
stamp = ts.utc(2021, 6, 27, 0, mins, secs)
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)
ys.append(sat)
zs.append(sat)
``````

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,
2*np.pi / P,
)
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)
# print((path - path))

#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.