RosettaCodeData/Task/Total-circles-area/MATLAB/total-circles-area.m
2015-02-20 09:02:09 -05:00

45 lines
1.6 KiB
Matlab

function res = circles()
tic
%
% Size of my grid -- higher values => higher accuracy.
%
ngrid = 5000;
xc = [1.6417233788 -1.4944608174 0.6110294452 0.3844862411 -0.2495892950 1.7813504266 -0.1985249206 -1.7011985145 -0.4319462812 0.2178372997 -0.6294854565 1.7952608455 1.4168575317 1.4637371396 -0.5263668798 -1.2197352481 -0.1389358881 1.5293954595 -0.5258728625 -0.1403562064 0.8055826339 -0.6311979224 1.4685857879 -0.6855727502 0.0152957411];
yc = [1.6121789534 1.2077959613 -0.6907087527 0.2923344616 -0.3832854473 1.6178237031 -0.8343333301 -0.1263820964 1.4104420482 -0.9499557344 -1.3078893852 0.6281269104 1.0683357171 0.9463877418 1.7315156631 0.9144146579 0.1092805780 0.0030278255 1.3782633069 0.2437382535 -0.0482092025 0.7184578971 -0.8347049536 1.6465021616 0.0638919221];
r = [0.0848270516 1.1039549836 0.9089162485 0.2375743054 1.0845181219 0.8162655711 0.0538864941 0.4776976918 0.7886291537 0.0357871187 0.7653357688 0.2727652452 1.1016025378 1.1846214562 1.4428514068 1.0727263474 0.7350208828 1.2472867347 1.3495508831 1.3804956588 0.3327165165 0.2491045282 1.3670667538 1.0593087096 0.9771215985];
r2 = r .* r;
ncircles = length(xc);
%
% Compute the bounding box of the circles.
%
xmin = min(xc-r);
xmax = max(xc+r);
ymin = min(yc-r);
ymax = max(yc+r);
%
% Keep a counter.
%
inside = 0;
%
% For every point in my grid.
%
for x = linspace(xmin,xmax,ngrid)
for y = linspace(ymin,ymax,ngrid)
if any(r2 > (x - xc).^2 + (y - yc).^2)
inside = inside + 1;
end
end
end
box_area = (xmax-xmin) * (ymax-ymin);
res = box_area * inside / ngrid^2;
toc
end