73 lines
1.7 KiB
Julia
73 lines
1.7 KiB
Julia
include("utils/point.jl")
|
|
|
|
# list of things we want to iterate over
|
|
domain = Iterators.filter(
|
|
a -> a[3] > a[2] && a[2] > a[1] && gcd(collect(a)) == 1,
|
|
Iterators.product(1:100, 1:100, 1:100)
|
|
)
|
|
|
|
# law of cosines, determines angle based on the side lengths
|
|
inv_law_cosines(a, b, c) = acos((a * a + b * b - c * c) / (2.0 * a * b))
|
|
|
|
function circumcenter(A, B, C)
|
|
magA = magnitude(A)
|
|
magB = magnitude(B)
|
|
magC = magnitude(C)
|
|
D = 2.0 * (A.x * (B.y - C.y) + B.x * (C.y - A.y) + C.x * (A.y - B.y))
|
|
Ux = (magA * (B.y - C.y) + magB * (C.y - A.y) + magC * (A.y - B.y)) / D
|
|
Uy = (magA * (C.x - B.x) + magB * (A.x - C.x) + magC * (B.x - A.x)) / D
|
|
Point(Ux, Uy)
|
|
end
|
|
|
|
function lerp(src, dst, dist)
|
|
total = distance(src, dst)
|
|
ratio = dist / total
|
|
dx = src.x + (dst.x - src.x) * ratio
|
|
dy = src.y + (dst.y - src.y) * ratio
|
|
Point(dx, dy)
|
|
end
|
|
|
|
function compute_d(A, B, C, r1, r2, r3)
|
|
ABmid = lerp(A, B, r1)
|
|
BCmid = lerp(B, C, r2)
|
|
CAmid = lerp(C, A, r3)
|
|
D = circumcenter(ABmid, BCmid, CAmid)
|
|
|
|
k1, k2, k3 = 1.0 / r1, 1.0 / r2, 1.0 / r3
|
|
s = k1 + k2 + k3
|
|
s2 = k1 * k2 + k2 * k3 + k1 * k3
|
|
e1 = s + 2.0 * sqrt(s2)
|
|
e2 = s - 2.0 * sqrt(s2)
|
|
re = 1.0 / max(e1, e2)
|
|
|
|
EAB = inv_law_cosines(r1 + re, r1 + r2, r2 + re)
|
|
Ex = (r1 + re) * cos(EAB)
|
|
Ey = (r1 + re) * sin(EAB)
|
|
E = Point(Ex, Ey)
|
|
|
|
distance(D, E)
|
|
end
|
|
|
|
function calc(args)
|
|
r1 = args[1]
|
|
r2 = args[2]
|
|
r3 = args[3]
|
|
|
|
a = r2 + r3
|
|
b = r1 + r3
|
|
c = r1 + r2
|
|
|
|
A = Point(0.0, 0.0)
|
|
B = Point(Float64(c), 0.0)
|
|
|
|
CAB = inv_law_cosines(b, c, a)
|
|
Cx = b * cos(CAB)
|
|
Cy = b * sin(CAB)
|
|
C = Point(Cx, Cy)
|
|
|
|
compute_d(A, B, C, r1, r2, r3)
|
|
end
|
|
|
|
l = map(calc, domain)
|
|
print(round(sum(l) / length(l), digits=8))
|