import numpy as np
import matplotlib.pyplot as plt

# Set up figure font sizes for presentation
import matplotlib.pylab as pylab
from scipy.interpolate import interp1d

params = {'legend.fontsize': 22,
          'figure.figsize': (15, 10),
         'axes.labelsize': 22,
         'axes.titlesize':30,
          'lines.markersize':15,
         'xtick.labelsize':22,
         'ytick.labelsize':22}
pylab.rcParams.update(params)


# Set up the time grid
t_final = 10.0    # end point
h       = 0.1    # time step
t_vec   = np.arange(0, t_final, h) # create evenly spaced points
t_vec   = np.append(t_vec, t_vec[-1] + h)  # add endpoint
n = len(t_vec) - 1  # number of intervals, or Euler steps

A = np.zeros(n+1)
P = np.zeros(n)
Q = np.zeros(n)

# interpoleingeksempel
P2 = [1, 2, 0, 2, 1, 0, 0, 0.5, 0.4, 2, 0 ]
t2 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
P2_interp = interp1d(t2, P2, kind = 'cubic')


A[0] = 0.0
p = 1.0
K = 1.0
Amin = 0.5

for k in range(n):
    P[k] = 0.0
    if t_vec[k] < 1.0:
        P[k] = p
    Q[k] = K * max(A[k] - Amin, 0.0)

#    A[k+1] = A[k] + h * (P[k] - Q[k] )
    A[k+1] = A[k] + h * (P2_interp(t_vec[k]) - Q[k] )


plt.plot(t_vec, A, 'o-')
plt.plot(t_vec[0:n], P2_interp(t_vec[0:n]))
plt.plot(t_vec[0:n], Q)

plt.show()

VS = h * np.sum(P - Q)
HS = A[-1] - A[0]
print VS - HS