# import some packages
from __future__ import print_function
import numpy as np

########################
# define some functions:
########################

# right deriv
def Dp(u, dx, du):
    for i in range(0, len(u)-1):
        du[i] = (u[i+1] - u[i])/dx 

# left deriv
def Dm(u, dx, du):
    for i in range(1, len(u)):
        du[i] = (u[i] - u[i-1])/dx 

# centered deriv
def D0(u, dx, du):
    for i in range(1, len(u)-1):
        du[i] = (u[i+1] - u[i-1])/(2.0*dx) 

# RHS of $\partial_t u$
def eval_rhs(u, dx, du):
    D0(u, dx, du)
    return -du

# u(t+dt) = u(t) + \partial_t u * dt
def calc_unew(u, rhs, dt):
    return u + rhs * dt

# print columns with data at time t
def pr_timeframe(t, x, u):
    print("# time =", t)
    for i in range(0, len(u)):
        print(x[i], u[i])
    print()

###############
# main program:
###############

# grid: 11 points from 0 to 1, i.e. x[0]=0, ..., x[10]=1
x = np.linspace(0, 1, 11)
dx = x[1]-x[0]  # grid spacing
dt = 0.5*dx     # time step

# du will contain deriv of u, initialze to 0
du = np.zeros(len(x))

# initial u
t = 0.0
u = np.sin(x)
pr_timeframe(t, x, u)

timesteps = 20

# loop over time steps, and print results
for n in range(0, timesteps):
    # make time step, using simple Euler method
    t = t + dt
    rhs = eval_rhs(u, dx, du)
    u = calc_unew(u, rhs, dt)
    # print columns with data
    pr_timeframe(t, x, u)