Простейшая задача n тел пишется за 5 минут, она состоит всего из трёх циклов. На более продвинутый метод Barnes-Hut может уйти часа два. Проще сделать самому, если число тел не слишком велико.

Ну 1-3 тыс. квадратных звезд + по 1-10 планет на звезду. И таких 2-4 голактеко.

Это немного, можно даже не особо заморачиваться с оптимизациями. Правда у меня есть ощущение, что моделировать одновременно и звёзды и планеты не выйдет: точности флотов может не хватить.

А какая цель? Чтобы получилось похоже, или чтобы получилось как следует?

Эх =( Чтоб епично! А если серьезно, то для начала хоть как нибудь, ибо опыта работы с 3D нет вообще. Для начала и «как следует» будет ок.

Моделировать планеты можно отдельно: на эволюцию системы звёзд они, можно сказать, не влияют (масса на много порядков меньше). Т.е. посчитать сначала положение звёзд, затем перейти в с.к. с меньшим масштабом и посчитать движение планет в поле ближайших звёзд.

Disclaimer: я в этом вопросе диванный теоретик, и с такой задачей не сталкивался.

Считать хочется слабо, поэтому и ищу движок в котором гравитация уже есть — останется всё только настроить.

В игровых движках того, о чём я писал в предыдущем посте, точно нет. Не думаю даже, что там и кулоновские взаимодействия-то есть: слишком специфичная задача. В BGE (blender game engine?) используется физ. движок bullet, емнип. Попробуй погуглить «bullet engine n body», «bullet engine coulomb» и т.д.

попробуй вот это - может поможет осилить.

tt=0;
close all;

M0=[1e16 3e14 1e14]';
X0=     [0 0;
        100 50;
        200 50]';

function retval = pot(x, t)

        M0=[1e16 3e14 1e14]';
        X0=     [0 0;
                100 50;
                200 50]';

        G=6.67e-11;

        retval=[0 0]';
        for l=1:columns(X0)
                retval+=G*M0(l)/norm(x-X0(:,l))^2*(X0(:,l)-x)/norm(x-X0(:,l));
        end
end

#x''=f(x,t), x'(t0)=v_0, x(t0)=x_0; здесь x - это радиус-вектор ракеты.
#
#раскладываем в систему уравнений:
#
#w: это [x, u, y, v]
#x'=u; u'=производная f по x;
#y'=v; v'=производная f по y;

function retval = f(w, t)
        retval=[];

        retval(2)=pot(w([1 3]), t)(1);
        retval(4)=pot(w([1 3]), t)(2);
        retval(1)=w(2);
        retval(3)=w(4);
end


#начальная скорость
v_0=[-10 30]';

#начальное положение
x_0=[200 0]';

w0=[x_0(1) v_0(1) x_0(2) v_0(2)]';

t=0:0.005:30;

z=lsode("f", w0, t);

hold on;
plot(z(:,1), z(:,3));
plot(X0(1,:), X0(2,:), 'bo');
plot(x_0(1), x_0(2), 'r^');

гугл рулит

http://en.wikipedia.org/wiki/Physics_engine#Real-time_physics_engines

Ну? Т. е. BGE сможет справиться с задачей?

Всем благодарностей, буду пробовать.

ты плохо описал свою задачу.

изучи возможности движков и посмотри, подходят ли они под нее.

Оптимально он справиться с задачей не сможет. Большое количество активных объектов очень сильно влияет на производительность Bullet (даже без проверки столкновений).

Самый быстрый вариант — решение задачи n-body, как уже подсказали.