monte_carlo_numpy_matplotlib_square_circle.py

# Heitetään pisteitä satunnaisesti neliölle
# Osa pisteistä osuu neliön sisällä olevalle neljännesympyrälle
# Tulostetaan heitot
# Lasketaan tuloksista Piin arvo
# Juhani Kaukoranta 29.3.2019
import numpy as np
import matplotlib.pyplot as plt
import time

# heitetään satunnaispisteitä ympyrän jonka keskellä on neliö
# neliön koko 400*400 pikseliä
# neliön vasen ylänurkka (100,500)
# neliö oikea alanurkka (500,500)
# keskipiste (300,300)

def monte_carlo_pi(n):
    # pisteiden x- ja y-koordinaatit satunnaisvektoreina
    # ja x- ja ykoordinaattien neliöt
    x= np.random.randint(100,501,n)
    y= np.random.randint(100,501,n)
    # circle-ehdon täyttävät pisteet ympyrän sisällä
    circle = (x - 300)**2 + (y - 300)**2 <= 40000
    nocircle = (x - 300)**2 + (y - 300)**2 > 40000
    xout = x[nocircle]
    yout = y[nocircle]
    xcircle = x[circle]
    ycircle = y[circle]
    # Pii = 4 * pisteet ympyrän sisällä/kaikki pisteet
    sisalla = np.sum(circle)
    ulkona = n - sisalla
    Pii = 4.0*sisalla/n
    # plottaus
    plt.plot(xcircle, ycircle,'r.')
    plt.plot(xout,yout,'b.')
    plt.axis([100, 500, 100, 500])
    plt.xlabel("Satunnaispisteitä "+str(n)+", sisällä "+str(sisalla)+" ja ulkona "+str(ulkona))
    plt.text(250,300, 'Pii '+str(Pii),family="DejaVu Sans",size="20" )
    plt.show()
    # accum sisältää ympyräneljännekseen osuvien pisteiden lukumäärän
    # Piin arvon laskeminen
    return Pii

pisteita = int(input("Kuinka monta pistettä arvotaan "))
time0 = time.perf_counter()
print("Pii = ",monte_carlo_pi(pisteita))
time1 = time.perf_counter()
print("Aikaa kului ",time1-time0," sekuntia laskentaan, piirtämiseen ja kuvan katsomiseen")