two_balls_create_pi.jl

two_balls_create_pi.jl

using LinearAlgebra

using Plots
using ProgressMeter: @showprogress

Base.@kwdef mutable struct Ball
    x::Float64
    v::Float64
    r::Float64 = 1.0
    m::Float64 = 1.0
end

Base.@kwdef struct Wall
    xmin::Float64 = -10
    xmax::Float64 = 10
    ground::Float64 = -1
end

@recipe function f(w::Wall)
    seriescolor := :black
    label := :none
    @series begin
        [w.xmin, w.xmax], [w.ground, w.ground]
    end
    @series begin
        [w.xmax, w.xmax], [w.ground, 5.0]
    end
end

@recipe function f(b::Ball)
    seriestype := :shape
    linecolor --> :black
    aspect_ratio --> :equal
    fillalpha --> 0.25
    label --> :none

    θ = 0:0.1:2π
    b.x .+ b.r .* cos.(θ), 0.0 .+ b.r .* sin.(θ)
end

@recipe function f(geomobjs::AbstractArray{<:Ball})
    for obj in geomobjs
        @series begin
            obj
        end
    end
end

@recipe f(geomobjs::Ball...) = collect(geomobjs)

function update!(b::Ball, Δt = 1.0)
    b.x += Δt * b.v
end

function creategif()
    M = 1000
    if M == 100
        nsteps = 800
        Δt = 0.5
    end

    if M == 1000
        nsteps = 3500
        Δt = 0.1
    end

    if M == 10000
        nsteps = 40000
        Δt = 0.01
    end

    if M == 100000
        nsteps = 120000
        Δt = 0.003
    end

    if M == 1000000
        nsteps = 120000
        Δt = 0.0025
    end

    if M == 10000000
        nsteps = 2400000
        Δt = 0.000125
    end
    b2 = Ball(x = -0.5, v = 0.2, r = 1, m = M)
    b1 = Ball(x = 2.5, v = 0, r = 1, m = 1)
    balls = [b2, b1]
    w = Wall(xmax = 20)
    c = 0
    accumΔt = 0.0
    plts = [plot(w)]
    @showprogress for s = 1:nsteps
        accumΔt += Δt
        for b in balls
            update!(b, Δt)
        end
        # これは玉ドン
        if b1.x - b1.r < b2.x + b2.r
            c += 1
            m1 = b1.m
            m2 = b2.m
            v1 = b1.v
            v2 = b2.v
            while (b1.r + b2.r) > abs(b1.x - b2.x)
                b1.x -= 0.001 * v1
                b2.x -= 0.001 * v2
            end
            v1_next = v1 - 2m2 / (m1 + m2) * dot(v1 - v2, b1.x - b2.x) * (b1.x - b2.x) / dot(b1.x - b2.x, b1.x - b2.x)
            v2_next = v2 - 2m1 / (m2 + m1) * dot(v2 - v1, b2.x - b1.x) * (b2.x - b1.x) / dot(b2.x - b1.x, b2.x - b1.x)
            b1.v = v1_next
            b2.v = v2_next
        end

        if b1.x ≥ w.xmax - b1.r
            c += 1
            b1.v *= -1
            b1.x = 2(w.xmax - b1.r) - b1.x
        end
        if accumΔt ≥ 1.0
            p = plot(w)
            plot!(p, balls)
            plot!(title = "M=$M, c=$c, c/sqrt(M) = $(c / sqrt(M))")
            push!(plts, p)
            accumΔt = 0.0
        end
    end
    @show length(plts)
    anim = @animate for p in plts
        plot(p)
    end
    gif(anim, "ball.gif")
end

Screen-Recording-2021-11-24-at-22.33.23.mp4