Самоучитель бильярда: Книга: «Бильярд. Самоучитель» — Николай Здобников. Купить книгу, читать рецензии | ISBN 985-6751-56-X

Система и алгоритм

Бильярд состоит из трех частей:

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

Важно, чтобы эти три вещи, «частицы», «препятствия» и «процессы» оставались независимыми . Это позволит расширять наш код!

5

.

Алгоритм

Для имитации бильярда нужно следовать очень простому алгоритму. Предполагая, что частица p уже инициализирована в бильярде, шаги

• Вычислить столкновения p со всеми препятствиями в бильярде
• Найти первое по времени столкновение и препятствие, которому оно соответствует
• Промыть и повторить!

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

Эта независимость возможна только при многократной отправке. Без него, так или иначе, этот метод был бы «принадлежит» либо к «классу» Препятствий , либо к «классу» Частиц . Однако это помешало бы расширяемости! Отсутствие необходимости «присваивать» процесс определенному «классу» — огромное преимущество Джулии.

6

.

Базовая структура s нам понадобится

Все векторы, которые мы будем использовать, двумерны. Таким образом, мы можем воспользоваться Статические массивы . Для упрощения я определю SV как 2D-вектор Float (все векторы, которые мы будем использовать, будут этого типа)

 с использованием StaticArrays. 
const SV = SVector{2,Float64} 

Теперь мы определяем частицу struct . Для расширяемости лучше всего определить его на основе абстрактного типа.

 абстрактный тип AbstractParticle конец

изменяемая структура Particle <: AbstractParticle
    поз::СВ
    вел::СВ
конец
Частица(x0, y0, φ0) = Particle(SV(x0, y0), SV(cos(φ0), sin(φ0))) 

Затем мы определяем препятствие struct s. Опять же, чтобы обеспечить расширяемость, абстракцию и обобщение, лучше всего определить все препятствия как подтип абстрактного типа.

 абстрактный тип Конец препятствия

структура Стена <: Препятствие
    сп::СВ
    эп::СВ
    нормальный :: СВ
конец

структура диска <: препятствие
    c::СВ
    р::Float64
end 

"Бильярд" в дальнейшем будет просто Tuple из Obstacle подтипов. В частности, тип бильярда NTuple{N, Препятствие} где {N} , поэтому удобно определить:

 const Billiard = NTuple{N, Препятствие} где N 

7

.

Столкновения

Как было сказано во введении, наиболее важной частью является функция, которая находит столкновения между частицами и препятствиями. Вот как это выглядит для столкновения стандартной частицы со стеной:

 с помощью LinearAlgebra: точка, нормализация

"""
    столкновение(p::AbstractParticle, o::Obstacle) → t, cp
Найдите столкновение (если есть) между данной частицей и препятствием.
Возвращает время до столкновения и предполагаемую точку столкновения `cp`.
"""
Коллизия встроенной функции @(p::Particle, w::Wall)
    n = normalvec(w, p.pos)
    denom = точка (p.vel, n)
    если номинал ≥ 0,0
        возврат без столкновений()
    еще
        t = точка (w. sp - p.pos, n)/denom
        вернуть t, p.pos + t * p.vel
    конец
конец

normalvec(w::Стена, позиция) = w.normal 

а вот та же функция, но вместо коллизий с дисками:

 @inline function Collision(p::Particle, d::Disk)
    dotp = точка (p.vel, normalvec (d, p.pos))
    dotp ≥ 0.0 && вернуть отсутствие столкновений()

    dc = p.pos - d.c
    B = dot(p.vel, dc) # указывает на центр круга: B < 0
    C = dot(dc, dc) - d.r*d.r #нахождение вне круга: C > 0
    Δ = В*В - С

    Δ ≤ 0,0 && вернуть без столкновений()
    sqrtD = sqrt(Δ)
    # Ближайшая точка:
    t = -B - sqrtD
    вернуть t, p.pos + t * p.vel
конец

normalvec(d::Disk, pos) = normalize(pos - d.c) 

Вы можете видеть, что есть случаи, когда коллизии невозможны или они происходят в обратном направлении во времени. По соглашению, это значение, которое мы тогда возвращаем:

 @inline nocollision() = (Inf, SV(0. 0, 0.0)) 

next_collision — полезная функция, которая находит «истинное» следующее столкновение. Он просто обходит препятствия в бильярде. Он просто проверяет, какое препятствие имеет наименьшее время столкновения:

 function next_collision(p::AbstractParticle, bd)
    j, ct, cp = 0, Inf, SV(0,0, 0,0)
    для i в каждом индексе (bd)
        t, c = столкновение (p, bd [i])
        если t < ct
            дж = я
            кт = т
            ср = с
        конец
    конец
    вернуть j, ct, ср
конец 

8

.

Эволюция частицы в бильярде

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

 распространение!(p::Particle, pos, t) = (p.pos = pos) 

Это встроенная функция (примечание ! в конце)

Мы также должны определить функцию, которая выполняет зеркальное отражение, т.е. изменяет скорость частицы (после столкновения)

 function resolvecollision!(p::AbstractParticle, o::Obstacle)
    n = normalvec(o, p.pos)
    p.vel = p.vel - 2*точка(n, p.vel)*n
end 

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

Для удобства стоит определить следующую функцию:

 """
    подпрыгнуть! (р, бд)
Эволюционируйте частицу за одно столкновение (на месте).
"""
@inline function bounce!(p::AbstractParticle, bd)
    i::Int, tmin::Float64, cp::SV = next_collision(p, bd)
    если tмин != Inf
        распространять!(p, cp, tmin)
        разрешить столкновение!(p, bd[i])
    конец
    возврат i, tmin, p. pos, p.vel
end 

Затем мы можем использовать его внутри более крупной функции, которая вызывает bounce! до указанного времени:

 """
    временной ряд!(p::AbstractParticle, bd, n) -> xt, yt, t
Эволюционируйте частицу в бильярде `bd` для `n` столкновений
и вернуть временной ряд позиций `xt, yt` вместе с вектором времени `t`.
"""
временной ряд функций!(p::AbstractParticle, bd, n::Int)

    т = [0,0]; xt = [p.pos[1]]; yt = [p.pos[2]]; с = 0

    в то время как c < n
        предыдущая = п.поз.; prevvel = p.vel
        i, ct = отскок!(p, bd)
        xs, ys = экстраполировать (p, prevpos, prevvel, ct)
        толчок! (т, кт)
        добавить! (xt, xs)
        добавить!(yt, ys)
        с += 1
    конец

    вернуть xt, yt, t
конец 

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

 экстраполировать (p::Particle, prevpos, prevvel, ct) = p. pos 

Почему существует эта функция экстраполировать ? Смотрите ниже, когда мы расширяем наш код для магнитных частиц!

9

.

Запустив код

давайте проверим это прямо сейчас! Мы создадим знаменитый Синайский бильярд

 х, у, г = 1,0, 1,0, 0,3
зр = [0,0,у]; эп = [0,0, 0,0]; п = [х, 0,0]
leftw = Стена (sp, ep, n)
зр = [х, 0,0]; эп = [х, у]; п = [-х, 0,0]
rightw = Стена (sp, ep, n)
зр = [х, у]; эп = [0,0, у]; п = [0,0,-у]
topw = Стена (sp, ep, n)
сп = [0,0,0,0]; эп = [х, 0,0]; п = [0,0, у]
botw = Стена (sp, ep, n)
диск = Диск([x/2, y/2], r)

bd = (botw, rightw, topw, leftw, disk)
bd isa Billiard 

, а также инициализировать частицу

 p = Particle(0.1, 0.1, 2π*rand()) 

и развить ее

 xt, yt, t = временные ряды!(p, bd, 10) 

10

.

Построение графика

Давайте определим несколько простых методов построения графика и построим результат!

 с использованием PyPlot
импортировать PyPlot: сюжет

const EDGECOLOR = (0,0.6,0)
график функции (d :: Disk)
    цвет лица = (КРАЙЦВЕТ..., 0,5)
    круг1 = PyPlot.plt[:Круг](d.c, d.r;
        edgecolor = EDGECOLOR, facecolor = facecolor, lw = 2.0)
    PyPlot.gca()[:add_artist](круг1)
конец
график функции (w :: Стена)
    PyPlot.plot([w.sp[1],w.ep[1]],[w.sp[2],w.ep[2]];
        цвет = ЦВЕТ КРАЯ, lw = 2,0)
конец
График функции (bd:: Billiard)
    для o ∈ bd; участок (о); конец
    gca()[:set_aspect]("равно")
конец

фигура(); участок (бд)
gcf() 

Отлично! Теперь посмотрим на орбиту

 figure(); участок (бд)
xt, yt, t = временные ряды!(p, bd, 10)
сюжет (хт, ут)
gcf() 

Постройте 10 орбит!

 цифра(); участок (бд)
p = частица (0,1, 0,5, 2π*rand())
для j в 1:10
    xt, yt = временные ряды!(p, bd, 20)
    график (xt, yt, альфа = 0,5)
конец
gcf() 

11

.

Витрина 1: производительность и метапрограммирование

 с использованием BenchmarkTools

p = частица (0,1, 0,1, 2π*rand())
@btime bounce!($p, $bd) 

Уже очень быстро разнести частицу за одно столкновение, однако есть некоторые выделения (хотя функция на месте).

Эти выделения происходят из-за нестабильности типа в next_collision , поскольку кортеж содержит элементы разных типов. Однако с помощью метапрограммирования легко решить эту нестабильность типа, потому что Кортеж имеет известный размер!

В следующем определении мы используем метапрограммирование, чтобы «развернуть» цикл

 @сгенерированная функция next_collision(p::AbstractParticle, bd::Billiard)
    L = length(bd.types) # обратите внимание, что здесь bd означает ТИП bd!
    выход = :(ind = 0; tmin = Inf; cp = SV(0,0, 0,0))
    для j=1:L
        push!(out. args, цитата
                            пусть х = bd[$j]
                                tcol, pcol = столкновение (p, x)
                                # Установить минимальное время:
                                если tкол < tмин
                                  tмин = tкол
                                  индекс = $j
                                  ср = пкол
                                конец
                            конец
                        конец
                        )
    конец
    push!(out.args, :(return ind, tmin, cp))
    вернуться
конец


@btime bounce!($p, $bd) 

Этот номер безумный !!! Обратите внимание, что этот код не зависит от бильярда! Вы можете пройти любой набор препятствий, и это все равно будет так же эффективно!!! Время отказов! масштабируется линейно с количеством препятствий в бильярде.

12

.

Витрина 2: Расширяемость

Допустим, мы хотим добавить еще одно препятствие к этому «бильярдному набору», который мы делаем. Мы должны переписать все для него? Неа! В конце концов нам нужно расширить только два метода! Только два !

Чтобы показать это, давайте создадим эллипс как препятствие с полуосями a, b

 struct Ellipse <: Obstacle
    c::СВ
    а::поплавок64
    б::поплавок64
end 

Нам нужно расширить следующие методы:

 normalvec
столкновение 

Да!!! Только два! Итак, приступим! normalvec довольно просто:

 function normalvec(e::Ellipse, pos)
    х₀, у₀ = поз.
    ч, к = эк
    return normalize(SV((x₀-h)/(e.a*e.a), (y₀-k)/(e.b*e.b)))
конец

используя LinearAlgebra: норма
функция столкновения (p::Particle, e::Ellipse)
    dotp = точка (p. vel, normalvec (e, p.pos))
    dotp ≥ 0.0 && вернуть отсутствие столкновений()

    а = э.а; б = е.б
    pc = p.pos - ec
    μ = p.vel[2]/p.vel[1]
    ψ = пк[2] - мк*пк[1]

    номинал = а*а*μ*μ + b*b
    Δ² = номинал - ψ*ψ
    Δ² ≤ 0 && возврат без столкновений()
    Δ = квадрат (Δ²); f1 = -a*a*µ*ψ; f2 = b*b*ψ # просто множители
    I1 = SV(f1 + a*b*Δ, f2 + a*b*μ*Δ)/деномин
    I2 = SV(f1 - a*b*Δ, f2 - a*b*µ*Δ)/деномин

    d1 = норма (pc - I1); d2 = норма (pc - I2)
    вернуть d1 < d2? (d1, I1 + д.к) : (d2, I2 + д.к)
конец 

Хорошо, теперь давайте создадим биллиард с эллипсом и диском, ради удовольствия

 el = Ellipse([0.4, 0.2 ], 0.3, 0.1)
ди = Диск ([0,6, 0,7], 0,25)

bd2 = Billiard((bd[1:4]..., el, di)) 

и построить его

 function plot(e::Ellipse)
    цвет лица = (КРАЙЦВЕТ..., 0,5)
    эллипс = PyPlot.matplotlib[:patches][:Ellipse](e.c, 2e.a, 2e.b;
        edgecolor = EDGECOLOR, facecolor = facecolor, lw = 2. 0)
    PyPlot.gca()[:add_artist](эллипс)
конец

фигура(); сюжет(bd2)
gcf() 

Теперь мы готовы создать частицу в этом совершенно новом бильярде:

 p = Particle(0.1, 0.1, 2π*rand())
xt, yt, t = временной ряд!(p, bd2, 20)
фигура(); сюжет(bd2)
сюжет (хт, ут)
gcf() 

нарисуйте еще кучу!

 цифра(); сюжет(bd2)
для j в 1:10
    p = частица (0,1, 0,1, 2π*rand())
    xt, yt = временной ряд!(p, bd2, 20)
    график (xt, yt, альфа = 0,5)
конец
gcf() 

13

.

Витрина 3: Снова расширяемость.

Итак, оказалось почти тривиально добавить дополнительное препятствие в наш код. А как насчет лишней частицы?

Я не говорю о еще одном экземпляре Частицы . Я говорю о новой частице типа , которая движется по-другому.

В этой части мы создадим этот новый тип, MagneticParticle , который движется по кругу, а не по прямой линии! Но сколько функций нам нужно определить? Если вы уже определили тип MagneticParticle , тогда столько:

 столкновения # для каждого препятствия, которое мы хотим поддерживать
распространять!
экстраполировать 

и да, вот и все. Может быть трудно поверить, что это занимает так мало, но это правда!!!

13.1

.

Тип

 изменчивая структура MagneticParticle <: AbstractParticle
    поз::СВ
    вел::СВ
    ω :: Float64
конец
MagneticParticle(x0, y0, φ0, ω) = MagneticParticle(SV(x0, y0), SV(cos(φ0), sin(φ0)), ω) 

Эта частица движется по окружности с угловой скоростью ω .

13,2

.

Продление столкновения

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

Вот столкновение со стеной: 92 -4*а*в Δ ≤ 0,0 && вернуть без столкновений() u1 = (-b - sqrt(Δ))/2a u2 = (-b + sqrt(∆))/2a cond1 = 0,0 ≤ u1 ≤ 1,0 услов2 = 0,0 ≤ и2 ≤ 1,0 θ, I = без столкновений () если условие1 || условие2 dw = w. ep - w.sp для (u, cond) в ((u1, cond1), (u2, cond2)) Y = w.sp + u*dw если условие φ = реальный угол (p, w, Y) φ < θ && (θ = φ; I = Y) конец конец конец возврат θ*pr, I конец

и вот столкновение с диском: 92) I1 = СВ( pc[1] + a*(p1[1] - pc[1])/d + h*(p1[2] - pc[2])/d, pc[2] + a*(p1[2] - pc[2])/d - h*(p1[1] - pc[1])/d ) I2 = СВ( pc[1] + a*(p1[1] - pc[1])/d - h*(p1[2] - pc[2])/d, pc[2] + a*(p1[2] - pc[2])/d + h*(p1[1] - pc[1])/d ) θ1 = реальный угол (p, o, I1) θ2 = реальный угол (p, o, I2) вернуть θ1 < θ2 ? (θ1*rc, I1) : (θ2*rc, I2) end

Функции cyclotron и realangle являются вспомогательными функциями. Первый находит центр и радиус циклотрона, описываемого частицей.

 циклотрон(p) = (p. 2))
    d2r > 2 && (d2r = 2,0)
    θпростое = acos (1,0 - d2r)

    ПИ = я - Р0
    сторона = (PI[1]*PC[2] - PI[2]*PC[1])*ω

    сторона < 0 && (θпростое = 2π-θпростое)
    возврат
конец