import numpy as np
from numpy import sin, cos, tan, pi
from pprint import pprint
import sys

from tklib.tkutils import replaceif


# numerical differentiation by two-point forward difference method
def diff2forward(func, x, h):
    return (func(x+h) - func(x)) / h;

def diff2backward(func, x, h):
    return (func(x) - func(x-h)) / h;

def diff3(func, x, h):
    return (func(x+h) - func(x-h)) / 2.0 / h;

def diff5(func, x, h):
    return (-func(x+2.0*h) + 8.0*func(x+h) - 8.0*func(x-h) + func(x-2.0*h)) / 12.0 / h;

def diff7(func, x, h):
    return (func(x+3.0*h) - 9.0*func(x+2.0*h) + 45.0*func(x+h) 
            - 45.0*func(x-h) + 9.0*func(x-2.0*h) - func(x-3.0*h)) / 60.0 / h;

def diff_richardson(func, x, h0):
# create 2D array by specifing the size to be nhmax * nhmax
    D = []
    for i in range(nhmax):
        D.append([])
        for j in range(nhmax):
            D[i].append([])

    D[0][0] = (func(x0+h0) - func(x0-h0)) / 2.0 / h0
    print("D{}{}: {}".format(0, 0, D[0][0]))

    for m in range(1, nhmax+1):
        h0 *= 0.5
        D[0][m] = (func(x0+h0) - func(x0-h0)) / 2.0 / h0
        print("D0{}: {}".format(m, D[0][m]))

        for k in range(m-1, -1, -1):
            m1 = m - k
            D[m1][k] = (4.0**m1 * D[m1-1][k+1] - D[m1-1][k]) / (4.0**m1 - 1)
            print("D{}{}: {}".format(m1, k, D[m1][k]))
        diff = D[m][0] - D[m-1][0]
        print("  diff = D[{}][0] - D[{}][0] = {} - {} = {}".format(m, m-1, D[m][0], D[m-1][0], diff))

        if abs(diff) < eps:
            return (D[m][0] + D[m-1][0]) /2.0
    return (D[m][0] + D[m-1][0]) /2.0

