# Path: DiscEvolution/FV_interpolation.py
# def compute_centroids(xe, m):
#     '''First order upwind reconstruction
# 
#     args:
#         xe : Cell edge locations
#         m  : Index of radial scaling giving volume, V = x^{m+1} / m+1.
#     '''
#     return ((m + 1) * np.diff(xe**(m+2))) / ((m + 2) * np.diff(xe**(m+1)))
# 
# def construct_FV_edge_weights(xi, m, iL, iR, max_deriv=None, dtype='f8'):
#     '''Solves for the finite-volume interpolation weights.
# 
#     This code follows the methods outlined in Mignone (2014, JCoPh 270 784) to
#     compute the weights needed to reconstruct a function and its derivatives
#     to edges of computational cells. The polynomial is reconstructed to
#     reproduce the averages of the cell and its neighbours.
# 
#     Note that the polynomial computed near the domain edges will be lower 
#     order, which also reduces the number of derivatives available
# 
#     args:
#         xi : locations of cell edges
#         m  : Index of radial scaling giving volume, V = x^{m+1} / m+1.
#         iL : Number of cells to the left of the target cell to use in the 
#              interpolation.
#         iR : Number of cells to the right of the target cell to use in the 
#              interpolation.
#         
#         max_deriv : maximum derivative level to calculate. If not specified,
#                     return all meaningful derivatives.
#     
#         dtype : data type used for calculation, default = 'f8'
# 
#     returns:
#         wp, wm : The weights used for reconstructing the function and its 1st
#                  iL+iR derivatives to the left and right of the cell edges. 
#                  The shape is [len(xi)-1,max_deriv+1, 1+iL+iR]
#     '''
#     # Order of the polynomial
#     order = 1 + iL + iR
# 
#     if max_deriv is None:
#         max_deriv = order - 1
#     elif max_deriv > order - 1:
#         raise ValueError("Maximum derivative must be less than the order of the"
#                          " polynomial fitted")
#     
#     # Setup the beta matrix of Mignone:
#     beta =  _construct_volume_factors(xi, m, order, dtype)
# 
#     # The matrix of extrapolations to the RHS of cell
#     wp = _solve_FV_matrix_weights(xi[1:], iL, iR, beta, max_deriv, dtype)
#         
#     # The matrix of extrapolations to the LHS of cell
#     wm = _solve_FV_matrix_weights(xi[:-1], iR, iL, beta, max_deriv, dtype)
# 
#     return wp, wm


import numpy as np
    import matplotlib.pyplot as plt
from scipy.special import erf
from .FV_interpolation import compute_centroids, construct_FV_edge_weights
        dQ = Q[...,1:] - Q[...,:-1]
        QF, QB = dQ[...,1:], dQ[...,:-1]

        cF, cB = self._cF, self._cB

        num = QF*QB*(cF*QB + cB*QF)
        den = (QB*QB + (cF + cB - 2)*QB*QF + QF*QF) + 1e-300

        dQ_lim = np.where(QB*QF > 0, num/den, 0) / self._dxe[1:-1]

        # Reconstruct the face states
        Qp = Q[...,1:-1] + dQ_lim * (self._dxp[1:-1] - v_edge[1:  ]*dt/2.)
        Qm = Q[...,1:-1] + dQ_lim * (self._dxm[1:-1] - v_edge[ :-1]*dt/2.)

        return np.where(v_edge[1:-1] > 0, Qp[...,:-1], Qm[...,1:])


class Weno3(object):
    '''Third-order upwind WENO reconstruction.

    Uses the geometrically consistent formulation of Mignone (2014).

    args:
        xe : Cell edge locations
        m  : Index of radial scaling giving volume, V = x^{m+1} / m+1.
    '''
    STENCIL = 2
    def __init__(self, xe, m, coeff=1.0):

        self._xe = xe

# Path: DiscEvolution/FV_interpolation.py
# def compute_centroids(xe, m):
#     '''First order upwind reconstruction
# 
#     args:
#         xe : Cell edge locations
#         m  : Index of radial scaling giving volume, V = x^{m+1} / m+1.
#     '''
#     return ((m + 1) * np.diff(xe**(m+2))) / ((m + 2) * np.diff(xe**(m+1)))
# 
# def construct_FV_edge_weights(xi, m, iL, iR, max_deriv=None, dtype='f8'):
#     '''Solves for the finite-volume interpolation weights.
# 
#     This code follows the methods outlined in Mignone (2014, JCoPh 270 784) to
#     compute the weights needed to reconstruct a function and its derivatives
#     to edges of computational cells. The polynomial is reconstructed to
#     reproduce the averages of the cell and its neighbours.
# 
#     Note that the polynomial computed near the domain edges will be lower 
#     order, which also reduces the number of derivatives available
# 
#     args:
#         xi : locations of cell edges
#         m  : Index of radial scaling giving volume, V = x^{m+1} / m+1.
#         iL : Number of cells to the left of the target cell to use in the 
#              interpolation.
#         iR : Number of cells to the right of the target cell to use in the 
#              interpolation.
#         
#         max_deriv : maximum derivative level to calculate. If not specified,
#                     return all meaningful derivatives.
#     
#         dtype : data type used for calculation, default = 'f8'
# 
#     returns:
#         wp, wm : The weights used for reconstructing the function and its 1st
#                  iL+iR derivatives to the left and right of the cell edges. 
#                  The shape is [len(xi)-1,max_deriv+1, 1+iL+iR]
#     '''
#     # Order of the polynomial
#     order = 1 + iL + iR
# 
#     if max_deriv is None:
#         max_deriv = order - 1
#     elif max_deriv > order - 1:
#         raise ValueError("Maximum derivative must be less than the order of the"
#                          " polynomial fitted")
#     
#     # Setup the beta matrix of Mignone:
#     beta =  _construct_volume_factors(xi, m, order, dtype)
# 
#     # The matrix of extrapolations to the RHS of cell
#     wp = _solve_FV_matrix_weights(xi[1:], iL, iR, beta, max_deriv, dtype)
#         
#     # The matrix of extrapolations to the LHS of cell
#     wm = _solve_FV_matrix_weights(xi[:-1], iR, iL, beta, max_deriv, dtype)
# 
#     return wp, wm


# Path: DiscEvolution/reconstruction.py
import numpy as np
    import matplotlib.pyplot as plt
from scipy.special import erf
from .FV_interpolation import compute_centroids, construct_FV_edge_weights

        dQ = Q[...,1:] - Q[...,:-1]
        QF, QB = dQ[...,1:], dQ[...,:-1]

        cF, cB = self._cF, self._cB

        num = QF*QB*(cF*QB + cB*QF)
        den = (QB*QB + (cF + cB - 2)*QB*QF + QF*QF) + 1e-300

        dQ_lim = np.where(QB*QF > 0, num/den, 0) / self._dxe[1:-1]

        # Reconstruct the face states
        Qp = Q[...,1:-1] + dQ_lim * (self._dxp[1:-1] - v_edge[1:  ]*dt/2.)
        Qm = Q[...,1:-1] + dQ_lim * (self._dxm[1:-1] - v_edge[ :-1]*dt/2.)

        return np.where(v_edge[1:-1] > 0, Qp[...,:-1], Qm[...,1:])


class Weno3(object):
    '''Third-order upwind WENO reconstruction.

    Uses the geometrically consistent formulation of Mignone (2014).

    args:
        xe : Cell edge locations
        m  : Index of radial scaling giving volume, V = x^{m+1} / m+1.
    '''
    STENCIL = 2
    def __init__(self, xe, m, coeff=1.0):

        self._xe = xe

        self._xc = xc = compute_centroids(xe, m)