from IPython.display import display
from sympy import init_printing
from sympy import symbols, as_finite_diff, solve, latex
from sympy import Function, Eq

fg, f0, f1, f2 = symbols('f_g, f_0, f_1, f_2')
z, h = symbols('z, h')
a, b = symbols('a, b')
f = Function('f')

init_printing()

"""
Explanation: Derivation of the cauchy BC
We would like to derive the cauchy BC, which for a field $f$ reads
$$
f(0) = a\
\partial_z f\big|_0 = b
$$
using a fourth order approximation for the derivative.
WARNING: This scheme is only first order convergent
End of explanation
"""


extraPolate = Eq(fg, 16*a/5 - 3*f0 + f1 - f2/5)
display(extraPolate)

"""
Explanation: Extrapolation of $f(0) = a$ to the ghost point yields (see ghost4thOrder for calculation) yields
End of explanation
"""


eq1 = Eq(0, extraPolate.rhs - extraPolate.lhs)
display(eq1)

"""
Explanation: Which can be rewritten to
End of explanation
"""


deriv = as_finite_diff(f(z).diff(z), [z-h/2, z+h/2, z+3*h/2, z+5*h/2])
deriv = Eq(b ,deriv.subs([(f(z-h/2),  fg),\
                          (f(z+h/2),   f0),\
                          (f(z+3*h/2), f1),\
                          (f(z+5*h/2), f2),\
                         ]).together())
display(deriv)

"""
Explanation: Furthermore a 4th order FD of $\partial_z f\big|_0 = b$ reads
End of explanation
"""


eq2 = Eq(0, deriv.rhs - deriv.lhs)
display(eq2)

"""
Explanation: Which can be rewritten to
End of explanation
"""


full = Eq(eq1.rhs, eq2.rhs)
display(full)

fullSolvedForFg = Eq(fg, solve(full, fg)[0].collect(symbols('f_0, f_1, f_2, h'), exact=True).simplify())
display(fullSolvedForFg)
print(fullSolvedForFg)

print(latex(fullSolvedForFg))

"""
Explanation: Thus
End of explanation
"""


