Решение задачи N тел

Итак, приступим к написанию программы.  Для каждой звезды нам известны: масса, положение и скорость в начальный момент времени, также необходимо будет рассчитывать силы (или ускорения), которые также являются векторными величинами. Зададим их.

<span style="color: #0000ff;">#define N 1000</span>
<span style="color: #0000ff;">float mass[N];</span>
<span style="color: #0000ff;">vector_3d s[N];</span>
<span style="color: #0000ff;">vector_3d v[N];</span>
<span style="color: #0000ff;">vector_3d a[N];</span>

Реализацию класса vector_3d мы уже рассмотрели в http://cpblog.ru/algorithms/physical-modeling/clas.html.

Далее в каждый момент времени нам необходимо будет рассчитывать силы, действующие на данную «звезду». Для этого напишем специальную функцию.

<span style="color: #0000ff;">void CalculateAcceleration()</span>
<span style="color: #0000ff;">{</span>
<span style="color: #0000ff;">vector_3d tmp, tmpNorm;</span>
<span style="color: #0000ff;">for(int i = 0 ; i &lt; N ; ++i)</span>
<span style="color: #0000ff;">{</span>
<span style="color: #0000ff;">a[i] = vector_3d(0,0,0);</span>
<span style="color: #0000ff;">for(int j = 0 ; j &lt; N ; ++j)</span>
<span style="color: #0000ff;">if(i!=j)</span>
<span style="color: #0000ff;">{</span>
<span style="color: #0000ff;">tmp = s[j] - s[i];</span>
<span style="color: #0000ff;">tmpNorm = tmp.Normalize();</span>
<span style="color: #0000ff;">a[i]+=tmpNorm* G * mass[j] / (tmp.NormSquared()+0.001f);</span>
<span style="color: #0000ff;">}</span>
<span style="color: #0000ff;">a[i] /= mass[i];</span>
<span style="color: #0000ff;">}</span>
<span style="color: #0000ff;">}</span>

При расчете ускорения  (а точнее  в данной строчке мы рассчитываем именно силу)

<span style="color: #0000ff;">a[i]+=tmpNorm* G * mass[j] / (tmp.NormSquared()+0.01f);</span>

Мы прибавляем в знаменателе 0.001f , чтобы избежать деления на ноль. Так как для расчета ускорения нам необходимо силу, действующую на тело, поделить на массу этого тела, то и получаем

<span style="color: #0000ff;">a[i] /= mass[i];</span>

Теперь собственно функция для рассчета перемещение за промежуток времени t.

<span style="color: #0000ff;">void NBodySim(float t)</span>
<span style="color: #0000ff;">{</span>
<span style="color: #0000ff;">CalculateAcceleration();</span>
<span style="color: #0000ff;">for(int i = 0 ; i &lt; N; ++i)</span>
<span style="color: #0000ff;">{</span>
<span style="color: #0000ff;">v[i] += a[i] * t;</span>
<span style="color: #0000ff;">s[i] += v[i] * t;</span>
<span style="color: #0000ff;">}</span>
<span style="color: #0000ff;">}</span>

Как уже говорилось в, мы решаем задачу численно. В данном случае мы используется алгоритм Эйлера.

Из условия задачи следует, что в начальный момент времени нам известны положение и скорость каждой «звезды». Поэтому предлагаю написать специальную функцию для задания начальных условий.

<span style="color: #0000ff;">void init()</span>
<span style="color: #0000ff;">{</span>
<span style="color: #0000ff;">for(int i = 0 ; i &lt; N; ++i)</span>
<span style="color: #0000ff;">{</span>
<span style="color: #0000ff;">mass[i] = rand();</span>
<span style="color: #0000ff;">s[i].x = rand() / 10000.0f;</span>
<span style="color: #0000ff;">s[i].y = -rand() / 10000.0f;</span>
<span style="color: #0000ff;">s[i].z = rand() / 10000.0f;</span>
<span style="color: #0000ff;">v[i].x = rand() / 10000.0f;</span>
<span style="color: #0000ff;">v[i].y = rand() / 10000.0f;</span>
<span style="color: #0000ff;">v[i].z = -rand() / 10000.0f;</span>
<span style="color: #0000ff;">}</span>
<span style="color: #0000ff;">}</span>

Мы уже рассматривали метод rand()(метод для генерации псевдослучайных чисел) в http://cpblog.ru/algorithms/general-algorithms/linejnyj-kongruentnyj-metod.html .

Почти готово. Далее необходимо решить, как мы будем отображать наши данные. За основу в данном случае взят пример из книги по OpenGL Игоря Тарасова http://opengl.org.ru/books/open_gl/. Скажу лишь только, что в моей реализации «звезды» имеют разный размер в зависимости от массы и разный цвет в зависимости от скорости.

Исходный код можно посмотреть тут NBodySim. Для того, чтобы все работало правильно, необходимо скачать и подключить библиотеки opengl32.lib, glu32.lib, glaux.lib.

2 Коммент.

  1. sveticksmile

    хорошо)))

  2. А промежуток времени t за одну итерацию, это сколько у вас?

Оставить комментарий

Вы должны быть зарегистрированы чтобы комментировать.