Se precisar de uma estimativa aproximada, você pode amostrar as posições dos objetos estelares a partir de sua aceleração, usando a lei de Newton. Uma imagem completa é desenhada nesta página da Wikipedia, mas basicamente, dados N objetos estelares na posição P (i), com suas respectivas massas 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} $$

Você pode então derivar a aceleração em velocidade e depois em posição usando um delta de tempo pequeno o suficiente. Com posições iniciais e velocidade (que são as mais complicadas e também mais engraçadas de escolher), você pode simular um sistema de corpo N bruto.

No lado da diversão e da simulação, é isso que Universe Sandbox visa apresentar (antecipadamente, desculpe por me vincular a uma loja comercial, não sou relacionado a este estúdio).

Os pontos lagrangianos também são interessantes de se olhar durante a simulação.

No entanto, para simplificar as coisas, você pode pensar que quase todos os objetos estelares "próximos" do sistema binário teriam sido comidos pelo casal massivo ao longo do tempo, e então considere apenas o baricentro do par de estrelas.

EDITAR: embora OP tenha sido claro, minha resposta é para N objetos. A simplificação para N = 2, com A e B sendo as posições dos dois objetos, resulta em: $$ \ vec {acc_A} = \ frac {G.M_B} {AB ^ 3} \ vec {AB} $$$ $ \ vec {acc_B} = \ frac {G.M_A} {BA ^ 3} \ vec {BA} $$ E com um processo iterativo, uma vez que a aceleração é calculada para a etapa k, usando condições iniciais conhecidas para velocidade e posição: $ $ \ 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: bem, outra estimativa grosseira, para a duração do período (que eu entendi mal como "trajetória" da questão ) Eu usaria coordenadas cilíndricas como referência, para sua representação espacial limpa e para a conversão de ângulo único:

$$ \ 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) $$

... com precaução adequada. Apresente parte dessa conversão a cada iteração e espere que theta complete uma rotação.

(antecipadamente, desculpe pela notação estranha, que eu tentei tornar consistente em minhas postagens anteriores)