Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
@terasakisatoshi
Copy link
Author

terasakisatoshi commented Nov 24, 2021

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment