import matplotlib.pyplot as plt
import numpy as np
import random
import scipy.special as ss

def plot(N = 10000, b = 2, y_max = 2 / (np.sqrt(np.pi)) ):
    def erf_func(x):
        return ((2 / np.sqrt(np.pi)) * np.exp(-x ** 2))


    Integral = np.zeros(N)
    N_I = 0

    for i in range(N):
        # (x,y) er et tilfeldig punkt innenfor kvadratet
        # (0,b) x (0,y_max)
        x = b * random.random()  # kanskje gange med b her, da det er randbetingelsen
        y = y_max * random.random()

        if y < erf_func(x):
            N_I += 1

        # totalt areal er (y_max * b)
        # N_I er totalt antall treff innenfor integralet
        # (i+1) er antall forsøk
        I = (y_max * b) * (N_I / (i + 1))
        Integral[i] = I

    plt.plot([i + 1 for i in range(N)], Integral, color="green", label = "Calculated value of erf("+str(b)+")")
    plt.plot([0, N], [ss.erf(b), ss.erf(b)], label = "Exact value of erf("+str(b)+")")
    plt.ylabel("Value of errorfunction")
    plt.xlabel("Number of trials")
    plt.legend()
    plt.ylim(ss.erf(b)-(b/25), ss.erf(b) + (b/25))
    plt.title("Calculated value of erf("+str(b)+") using Monte Carlo method vs. the exact value")