import numpy as np
import scipy.integrate

def _bd_mass_integrand(rho, z, bd):
    rho = np.array([rho])
    z   = np.array([z])
    psi      = bd.eval_psi(rho, z, [0, 0])
    dpsi_rho = bd.eval_psi(rho, z, [1, 0])
    dpsi_z   = bd.eval_psi(rho, z, [0, 1])
    ret = rho * ((dpsi_rho**2) + (dpsi_z**2)) / (psi ** 2)
    return ret[0, 0]

def calc_mass(bd, tol = 1e-4, bound = 200):
    """
    Calculate the 
    """
    val, err = scipy.integrate.dblquad(_bd_mass_integrand, 0, bound,
                                       lambda x: 0, lambda x: bound,
                                       args = (bd,), epsrel = tol, epsabs = tol)
    return 2 * val, err
