#!/usr/bin/env python
"""
Try an analytic solution to Burgers equation.

author: Bob Wimmer
date: January 1, 2013
email: wimmer@physik.uni-kiel.de

"""

a = 1.
mu = 0.1

from scipy.integrate import quad
from numpy import inf as inf
from numpy import exp, linspace


def phi(xi):
    """
    define the initial condition. Currently, this returns a ramp.
    """
    if xi < 0.: r = 0.
    elif xi < 1.: r= 1.-xi/1.
    else: r= 0.
    return r

def G(theta,x,t):
    tmp1 = ((x-theta)**2)/2./t
    tmp2,aux = quad(phi,0.,theta)
    return tmp1 + tmp2

def denom_integrand(theta,x,t):
    return exp(-(a*G(theta,x,t)/(2*mu)))

def num_integrand(theta,x,t):
    return (x-theta)/t*exp(-(a*G(theta,x,t)/(2*mu)))

def u(x,t):
    tmp1 = quad(denom_integrand,-inf,inf,args=(x,t))
    tmp2 = quad(num_integrand,-inf,inf,args=(x,t))
    return tmp2[0]/tmp1[0]

#this allows us to use this file as a module, but also to call it as a standalone script
if __name__ == "__main__":
    t = 0.0
    for x in linspace(0.,1.,11):
        print 'x, t, u(x,t) = ', x, t, phi(x)
    print '---------------------------------'
    t = 0.001
    for x in linspace(0.,1.,11):
        print 'x, t, u(x,t) = ', x, t, u(x,t)
    print '---------------------------------'
    t = 0.1
    for x in linspace(0.,1.,11):
        print 'x, t, u(x,t) = ', x, t, u(x,t)
    print '---------------------------------'
    t = 0.2
    for x in linspace(0.,1.,11):
        print 'x, t, u(x,t) = ', x, t, u(x,t)
    print '---------------------------------'