Compute collisions/encounters between satellites in a catalog

cola.py

#!/usr/bin/env python3
import numpy as np
from sgp4.earth_gravity import wgs84
from sgp4.io import twoline2rv
from astropy.time import Time

class TwoLineElement:
    """TLE class"""

    def __init__(self, tle0, tle1, tle2):
        """Define a TLE"""

        self.tle0 = tle0
        self.tle1 = tle1
        self.tle2 = tle2
        if tle0[:2] == "0 ":
            self.name = tle0[2:].rstrip()
        else:
            self.name = tle0.rstrip()
        if tle1.split(" ")[1] == "":
            self.id = int(tle1.split(" ")[2][:4])
        else:
            self.id = int(tle1.split(" ")[1][:5])

def read_tles(fname):
    with open(fname, "r") as fp:
        lines = fp.readlines()

        tles = []
        for i in range(0, len(lines), 3):
            tles.append(TwoLineElement(lines[i], lines[i+1], lines[i+2]))


    return tles

if __name__ == "__main__":
    # Read TLEs
    tles = read_tles("/data2/tle/catalog.tle")

    # Time
    t = Time("2019-11-16T10:00:00", format="isot", scale="utc")

    # Compute positions
    p = np.zeros(len(tles)*3).reshape(len(tles), 3).astype("float32")
    for i, tle in enumerate(tles):
        sat = twoline2rv(tle.tle1, tle.tle2, wgs84)
        
        p[i], v = sat.propagate(t.datetime.year, t.datetime.month, t.datetime.day,
                                t.datetime.hour, t.datetime.minute, t.datetime.second)
        

    # Compute ranges
    dp = (p[:, np.newaxis]-p[np.newaxis, :])
    r = np.sqrt(np.sum(dp**2, axis=2))

    idxs = np.where((r>0.0) & (r<20.0))
    idx1, idx2 = idxs[0], idxs[1]
    
    for i in range(len(idx1)):
        j, k = idx1[i], idx2[i]
        if j<k:
            print(r[j, k], tles[j].name, tles[k].name, tles[j].id, tles[k].id)