<jupyter_start><jupyter_text>Plot the orbit of S2 from the precession rate found in [this article](https://www.aanda.org/articles/aa/abs/2020/04/aa37813-20/aa37813-20.html)<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from astropy import units, constants
%matplotlib inline
#Table E.1
fSP = 1.1
fRS = 1.
MSGa = 4.261e6*units.solMass
R0 = 8246.7*units.pc
a = 125.058*units.mas
e = 0.884649
i = (134.567*units.degree).to(units.radian).value
omega = (66.263*units.degree).to(units.radian).value #argument of periapsisi
OMEGA = (228.171*units.degree).to(units.radian).value #longitude of ascending node
P = 16.0455*units.yr
tPeri = 2018.379*units.yr
x0 = -0.9*units.mas
y0 = 0.07*units.mas
vx0 = 0.08*units.mas/units.yr
vy0 = 0.0341*units.mas/units.yr
vz0 = -1.6*units.km/units.s

#calculate dphi using Eq. C.3
# RS = (2.*constants.G*MSGa/constants.c**2.).decompose() #Schwarszchild radius
# RSangle = np.arctan2(RS, R0).to(units.mas)
# dphi = 3.*fSP*np.pi*RSangle/(a*(1. - e**2.)) #I think this is the angular change in the location of periastron per orbit? but what are the units? radians per orbit?
# #is dphi meant to be applied to the argument of periapsis? https://en.wikipedia.org/wiki/Apsidal_precession
# print(RS)
# print((dphi*units.radian).to(units.arcmin))
#I get about 13.45 arcmin/orbit; and I just noticed that the paper quotes 12 arcmin/orbit in the abstract.  Well, it was a good excersize :)

dphi = (12*units.arcmin).to(units.radian).value

#how many orbits do we need to get all the way back around?
print("N orbits to repeat = ", 2.*np.pi/dphi)
def createOrbit(semi, ecc, inc, lan, ap, tperi, period, dphi, times):
    def getPQ(b):
        c = np.cos(b)
        s = np.sin(b)

        P = np.empty(3);
        P[0] = -1.*c[2]*c[1] + s[2]*c[0]*s[1]
        P[1] = -1.*c[2]*s[1] - s[2]*c[0]*c[1]
        P[2] = -1.*s[2]*s[0]
        Q = np.empty(3);
        Q[0] = s[2]*c[1] + c[2]*c[0]*s[1]
        Q[1] = s[2]*s[1] - c[2]*c[0]*c[1]
        Q[2] = -1.*s[0]*c[2]
        return P, Q
    
    b = np.array([-1.*inc, lan, ap])
    dct = {"x":[], "y":[], "z":[], "t":[]}
    for t in times:
        tdiff = t - tperi
        phase = (tdiff % period)/period
        
        b[2] = ap + dphi*(tdiff/period)
        P, Q = getPQ(b)
        
        E = ((2.*phase*np.pi) % (2.*np.pi)).value
        pos = semi*( (np.cos(E) - ecc)*P + np.sqrt(1.0 - ecc**2.)*np.sin(E)*Q )
        
        dct["x"].append(pos[0].value)
        dct["y"].append(pos[1].value)
        dct["z"].append(pos[2].value)
        dct["t"].append(t.value)

    
    return pd.DataFrame(dct)

Norbits = 100
Nvals = 100*Norbits
times = np.linspace(2020, 2020 + Norbits*P.value, Nvals)*units.yr

orbit = createOrbit(a, e, i, OMEGA, omega, tPeri, P, dphi, times)
orbit
f,ax = plt.subplots(1,3, figsize=(20, 5))

ax[0].plot(orbit["x"], orbit["y"])
ax[0].set_xlabel('x (mas)')
ax[0].set_ylabel('y (mas)')

ax[1].plot(orbit["x"], orbit["z"])
ax[1].set_xlabel('x (mas)')
ax[1].set_ylabel('z (mas)')

ax[2].plot(orbit["z"], orbit["y"])
ax[2].set_xlabel('z (mas)')
ax[2].set_ylabel('y (mas)')

for aax in ax:
    aax.set_aspect('equal')
    aax.set_xlim(-200, 200)
    aax.set_ylim(-200, 200)
output = orbit
output['foo1'] = np.linspace(0, len(orbit.index), len(orbit.index), dtype=int)
output['foo2'] = np.linspace(0, len(orbit.index), len(orbit.index), dtype=int)
output['foo3'] = np.linspace(0, len(orbit.index), len(orbit.index), dtype=int)
output['foo4'] = np.linspace(0, len(orbit.index), len(orbit.index), dtype=int)
output['foo5'] = np.linspace(0, len(orbit.index), len(orbit.index), dtype=int)

output.to_csv("S2orbit.raw", sep=' ', index=False, header=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20237132848288872, 'fraction_numerical': 0.05739692805173808, 'mean_word_length': 3.4885126964933493, 'pattern_counts': {'":': 4, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '69475', 'license': 'MIT', 'n_tokens_mistral': 1689, 'n_tokens_neox': 1564, 'n_words': 328, 'path': 'rawdata/calculateS2orbit.ipynb', 'repo_name': 'ageller/UniviewS2Precession', 'size': '3711'}