def b_vdm(B, lat):
    """
    Converts a magnetic field value (input in units of tesla) to a virtual
    dipole moment (VDM) or a virtual axial dipole moment (VADM); output
    in units of Am^2)

    Parameters
    ----------
    B: local magnetic field strength in tesla
    lat: latitude of site in degrees

    Returns
    ----------
    V(A)DM in units of Am^2

    Examples
    --------
    >>> pmag.b_vdm(33e-6,22)*1e-21

    71.58815974511788
    """
    # changed radius of the earth from 3.367e6 3/12/2010
    fact = ((6.371e6)**3) * 1e7
    colat = np.radians(90. - lat)
    return fact * B / (np.sqrt(1 + 3 * (np.cos(colat)**2)))