import numpy as np
import matplotlib.pyplot as plt #Will be used to create the relevant plots

#Parameters
t_end = 150

Na0 = 1e4
Nb0 = 5e3
Nc0 = 0

V = 5

k_f = 1e-3
k_r = 5e-2

#packing in to vectors
k = [k_f, k_r]
N0 = np.array([Na0, Nb0, Nc0])

def rates(c):
    r1 = k[0]*c[0]*c[1]
    r2 = k[1]*c[2]

    return np.array([r1, r2])

#Execution code
eqil = False
t = 0
N = N0
N_list = [[x for x in N]]
t_list = [t]
while not eqil and t < t_end:
    c = N / V
    r = rates(c)
    p = r/sum(r)

    dN = N_list[9 * len(N_list)//10][0] - N_list[-1][0] #Change in Na over last 10% of steps
    dt = t_list[9 * len(t_list)//10] - t_list[-1] #Change in t over last 10% of steps
    dNdt = dN/dt #average over last 10% of steps

    if len(t_list) > 100 and abs(dNdt) < 1:
        eqil = True

    else:
        ran = (1 - 1e-16)*np.random.random() + 1e-16 #sikrer at 1 er med og at ikke 0 er med
        if ran <= p[0]:
            N += np.array([-1, -1, 1])
        else:
            N += np.array([1, 1, -1])

        t += -np.log(ran)/sum(r)

        N_list.append([x for x in N])
        t_list.append(t)

#Results and plotting
if t_list[-1] >= t_end:
    print('Forced stop at t =', t_end)
    print('N :', N)
    print('V :', V)
    print('t :', t)

Na, Nb, Nc = np.array(N_list).transpose()
plt.plot(t_list, Na/V, label=r'$N_a$')
plt.plot(t_list, Nb/V, label=r'$N_b$')
plt.plot(t_list, Nc/V, label=r'$N_c$')
plt.ylabel('c')
plt.xlabel('t')
plt.legend()
plt.show()