#!/usr/bin/python

# (C) Michael Puerrer, Sascha Husa

import matplotlib.pyplot as plt
import numpy as np

from math import pi

from scipy.interpolate import UnivariateSpline
from scipy import optimize

from mol_base import *

global h

#################

def advectRHS(u,t):
  [fi, pi] = u

  fiRHS = FD1periodic(fi,h)  # compare boundary conditions for derivative
  piRHS = FD1(fi,h)

  return np.array([fiRHS, piRHS])
  
#################


if __name__ == "__main__":  # now we can import this file as a module without running 'main', or execute it as a script
  
  
# names for the solution components, used for filenames
  solnames = ['fi', 'pi']  # global variable - we are lazy!
  
# grid parameters
  deltax = 0.05
  dt0    = 0.05

  h = deltax

  npoints = 200
  grid = np.arange(npoints) * deltax

  np.savetxt('grid.dat', grid)

# set up initial data
  A     = 1.
  xc    = 5.
  sigma = 1.

  Gaussian = lambda x: A*np.exp(-(x - xc)**2 / sigma)
  u0 = np.array([Gaussian(grid), np.zeros(len(grid))]) # [u0,v0] -- time-symmetric data

  np.savetxt('fi0.dat', u0[0])
  np.savetxt('pi0.dat', u0[1])

  nequations = 2
  tmax = 10
  nstepsmax = 100000
  
# how much to output
  tsample = 1
  xsample = 1
  
  param =[tmax, nstepsmax, nequations, npoints, deltax,tsample,xsample]

 
  dt1 = dt0 # set stepsize
  (tE1, sol1) = ODEint(u0, dt1, param, solnames, RK4Step, advectRHS)

#  dt2 = dt0/2 # set stepsize
#  (tE2, sol2) = ODEint(u0, dt2, param, RK4Step, myODErhs)

#  dt3 = dt0/4 # set stepsize
#  (tE3, sol3) = ODEint(u0, dt3, param, RK4Step, myODErhs)


#  diff12 = sol1                    - (sol2[::2])[:len(sol1)]
#  diff23 = (sol2[::2])[:len(sol1)] - (sol3[::4])[:len(sol1)]

#  ODEsave('diff12', tE1, diff12, '1-2')
#  ODEsave('diff23', tE1, diff23, '2-3')