import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Определение пределов переменновй величины
t = np.arange(0, 24, 0.1)

# Создание дифференциальной функции изменяемой величины
def victoria_rega(s, t):
    dsdt = k * E_sun / (np.sqrt(np.pi)) \
         * s * np.sqrt(s) * np.cos(np.pi/12 * (t - 12))
    return dsdt

# Определение начальных условий и параметров
s_0 = 1600 # Начальная площадь
E_sun = 1367 # Вт / м^2
k = 340 * 10**(-8) # Подбираемый коэффициент

# Решение поставленной задачи
solve = odeint(victoria_rega, s_0, t)

# Построение решения в виде графика функции
plt.plot(t, solve[:,0])
plt.plot([0, 24],[2500, 2500], label='Максимальная площадь')
plt.xlabel('Время, часы')
plt.ylabel('Площадь листа, см^2')
plt.title('Виктория-рега')
plt.legend()
plt.grid()

plt.show()