#!/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 ode_base import *

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

def myODErhs(u,t):
  [x,y] = u

  XRHS = -x
  YRHS = -y

  return np.array([XRHS,YRHS])
  
#################


if __name__ == "__main__":  # now we can import this file as a module without running 'main', or execute it as a script
  
  nequations = 2
  tmax       = 100
  nstepsmax  = 100000
  param = [tmax, nstepsmax, nequations]

  u0 = [1.,1.]
  dt0 = 2.785
    
  dt1 = dt0 # set stepsize
  (tE1, sol1) = ODEint(u0, dt1, param, RK4Step, myODErhs)
  print 'stopped at n = %d, t = %.4f, state vector = ' %(len(sol1), tE1) #, sol1[-1])   
  ODEsave('dt1', tE1, sol1, ['x','y'])
  ODEplot(tE1, sol1, ['x','y'])

  dt2 = dt0/2 # set stepsize
  (tE2, sol2) = ODEint(u0, dt2, param, RK4Step, myODErhs)
  print 'stopped at n = %d, t = %.4f, state vector = ' %(len(sol2), tE2) #, sol2[-1])   
  ODEsave('dt2', tE2, sol2, ['x','y'])
  ODEplot(tE2, sol2, ['x','y'])

  dt3 = dt0/4 # set stepsize
  (tE3, sol3) = ODEint(u0, dt3, param, RK4Step, myODErhs)
  print 'stopped at n = %d, t = %.4f, state vector = ' %(len(sol3), tE3) #, sol3[-1])   
  ODEsave('dt3', tE3, sol3, ['x','y'])
  ODEplot(tE3, sol3, ['x','y'])


  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')