Anonymous user
Total circles area: Difference between revisions
→{{header|Haskell}}: analytical solution
(→{{header|Haskell}}: analytical solution) |
|||
Line 326:
{{out}}
<pre>Approximated area: 21.564955642878786</pre>
Analytical solution, breaking down circles to non-intersecting arcs and assemble zero winding paths, then calculate their areas. Pro: precision doesn't depend on a step size, so no need to wait longer for a more precise result; Con: probably not numerically stable in marginal situations, which can be catastrophic.
<lang haskell>import Data.List (sort)
vcross (a,b) (c,d) = a*d - b*c
vdot (a,b) (c,d) = a*c + b*d
vadd (a,b) (c,d) = (a+c, b+d)
vsub (a,b) (c,d) = (a-c, b-d)
vminus (a,b) (c,d) = (a-c, b-d)
vlen2 x = vdot x x
vlen x = sqrt $ vlen2 x
vdist a b = vlen (a `vminus` b)
vscale s (x,y) = (x*s, y*s)
vnorm (x,y) = (x/l, y/l) where l = vlen (x,y)
vangle (x,y) = atan2 y x
anorm a | a > pi = a - pi * 2
| a < (-pi) = a + pi * 2
| otherwise = a
circle_cross (x0,y0,r0) (x1,y1,r1)
| d >= r0 + r1 || d <= abs(r0 - r1) = []
| otherwise = map anorm [ang - da, ang + da]
where
d = vdist (x0,y0) (x1,y1)
s = (r0 + r1 + d)/2
a = sqrt $ s * (s - d) * (s - r0) * (s - r1)
h = 2 * a / d
dr = (x1 - x0, y1 - y0)
dx = vscale (sqrt $ r0^2 - h^2) $ vnorm dr
ang = if r0^2 + d^2 > r1^2
then vangle dr
else pi + (vangle dr)
da = asin (h / r0)
arc_point (x,y,r) a = vadd (x,y) (r * cos a, r * sin a)
arc_start (c,(a0,a1)) = arc_point c a0
arc_mid (c,(a0,a1)) = arc_point c ((a0+a1)/2)
arc_end (c,(a0,a1)) = arc_point c a1
arc_center((x,y,r),_) = (x,y)
arc_area ((_,_,r),(a0,a1)) = r * r * (a1 - a0) / 2
split_circles cs = filter (not.in_any_circle) arcs where
arcs = concatMap csplit circs
circs = map f cs where
f c = (c, sort $ [-pi,pi] ++ (concatMap (circle_cross c) cs))
csplit (c, angs) = zip (repeat c) $ zip angs $ tail angs
in_circle ((x0,y0),c) (x,y,r) = c /= (x,y,r) && vdist (x0,y0) (x,y) < r
in_any_circle arc = or $ map (in_circle (arc_mid arc, c)) cs
where (c,_) = arc
make_paths arcs = join_arcs [] arcs where
join_arcs a [] = [a]
join_arcs [] (x:xs) = join_arcs [x] xs
join_arcs a (x:xs)
| vdist (arc_start (head a)) (arc_end (last a)) < 1e-3
= a : join_arcs [] (x:xs)
| vdist (arc_end (last a)) (arc_start x) < 1e-3
= join_arcs (a++[x]) xs
| otherwise = join_arcs a (xs++[x])
path_area arcs = area_ 0 [] arcs where
area_ a e [] = a + polyline_area e
area_ a e (arc:as) = area_ (a + arc_area arc) (e++[arc_center arc, arc_end arc]) as
polyline_area [a,b] = 0
polyline_area (a:b:c:rest) = tri_area a b c + polyline_area (a:c:rest) where
tri_area a b c = (b `vsub` a) `vcross` (c `vsub` b) / 2
circles_area = sum . map path_area . make_paths . split_circles
main = print $ circles_area [
( 1.6417233788, 1.6121789534, 0.0848270516),
(-1.4944608174, 1.2077959613, 1.1039549836),
( 0.6110294452, -0.6907087527, 0.9089162485),
( 0.3844862411, 0.2923344616, 0.2375743054),
(-0.2495892950, -0.3832854473, 1.0845181219),
( 1.7813504266, 1.6178237031, 0.8162655711),
(-0.1985249206, -0.8343333301, 0.0538864941),
(-1.7011985145, -0.1263820964, 0.4776976918),
(-0.4319462812, 1.4104420482, 0.7886291537),
( 0.2178372997, -0.9499557344, 0.0357871187),
(-0.6294854565, -1.3078893852, 0.7653357688),
( 1.7952608455, 0.6281269104, 0.2727652452),
( 1.4168575317, 1.0683357171, 1.1016025378),
( 1.4637371396, 0.9463877418, 1.1846214562),
(-0.5263668798, 1.7315156631, 1.4428514068),
(-1.2197352481, 0.9144146579, 1.0727263474),
(-0.1389358881, 0.1092805780, 0.7350208828),
( 1.5293954595, 0.0030278255, 1.2472867347),
(-0.5258728625, 1.3782633069, 1.3495508831),
(-0.1403562064, 0.2437382535, 1.3804956588),
( 0.8055826339, -0.0482092025, 0.3327165165),
(-0.6311979224, 0.7184578971, 0.2491045282),
( 1.4685857879, -0.8347049536, 1.3670667538),
(-0.6855727502, 1.6465021616, 1.0593087096),
( 0.0152957411, 0.0638919221, 0.9771215985)]</lang>
{{out}}
21.565036603856395
=={{header|Perl 6}}==
| |||