<p>I am helping my high schooler with a research paper that uses <a href="https://pyastronomy.readthedocs.io/en/latest/" rel="nofollow noreferrer">PyAstronomy</a> to fit the orbit of stars around Sag A. It isn't working because the angles (parameters Omega, omega, and i) seem to have a different reference point. It does work for Star 2, but when we try to do the same thing for Star 1, the initial inputs from the paper fit terribly.</p>
<p>Here are the initial inputs:</p>
<p>i = 119.14 degrees<br />
Omega = 342.04 degrees<br />
omega (w) = 122.3 degrees<br />
a = 595 mas<br />
e = 0.556<br />
tau = 2001.8 years<br />
per = 166 years</p>
<p>She is trying to replicate the results of this paper: An Update on Monitoring Stellar Orbits in the Galactic Center</p>
<p>One difference between the paper and my teen's program is that she is using Python and PyAstronomy, and they used Mathematica. PyAstronomy creates a KeplerEllipseModel object and specifies the axes as xz.</p>
<p>This is the output graph. The blue are the observations and the red is the ellipse created using the values from the paper.</p>
<p><a href="https://i.sstatic.net/8b3yL.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/8b3yL.png" alt="Star 1 orbiting Sag A" /></a></p>
<p>We did a grid search to look at varying the 3 degree parameters (Omega, omega, i), but didn't come close. We also looked at the graphs using the other two axes (xy and yz), also no joy. Does anyone have any ideas on how we can improve this initial fit? The end goal is to run MCMC to sample a^3 against t^2. I hope this makes sense, I am a programmer.</p>
<p>Edit: here is the code</p>
<pre><code>import pandas as pd  # get numbers from csv file
import numpy as np   # Numerical library
import matplotlib.pyplot as plt # Plot commands
from PyAstronomy.modelSuite import KeplerEllipseModel
from PyAstronomy import funcFit as fuf

df = pd.read_csv(&quot;Star1_Orbit_Data_Numbers.csv&quot;)
</code></pre>
<p>Here is the csv file with the data.
<a href="https://drive.google.com/file/d/13KAnsmtyQV5L3ncW_dmW_cdBsfc1lV25/view?usp=sharing" rel="nofollow noreferrer">https://drive.google.com/file/d/13KAnsmtyQV5L3ncW_dmW_cdBsfc1lV25/view?usp=sharing</a></p>
<p>(I double checked this is correct)</p>
<pre><code>def create_arrays(df_data):
  ''' Input is a dataframe with the data for one star.
      Output is 3 numpy arrays: 
        pos - position, alternating RA and DE for each date
        pos_unc - uncertainty for position
        date - date for measurement. 
      Length of pos and pos_unc is 2 * length of date
  '''
  pos = np.empty([2*len(df_data),])
  pos_unc = np.empty([2*len(df_data),])
  date_arr = np.empty([len(df_data),])
  for i in range(len(df_data)):
    pos[2*i] = df_data['RA'][i]
    pos[2*i+1] = df_data['DE'][i]
    pos_unc[2*i] = df_data['RA_error'][i]
    pos_unc[2*i+1] = df_data['DE_error'][i]
    date_arr[i] = df_data['Year'][i]
  return pos, pos_unc, date_arr

def initial_plot(star_num, pos, pos_unc, kem):
  # the data in blue points
  fig = plt.figure(figsize=(10,10))
  plt.title(&quot;S0-&quot; + str(star_num))
  plt.errorbar(pos[0::2], pos[1::2], xerr=pos_unc[0::2], yerr=pos_unc[1::2], fmt=&quot;b.&quot;, markersize=2.0)
  plt.xlabel(&quot;<span class="math-container">$\Delta$</span>RA (arcsec)&quot;)
  plt.ylabel(&quot;<span class="math-container">$\Delta$</span>Dec (arcsec)&quot;)
  plt.axis('equal')
  plt.plot(np.array([0.0]),np.array([0.0]), marker='*', color='black', markersize=5.0)

  #plot the initial guesses in red dashed lines
  guess_model = kem.evaluate(np.linspace(kem[&quot;tau&quot;], kem[&quot;tau&quot;]+kem[&quot;per&quot;], 1000))
  plt.plot(guess_model[0::2], guess_model[1::2], 'r--')
#   fig.savefig(&quot;Star_&quot; + str(star_num) + &quot;_initial_guesses.png&quot;)
  plt.show()

star = 'Star1'
ax1 = 'xz'
kem = KeplerEllipseModel(relevantAxes=ax1)
print(f&quot;Begin initial guesses for grid of {star} {ax1}&quot;)
kem[&quot;a&quot;] = 595
kem[&quot;e&quot;] = 0.556
kem[&quot;tau&quot;] = 2001.8
kem[&quot;per&quot;] = 166
kem[&quot;Omega&quot;] = 342.04
kem[&quot;w&quot;] = 122.3
kem[&quot;i&quot;] = 119.14

pos, pos_unc, date_arr = create_arrays(df)
initial_plot(1, pos, pos_unc, kem)
</code></pre>

CC BY-SA 4.0
2022-10-18T20:48:06.307Z
These might or might not be somewhat helpful and/or fun, I see some orbital elements mentioned in a few places; [Why does the closest approach of star S2 to Sgr A* not appear to be near the focus of its elliptical orbit?](https://astronomy.stackexchange.com/q/35836/7982) and [How did they make a video of the center of the galaxy, and what is it exactly that's flashing there?](https://astronomy.stackexchange.com/q/28940/7982) and [How will they know when to start taking the picture of the black hole at the center of the Milky Way?](https://astronomy.stackexchange.com/q/39354/7982)
uhoh
CC BY-SA 4.0
2022-10-18T20:49:29.423Z
Also, if you'd like to add some relativistic correction to the orbit calculation see answers to [How to calculate the planets and moons beyond Newtons's gravitational force?](https://space.stackexchange.com/q/23408/12102)
uhoh
CC BY-SA 4.0
2022-10-19T04:21:44.400Z
It might help to post the code, there are lots of places there could be issues.
Greg Miller
CC BY-SA 4.0
2022-10-19T10:26:45.360Z
I will post the code in the OP.
Emmy
CC BY-SA 4.0
2022-10-19T12:29:58.347Z
We think we might have found a solution. Rather than running the initial guesses used in the paper, we fixed the parameters except for the 3 angles and ran MCMC on that. There is a bit of an issue because you have to turn off the sanity check on the data inside fitEMCEE and run it for many iterations and a long burnin but it worked for star 1. We are going to do the rest of the stars and fingers crossed.
Emmy
<orbit><orbital-mechanics><python><sagittarius-a>
Kepler ellipse orbit fit issues
Emmy
# Kepler ellipse orbit fit issues

(Asked by: username_0 on 2022-10-18)

I am helping my high schooler with a research paper that uses PyAstronomy to fit the orbit of stars around Sag A. It isn't working because the angles (parameters Omega, omega, and i) seem to have a different reference point. It does work for Star 2, but when we try to do the same thing for Star 1, the initial inputs from the paper fit terribly.
Here are the initial inputs:
i = 119.14 degrees
Omega = 342.04 degrees
omega (w) = 122.3 degrees
a = 595 mas
e = 0.556
tau = 2001.8 years
per = 166 years
She is trying to replicate the results of this paper: An Update on Monitoring Stellar Orbits in the Galactic Center
One difference between the paper and my teen's program is that she is using Python and PyAstronomy, and they used Mathematica. PyAstronomy creates a KeplerEllipseModel object and specifies the axes as xz.
This is the output graph. The blue are the observations and the red is the ellipse created using the values from the paper.

We did a grid search to look at varying the 3 degree parameters (Omega, omega, i), but didn't come close. We also looked at the graphs using the other two axes (xy and yz), also no joy. Does anyone have any ideas on how we can improve this initial fit? The end goal is to run MCMC to sample a^3 against t^2. I hope this makes sense, I am a programmer.
Edit: here is the code
import pandas as pd  # get numbers from csv file
import numpy as np   # Numerical library
import matplotlib.pyplot as plt # Plot commands
from PyAstronomy.modelSuite import KeplerEllipseModel
from PyAstronomy import funcFit as fuf

df = pd.read_csv("Star1_Orbit_Data_Numbers.csv")

Here is the csv file with the data.
https://drive.google.com/file/d/13KAnsmtyQV5L3ncW_dmW_cdBsfc1lV25/view?usp=sharing
(I double checked this is correct)
def create_arrays(df_data):
  ''' Input is a dataframe with the data for one star.
      Output is 3 numpy arrays: 
        pos - position, alternating RA and DE for each date
        pos_unc - uncertainty for position
        date - date for measurement. 
      Length of pos and pos_unc is 2 * length of date
  '''
  pos = np.empty([2*len(df_data),])
  pos_unc = np.empty([2*len(df_data),])
  date_arr = np.empty([len(df_data),])
  for i in range(len(df_data)):
    pos[2*i] = df_data['RA'][i]
    pos[2*i+1] = df_data['DE'][i]
    pos_unc[2*i] = df_data['RA_error'][i]
    pos_unc[2*i+1] = df_data['DE_error'][i]
    date_arr[i] = df_data['Year'][i]
  return pos, pos_unc, date_arr

def initial_plot(star_num, pos, pos_unc, kem):
  # the data in blue points
  fig = plt.figure(figsize=(10,10))
  plt.title("S0-" + str(star_num))
  plt.errorbar(pos[0::2], pos[1::2], xerr=pos_unc[0::2], yerr=pos_unc[1::2], fmt="b.", markersize=2.0)
  plt.xlabel("$\Delta$RA (arcsec)")
  plt.ylabel("$\Delta$Dec (arcsec)")
  plt.axis('equal')
  plt.plot(np.array([0.0]),np.array([0.0]), marker='*', color='black', markersize=5.0)

  #plot the initial guesses in red dashed lines
  guess_model = kem.evaluate(np.linspace(kem["tau"], kem["tau"]+kem["per"], 1000))
  plt.plot(guess_model[0::2], guess_model[1::2], 'r--')
#   fig.savefig("Star_" + str(star_num) + "_initial_guesses.png")
  plt.show()

star = 'Star1'
ax1 = 'xz'
kem = KeplerEllipseModel(relevantAxes=ax1)
print(f"Begin initial guesses for grid of {star} {ax1}")
kem["a"] = 595
kem["e"] = 0.556
kem["tau"] = 2001.8
kem["per"] = 166
kem["Omega"] = 342.04
kem["w"] = 122.3
kem["i"] = 119.14

pos, pos_unc, date_arr = create_arrays(df)
initial_plot(1, pos, pos_unc, kem)

**Tags:** orbit, orbital-mechanics, python, sagittarius-a
> **Comments:**
> **username_1**: These might or might not be somewhat helpful and/or fun, I see some orbital elements mentioned in a few places; [Why does the closest approach of star S2 to Sgr A* not appear to be near the focus of its elliptical orbit?](https://astronomy.stackexchange.com/q/35836/7982) and [How did they make a video of the center of the galaxy, and what is it exactly that's flashing there?](https://astronomy.stackexchange.com/q/28940/7982) and [How will they know when to start taking the picture of the black hole at the center of the Milky Way?](https://astronomy.stackexchange.com/q/39354/7982)  
> **username_1**: Also, if you'd like to add some relativistic correction to the orbit calculation see answers to [How to calculate the planets and moons beyond Newtons's gravitational force?](https://space.stackexchange.com/q/23408/12102)  
> **username_2**: It might help to post the code, there are lots of places there could be issues.  
> **username_0**: I will post the code in the OP.  
> **username_0**: We think we might have found a solution. Rather than running the initial guesses used in the paper, we fixed the parameters except for the 3 angles and ran MCMC on that. There is a bit of an issue because you have to turn off the sanity check on the data inside fitEMCEE and run it for many iterations and a long burnin but it worked for star 1. We are going to do the rest of the stars and fingers crossed.  

---

## 0 answers

