¿Cómo puedo calcular los períodos orbitales en un sistema estelar binario?

Si necesita una estimación aproximada, puede muestrear las posiciones de los objetos estelares a partir de su aceleración, utilizando la ley de Newton. Se dibuja una imagen completa en esta página de Wikipedia , pero básicamente, dados N objetos estelares en la posición P (i), con su respectiva masa M (i):

$$M_i.\vec{acc_i} = -G\sum_{n \neq i}^N \frac {M_i.M_n.(\vec{pos_i}-\vec{pos_n})}{\begin{vmatrix} \vec{pos_i}-\vec{pos_n}\end{vmatrix}^3 }$$

Luego puede derivar la aceleración en velocidad, luego en posición usando un delta de tiempo lo suficientemente pequeño. Con las posiciones iniciales y la velocidad (que son las más complicadas y divertidas de elegir), puede simular un sistema de N-cuerpo aproximado.

En el lado de la diversión y la simulación, esto es lo que pretende presentar Universe Sandbox (de antemano, perdón por vincular a una tienda comercial, no estoy relacionado con este estudio).

Los puntos de Lagrangian también son interesantes de observar al simular.

Sin embargo, para simplificar las cosas, puede pensar que casi todos los objetos estelares "cercanos" al sistema binario habrían sido devorados por la pareja masiva a lo largo del tiempo, y luego considerar solo el baricentro del par de estrellas.

EDITAR: aunque OP fue claro, mi respuesta es para N objetos. La simplificación a N=2, siendo A y B las posiciones de los dos objetos, da como resultado:

$$\vec{acc_A} = \frac{G.M_B}{AB^3}\vec{AB}$$ $$\vec{acc_B} = \frac{G.M_A}{BA^3}\vec{BA}$$ Y con un proceso iterativo, una vez calculada la aceleración para el paso k, usando condiciones iniciales conocidas para velocidad y posición: $$\vec{spd_{k+1}} = \vec{acc_k}.\Delta{t} + \vec{spd_k}$$ $$\vec{pos_{k+1}} = \vec{spd_k}.\Delta{t} + \vec{pos_k}$$

EDIT2: bueno, otra estimación aproximada, para la duración del período (que malinterpreté como "trayectoria" de la pregunta). Usaría coordenadas cilíndricas como referencia, para su representación espacial limpia y para la conversión de un solo ángulo:

$$\left( \begin{array}{c} R_k \\ \theta_k \\ h_k \end{array} \right) = \left( \begin{array}{c} \sqrt{pos_{x,k}^2+pos_{y,k}^2} \\ \arctan \left( \frac{pos_{y,k}}{pos_{x,k}} \right) \\ pos_{z,k} \end{array} \right)$$

...con la debida precaución. Introduce parte de esa conversión en cada iteración y espera a que theta complete una rotación.

(De antemano, perdón por la notación extraña, que solo traté de hacer consistente con mis publicaciones anteriores)