import numpy as np
import matplotlib.pyplot as plt
from cycler import cycler

from plot_helpers import CO, bold_ticks


def fig7_4():
    # euler method for various M values
    # y' = f(t,y) with y(0) = y0 '

    ### this requires the subaxis and parseArgs files

    r = 10
    t0 = 0
    y0 = 0.01
    tmax = 5

    # calculate exact solution
    tt = np.linspace(t0, tmax, 100)
    a0 = (1 - y0) / y0
    exact = np.zeros(100)
    for it in range(100):
        exact[it] = 1 / (1 + a0 * np.exp(-r * tt[it]))

    plt.close("all")
    # get(gcf)
    fig, axes = plt.subplots(2, 1, figsize=(8.8, 5.0))
    fig.subplots_adjust(left=0.08, right=0.99, bottom=0.12, top=0.99, hspace=0.18)
    for ax in axes:
        ax.set_prop_cycle(cycler(color=CO))

    ax = axes[0]
    step = 3 / r
    m = round(tmax / step)
    t = np.linspace(t0, tmax, m + 1)
    k = t[1] - t[0]
    kr = k * r
    print(kr)
    y_euler = eulers(t, y0, k, m + 1)
    ax.plot(t, y_euler, "--bo", markersize=7, linewidth=1.5)
    ax.set_xlim(0, tmax)
    ax.set_ylim(0, 1.5)
    ax.plot(tt, exact, "r", linewidth=1.5)
    ax.text(2.1, 0.35, "kr = 3", fontsize=14, fontweight="bold")
    ax.set_xticks([0, 1, 2, 3, 4, 5])
    ax.grid(True)
    ax.set_ylabel("y-axis")
    bold_ticks(ax, 14)
    for spine in ax.spines.values():
        spine.set_visible(True)

    ax = axes[1]
    step = 1 / r
    m = round(tmax / step)
    t = np.linspace(t0, tmax, m + 1)
    k = t[1] - t[0]
    kr = k * r
    print(kr)
    y_euler = eulers(t, y0, k, m + 1)
    ax.plot(t, y_euler, "--bo", markersize=7, linewidth=1.5)
    ax.plot(tt, exact, "r", linewidth=1.5)
    ax.text(2.1, 0.3, "kr = 1", fontsize=14, fontweight="bold")
    ax.set_xlim(0, tmax)
    ax.set_ylim(0, 1.1)
    ax.set_xticks([0, 1, 2, 3, 4, 5])
    ax.grid(True)
    ax.set_xlabel("t-axis")
    ax.set_ylabel("y-axis")
    bold_ticks(ax, 14)
    for spine in ax.spines.values():
        spine.set_visible(True)

    plt.show()

    # plt.savefig('/Users/mark/Desktop/eulerS.eps')


# right-hand side of DE
def f(t, y):
    r = 10
    z = r * y * (1 - y)
    return z


# euler method
def eulers(t, y0, k, n):
    y = y0
    ypoints = np.empty(n)
    ypoints[0] = y0
    for i in range(1, n):
        yy = y + k * f(t[i - 1], y)
        ypoints[i] = yy
        y = yy
    return ypoints


if __name__ == "__main__":
    fig7_4()
