import numpy as np
import matplotlib.pyplot as plt

r_X = 1
r_Y = 0.95
K_X = 1e6

n_hours = 24
initial_population = 1000
pop_X = np.zeros(n_hours + 1)
pop_Y = np.zeros(n_hours + 1)

pop_X[0] = initial_population
for i in range(n_hours):
    pop_X[i + 1] = pop_X[i] + r_X * (1 - (pop_X[i] / K_X)) * pop_X[i]
    
print("Population of species X:", pop_X)

pop_Y[0] = initial_population
for i in range(n_hours):
    pop_Y[i + 1] = pop_Y[i] + r_Y * (1 - (pop_Y[i] / K_X)) * pop_Y[i]

plt.figure(figsize=(6,3))
plt.plot(pop_X)

plt.xlabel("time (hours)")
plt.ylabel("population")
plt.title("Species X")

plt.figure(figsize=(6,3))
plt.plot(pop_Y)

plt.xlabel("time (hours)")
plt.ylabel("population")
plt.title("Species Y")

data_X = np.loadtxt("data_exp_X.txt")
data_Y = np.loadtxt("data_exp_Y.txt")

plt.figure(figsize=(6,3))
plt.plot(pop_X, label="model")
plt.plot(data_X, label="experiment")

plt.xlabel("time (hours)")
plt.ylabel("population")
plt.title("Species X")
plt.legend()

plt.figure(figsize=(6,3))
plt.plot(pop_Y, label="model")
plt.plot(data_X, label="experiment")

plt.xlabel("time (hours)")
plt.ylabel("population")
plt.title("Species Y")
plt.legend()

plt.figure(figsize=(6,3))
plt.plot(data_X, color='#ff7f0e')