{{ message }}

Instantly share code, notes, and snippets.

# terasakisatoshi/two_balls_create_pi.jl

Created Nov 24, 2021
two_balls_create_pi.jl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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 commented Nov 24, 2021

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