How to correctly make a fake, counter-propagating TLE?

Question

In this answer I have said:

Do not make FAKE TLEs!

I would now like to make a fake TLE.

Specifically I'd like to generate a counter-propagating "twin" of an orbit by flipping inclination by 180° but the problem is that the TLE specification only allows inclination between 0° and 180° to avoid ambiguities.

So I need to further malign the fake TLE so that both the real and fake TLEs will yield the same position at the TLE's epoch, but at other times will move in opposite directions.

What is the correct way to do this unspeakable deed?

This turns out to be trickier than I'd thought, and so I'm at a standstill.

---

I'm looking at the question How presise could 2 satellites “thread the needle” in orbit? and trying to do some simple tests with SGP4 for fun.

Answer 1

Using the diagram from Wikpedia

When the orbit reverses in direction the ascending node becomes the descending node, and vice versa. Fortunately the ascending node is always 180 degrees around from the descending, so that's easy to calculate, add 180 degrees, and subtract 360 if the result is greater than 360. so

$an = (an + 180) \% 360$
$_{_{an += an < 180 ? 180 : -180}} $

The mean anomaly is a fudge which basically describes what fraction of the orbital period has elapsed since the object passed periapsis. when the direction reverses, this becomes the time until it reaches periapsis, in order to turn this into a time after passing periapsis we need to add a whole revolution on, so ma = -ma +360 ... or

$ma = 360 - ma$

The argument of periapsis is the angle between the ascending node and the position of periapsis, in the direction of travel. because we've reversed the direction of travel, we need to negate it ap = -ap. As mentioned already, the ascending node has shifted 180 degrees, so we need to add that on ap = 180 - ap. Because computers can be a bit funny if values unexpectedly go negative, we'll add another 360 degrees on just in case, then take it off again if it wasn't needed

$ap = (540 - ap) \% 360$
$_{_{ap = (ap < 180 ? 180 : 540) -ap}}$

The inclination is annoying. If it were the angle the satelite crosses the equator at the ascending node (per the diagram), it would remain the same. If it were the the highest and lowest latitude it reached, it would remain the same. But no. it is 0-90 for a prograde satelite and 90-180 for a retrograde satelite.

$i = 180 - i$

Having guessed these answers, I played around with http://orbitsimulator.com/formulas/OrbitalElements.html and was surprised to discover it seemed to agree.

Answer 2

This is a supplementary answer to confirm that the procedure suggested in @JCRM's answer seems to work for me as well.

I do not recommend you use the following script; it's here only to illustrate how the test was done. Oh, and never try to make a fake TLE!

TLE_ISS = """ISS (ZARYA)             
1 25544U 98067A   18343.93002315  .00002202  00000-0  40582-4 0  9991
2 25544  51.6408 222.9207 0005154 129.7472  46.7852 15.54076162145848"""

TLE_GPS = """GPS BIIR-5  (PRN 28)    
1 26407U 00040A   18343.33921716 -.00000070  00000-0  00000+0 0  9991
2 26407  56.3821 324.0376 0198177 274.5104 260.5068  2.00561560134849"""

TLE_TDRS = """TDRS 13                 
1 42915U 17047A   18342.83940260 -.00000286  00000-0  00000+0 0  9996
2 42915   6.3999 330.9155 0027242   6.2708 353.7662  1.00270146  4833"""

TLEs = (TLE_ISS, TLE_GPS, TLE_TDRS)

def numit(x, nleft, nright):
    """this is awful, don't ever do this!"""
    neg        = x < 0.
    x          = abs(x)
    whole      = int(abs(x))
    fraction   = x - whole
    left  = str(whole + 10**(nleft))[1:]
    if neg:
            left = '-'+left[1:] 
    right = str(fraction)[2:2+nright]
    return left + '.' + right

class Sat(object):
    def __init__(self, lines):
        self.ok   = False
        self.okok = False
        lines = lines.splitlines()
        if len(lines) == 3:
            name, L1, L2 = lines
        elif len(lines) == 2:
            L1, L2 = lines
        else:
            print "not enough lines"
            exit

        if len(L1)==69 and len(L2)==69:
            self.L1 = L1
            self.L2 = L2
            self.ok = True
            if name and len(name)>0:
                self.name = name
        else:
            print "lines not correct length"

        if self.ok:
            try:
                epoch_year = 1900 + int(self.L1[18:20])
                if epoch_year < 1957:
                    epoch_year       += 100
                epoch_decimal_days    = float(self.L1[20:32])
                epoch_datetime        = (datetime(epoch_year, 1, 1) +
                                         timedelta(epoch_decimal_days - 1)) # https://stackoverflow.com/a/34910712/3904031
                self.epoch_datetime   =  epoch_datetime

                revs_per_day          = float(self.L2[52:63])
                self.period_days      = 1./revs_per_day
                self.period_hours     = 24./revs_per_day
                self.period_minutes   = 60.*24./revs_per_day
                self.period_seconds   = 3600.*24./revs_per_day

                raanstr               = self.L2[17:25]
                mastr                 = self.L2[43:51]
                arpestr               = self.L2[34:42]
                incstr                = self.L2[ 8:16]

                self.raan             = float(raanstr)
                self.ma               = float(mastr)
                self.arpe             = float(arpestr)
                self.inc              = float(incstr)

# ascending node:  raan = (raan + 180)%360
# mean anomaly:  ma = 360 - ma
# argument of periapsis:  ap = (540 - ap)%360
# inclination:   inc = 180 - i
# https://space.stackexchange.com/questions/32725/how-to-correctly-make-a-fake-counter-propagating-tle

                self.raan_flip        = (self.raan + 180)%360.
                self.ma_flip          = (360. - self.ma)
                self.rpe_flip         = (540. - self.arpe)%360.
                self.inc_flip         = (180. - self.inc)

                L1f                   = self.L1[:]  # copy

                L2f                   = self.L2[:]  # copy
                L2f = L2f[:17] + numit(self.raan_flip, 3, 4) + L2f[25:]
                L2f = L2f[:43] + numit(self.ma_flip,   3, 4) + L2f[51:]
                L2f = L2f[:34] + numit(self.rpe_flip,  3, 4) + L2f[42:]
                L2f = L2f[:8]  + numit(self.inc_flip,  3, 4) + L2f[16:]


                if len(L1f)==69 and len(L2f)==69:
                    self.L1_flip = L1f
                    self.L2_flip = L2f
                    self.okok = True
                else:
                    print "flipped lines not correct length"

            except:
                print "trouble in river city"


import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Loader, Topos, EarthSatellite
from datetime import datetime, timedelta

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)] 
degs, rads = 180/pi, pi/180

load = Loader('~/Documents/fishing/SkyData')  # avoids multiple copies of large files
ts   = load.timescale()

data    = load('de421.bsp')
Earth   = data['Earth']
ts      = load.timescale()

sats = []
for TLE in TLEs:
    sat = Sat(TLE)
    if sat.ok:
        sats.append(sat)

print "len(sats): ", len(sats)
print "len([s for s in sats if s.okok]): ", len([s for s in sats if s.okok])

for sat in sats:

    sat.Th        = sat.period_hours
    sat.dhours    = np.arange(-1.2*sat.Th, 1.2*sat.Th, 0.001*sat.Th)

    epoch         = sat.epoch_datetime

    year, month, day  = epoch.year, epoch.month, epoch.day
    hours             = epoch.hour + sat.dhours
    minute, second    = epoch.minute, epoch.second + 1E-06*epoch.microsecond
    sat.times         = ts.utc(year, month, day, hours, minute, second)

    sat.satobj  = EarthSatellite(sat.L1,      sat.L2     )
    sat.satobjf = EarthSatellite(sat.L1_flip, sat.L2_flip)
    sat.geocen  = sat.satobj.at(sat.times)
    sat.geocenf = sat.satobjf.at(sat.times)
    sat.posn    = sat.geocen.position.km
    sat.posnf   = sat.geocenf.position.km
    # sat.path    = sat.geocen.subpoint()
    # sat.pathf   = sat.geocenf.subpoint()

if True:
    plt.figure()

    for j, sat in enumerate(sats[:3]):

        plt.subplot(3, 1, j+1)

        colors = '-b', '-g', '-r'

        h  = max(np.abs(sat.posn).max(), np.abs(sat.posnf).max())

        for i, (a, b, c) in enumerate(zip(sat.posn, sat.posnf, colors)):
            # plt.subplot(3, 1, i+1)
            plt.plot(sat.dhours, a, c)
            plt.plot(sat.dhours, b, c)
            plt.plot([      0,       0], [-h, h], '-k')
            plt.plot([-sat.Th, -sat.Th], [-h, h], '-k')
            plt.plot([ sat.Th,  sat.Th], [-h, h], '-k')
        plt.title(sat.name)
    plt.show()