Орбиты частиц

Пример показывает бесконечную 2D-анимацию: 100 частиц движутся вокруг центрального тяжёлого тела по законам ньютоновской гравитации. Для каждой частицы сохраняются последние 20 положений, поэтому за ней остаётся короткий полупрозрачный след.

Разбор

Идея

Центральное тело находится в точке (CX, CY). Для каждой частицы выбирается случайный радиус и угол, а начальная скорость направляется по касательной к окружности:

radius := 65.0 + 215.0 * rnd(1.0)
angle := 2.0 * PI * rnd(1.0)
speed := sqrt(MU / radius) * speedScale

x[i] := CX + radius * cos(angle)
y[i] := CY + radius * sin(angle)
vx[i] := -speed * sin(angle)
vy[i] := speed * cos(angle)

speedScale немного меняет скорость, поэтому орбиты не выглядят одинаковыми.

Шаг физики

На каждом кадре вычисляется ускорение к центру:

dx := x[i] - CX
dy := y[i] - CY
r2 := dx * dx + dy * dy
invR := 1.0 / sqrt(r2)
invR3 := invR * invR * invR
ax := -MU * dx * invR3
ay := -MU * dy * invR3

Затем скорость и координаты обновляются простым численным шагом:

vx[i] := vx[i] + ax * DT
vy[i] := vy[i] + ay * DT
x[i] := x[i] + vx[i] * DT
y[i] := y[i] + vy[i] * DT

Следы

Массивы tx и ty хранят историю последних TAIL положений каждой частицы. Перед записью новой позиции старые значения сдвигаются:

нц для k от TAIL - 1 до 1 шаг -1
    tx[i, k] := tx[i, k - 1]
    ty[i, k] := ty[i, k - 1]
кц
tx[i, 0] := x[i]
ty[i, 0] := y[i]

При отрисовке дальние точки хвоста рисуются слабее, а текущая позиция частицы — ярче.

Кадры

Программа рисует весь кадр на текущем листе, а затем вызывает:

новый лист(W, H, чёрный)

В онлайн-песочнице для Рисователя это граница кадра: уже нарисованный лист показывается в браузере, после чего подготавливается новый чёрный лист для следующего кадра. Поэтому бесконечный цикл не блокирует страницу и виден как анимация.

Структура программы

Пример разделён на три алгоритма:

инициализировать — задаёт начальные координаты, скорости и хвосты частиц
шаг_физики — пересчитывает координаты и обновляет историю хвостов
нарисовать — рисует центральное тело, хвосты, частицы и завершает кадр командой новый лист

Полная программа

| 100 частиц на орбитах вокруг центрального тяжелого тела.
| Каждый кадр рисуется целиком, а "новый лист" в конце кадра
| показывает его и подготавливает следующий черный лист.

использовать Рисователь

цел W = 900
цел H = 700
цел N = 100
цел TAIL = 20

вещ CX = 450.0
вещ CY = 350.0
вещ MU = 2400.0
вещ DT = 0.35
вещ PI = 3.141592653589793

алг
нач
    вещ таб x[0:N-1], y[0:N-1], vx[0:N-1], vy[0:N-1]
    вещ таб tx[0:N-1, 0:TAIL-1], ty[0:N-1, 0:TAIL-1]

    цел frame

    инициализировать(x, y, vx, vy, tx, ty)
    новый лист(W, H, чёрный)
    frame := 0

    нц пока да
        шаг_физики(x, y, vx, vy, tx, ty)
        нарисовать(frame, x, y, tx, ty)
        frame := frame + 1
    кц
кон

алг инициализировать(арг рез вещ таб x[0:N-1], арг рез вещ таб y[0:N-1], арг рез вещ таб vx[0:N-1], арг рез вещ таб vy[0:N-1], арг рез вещ таб tx[0:N-1, 0:TAIL-1], арг рез вещ таб ty[0:N-1, 0:TAIL-1])
нач
    цел i, k
    вещ radius, angle, speed, speedScale

    нц для i от 0 до N - 1
        radius := 65.0 + 215.0 * rnd(1.0)
        angle := 2.0 * PI * rnd(1.0)
        speedScale := 0.94 + 0.12 * rnd(1.0)
        speed := sqrt(MU / radius) * speedScale

        x[i] := CX + radius * cos(angle)
        y[i] := CY + radius * sin(angle)
        vx[i] := -speed * sin(angle)
        vy[i] := speed * cos(angle)

        нц для k от 0 до TAIL - 1
            tx[i, k] := x[i]
            ty[i, k] := y[i]
        кц
    кц
кон

алг шаг_физики(арг рез вещ таб x[0:N-1], арг рез вещ таб y[0:N-1], арг рез вещ таб vx[0:N-1], арг рез вещ таб vy[0:N-1], арг рез вещ таб tx[0:N-1, 0:TAIL-1], арг рез вещ таб ty[0:N-1, 0:TAIL-1])
нач
    цел i, k
    вещ dx, dy, r2, invR, invR3, ax, ay

    нц для i от 0 до N - 1
        dx := x[i] - CX
        dy := y[i] - CY
        r2 := dx * dx + dy * dy
        invR := 1.0 / sqrt(r2)
        invR3 := invR * invR * invR
        ax := -MU * dx * invR3
        ay := -MU * dy * invR3

        vx[i] := vx[i] + ax * DT
        vy[i] := vy[i] + ay * DT
        x[i] := x[i] + vx[i] * DT
        y[i] := y[i] + vy[i] * DT

        нц для k от TAIL - 1 до 1 шаг -1
            tx[i, k] := tx[i, k - 1]
            ty[i, k] := ty[i, k - 1]
        кц
        tx[i, 0] := x[i]
        ty[i, 0] := y[i]
    кц
кон

алг нарисовать(цел frame, вещ таб x[0:N-1], вещ таб y[0:N-1], вещ таб tx[0:N-1, 0:TAIL-1], вещ таб ty[0:N-1, 0:TAIL-1])
нач
    цел i, k
    цел hue, alpha

    перо(1, RGB(255, 230, 130))
    кисть(RGB(255, 230, 130))
    окружность(int(CX), int(CY), 12)

    нц для i от 0 до N - 1
        hue := mod(div(i * 360, N) + div(frame, 4), 360)
        нц для k от TAIL - 1 до 1 шаг -1
            alpha := 16 + (TAIL - k) * 8
            перо(1, HSLA(hue, 85, 60, alpha))
            кисть(HSLA(hue, 85, 60, alpha))
            окружность(int(tx[i, k]), int(ty[i, k]), 1)
        кц

        перо(1, HSL(hue, 90, 65))
        кисть(HSL(hue, 90, 65))
        окружность(int(x[i]), int(y[i]), 2)
    кц

    новый лист(W, H, чёрный)
кон

▶ Запустить пример