%% Berechnung von Pi
% Mittels Monte-Carlo Verfahren
% Idee geht auf die Phaenomena 1984 in Zuerich zuerueck

% Vers. 1.2 05.2016 (03.2001)
% Copyright W.Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016

%% Visualisierung
fh1=figure;                % Figure-Umgebung oeffnen
%set(fh1,'Position',[100 100 1060 850])
hold on
posax=get(gca,'position'); % Position der Achse auslesen
                           % Rechteck plotten
anrec=annotation('rectangle',posax,...
    'FaceAlpha',0.2,'FaceColor','red','EdgeColor','red');
                           % Ellipse plotten
anell=annotation('ellipse',posax);
halt=1;k=0;
                           % Button zum anhalten der Berechnung plotten
h = uicontrol('Style', 'pushbutton', 'String', 'halt',...
    'Units','Normalized','Position', [0.9 0.05 0.08 0.05], 'Callback', 'halt=0;');
xlim([-1,1]),ylim([-1,1])  % Groesse der Achsen

%% Zufallsberechnung ausfuehren
zupi=2*rand(100000,2)-1;
drin=zupi(:,1).^2+zupi(:,2).^2<1; % Radius^2 berechnen
anzahl=[1:length(drin)]';
% Testen ob der zufaellig gewaehlte Punkt innerhalb des Kreises liegt 
while halt && k<anzahl(end)
    k=k+1;
    if drin(k)
        plot(zupi(k,1),zupi(k,2),'.g','MarkerSize',6)
    else
        plot(zupi(k,1),zupi(k,2),'.r','MarkerSize',6)
    end
    pause(0.11)
end
%% Abschaetzung
pimontecarlo=4*cumsum(drin)./anzahl;
ah=axes('Position',[0.2 0.2 0.6 0.7]);
plot(pimontecarlo)
ylim([3 3.25])
set(ah,'YGrid','on')
hold off
