import math


# def RK4_step(f, y, t, dt, N):
#     dt /= N
#     for k in range(N):
#         k1 = f(y, t) * dt
#         k2 = f(y + k1 / 2, t + dt / 2) * dt
#         k3 = f(y + k2 / 2, t + dt / 2) * dt
#         k4 = f(y + k3, t + dt) * dt
#         y, t = y + (k1 + 2 * (k2 + k3) + k4) / 6, t + dt
#     return y

def RK4_step(f, y, t, dt, N):
    temp_y = [y]
    temp_t = [t]
    for k in range(N):
        k1 = f(y, t) * dt
        k2 = f(y + k1 / 2, t + dt / 2) * dt
        k3 = f(y + k2 / 2, t + dt / 2) * dt
        k4 = f(y + k3, t + dt) * dt
        y, t = y + (k1 + 2 * (k2 + k3) + k4) / 6, t + dt
        temp_y.append(y + (k1 + 2 * (k2 + k3) + k4) / 6)
        temp_t.append(t + dt)
    return temp_y, temp_t


def AB_3(mass_y, mass_t, f, h, N):
    y = 0
    while mass_y[-1] >= 40.0:
        y = mass_y[-1] + h * (23 / 12 * f(mass_y[-1], mass_t[-1]) - 4 / 3 * f(mass_y[-2], mass_t[
            -2]) + 5 / 12 * f(mass_y[-3], mass_t[-3]))
        mass_y.append(y)
        mass_t.append(mass_t[-1] + h)
    return mass_t[-1]


def AB_4(mass_y, mass_t, f, h, N):
    y = 0
    while mass_y[-1] >= 40.0:
        y = mass_y[-1] + h * (55 / 24 * f(mass_y[-1], mass_t[-1]) - 59 / 24 * f(mass_y[-2], mass_t[
            -2]) + 37 / 24 * f(mass_y[-3], mass_t[-3]) - 3 / 8 * f(mass_y[-4], mass_t[-4]))
        mass_y.append(y)
        mass_t.append(mass_t[-1] + h)
    return mass_t[-1]


def AB_5(mass_y, mass_t, f, h, N):
    y = 0
    while mass_y[-1] >= 40.0:
        y = mass_y[-1] + h * (
                1901 / 720 * f(mass_y[-1], mass_t[-1]) - 1387 / 360 * f(mass_y[-2], mass_t[
            -2]) + 109 / 30 * f(mass_y[-3], mass_t[-3]) - 637 / 360 * f(mass_y[-4], mass_t[
            -4]) + 251 / 720 * f(mass_y[-5], mass_t[-5]))
        mass_y.append(y)
        mass_t.append(mass_t[-1] + h)
    return mass_t[-1]


def euler(f, y0, a, b, h):
    t, y = a, y0
    while t <= b:
        print(t, y)
        t += h
        y += h * f(t, y)


def func(T, time):
    return -1 / 10 * math.log(3) * (T - 25)


# euler(func, 70, 0, 10, 0.000001)
mass_y, mass_t = RK4_step(func, 70, 0, 0.001, 2)
print(AB_3(mass_y, mass_t, func, 0.001, 10))

mass_y, mass_t = RK4_step(func, 70, 0, 0.001, 3)
print(AB_4(mass_y, mass_t, func, 0.001, 10))

mass_y, mass_t = RK4_step(func, 70, 0, 0.001, 4)
print(AB_5(mass_y, mass_t, func, 0.001, 10))