def dia_vgp(*args):  # new function interface by J.Holmes, SIO, 6/1/2011
    """
    Converts directional data (declination, inclination, alpha95) at a given
    location (Site latitude, Site longitude) to pole position (pole longitude,
    pole latitude, dp, dm)

    Parameters
    ----------
    Takes input as (Dec, Inc, a95, Site latitude, Site longitude)
    Input can be as individual values (5 parameters)
    or
    as a list of lists: [[Dec, Inc, a95, lat, lon],[Dec, Inc, a95, lat, lon]]

    Returns
    ----------
    if input is individual values for one pole the return is:
    pole longitude, pole latitude, dp, dm

    if input is list of lists the return is:
    list of pole longitudes, list of pole latitude, list of dp, list of dm
    """
    # test whether arguments are one 2-D list or 5 floats
    if len(args) == 1:  # args comes in as a tuple of multi-dim lists.
        largs = list(args).pop()  # scrap the tuple.
        # reorganize the lists so that we get columns of data in each var.
        (decs, dips, a95s, slats, slongs) = list(zip(*largs))
    else:
        # When args > 1, we are receiving five floats. This usually happens when the invoking script is
        # executed in interactive mode.
        (decs, dips, a95s, slats, slongs) = (args)

    # We send all incoming data to numpy in an array form. Even if it means a
    # 1x1 matrix. That's OKAY. Really.
    (dec, dip, a95, slat, slong) = (np.array(decs), np.array(dips), np.array(a95s),
                                    np.array(slats), np.array(slongs))  # package columns into arrays
    rad = old_div(np.pi, 180.)  # convert to radians
    dec, dip, a95, slat, slong = dec * rad, dip * \
        rad, a95 * rad, slat * rad, slong * rad
    p = np.arctan2(2.0, np.tan(dip))
    plat = np.arcsin(np.sin(slat) * np.cos(p) +
                     np.cos(slat) * np.sin(p) * np.cos(dec))
    beta = old_div((np.sin(p) * np.sin(dec)), np.cos(plat))

    # -------------------------------------------------------------------------
    # The deal with "boolmask":
    # We needed a quick way to assign matrix values based on a logic decision, in this case setting boundaries
    # on out-of-bounds conditions. Creating a matrix of boolean values the size of the original matrix and using
    # it to "mask" the assignment solves this problem nicely. The downside to this is that Numpy complains if you
    # attempt to mask a non-matrix, so we have to check for array type and do a normal assignment if the type is
    # scalar. These checks are made before calculating for the rest of the function.
    # -------------------------------------------------------------------------

    boolmask = beta > 1.  # create a mask of boolean values
    if isinstance(beta, np.ndarray):
        beta[boolmask] = 1.  # assigns 1 only to elements that mask TRUE.
    # Numpy gets upset if you try our masking trick with a scalar or a 0-D
    # matrix.
    else:
        if boolmask:
            beta = 1.
    boolmask = beta < -1.
    if isinstance(beta, np.ndarray):
        beta[boolmask] = -1.  # assigns -1 only to elements that mask TRUE.
    else:
        if boolmask:
            beta = -1.

    beta = np.arcsin(beta)
    plong = slong + np.pi - beta
    if (np.cos(p) > np.sin(slat) * np.sin(plat)).any():
        boolmask = (np.cos(p) > (np.sin(slat) * np.sin(plat)))
        if isinstance(plong, np.ndarray):
            plong[boolmask] = (slong + beta)[boolmask]
        else:
            if boolmask:
                plong = slong + beta

    boolmask = (plong < 0)
    if isinstance(plong, np.ndarray):
        plong[boolmask] = plong[boolmask] + 2 * np.pi
    else:
        if boolmask:
            plong = plong + 2 * np.pi

    boolmask = (plong > 2 * np.pi)
    if isinstance(plong, np.ndarray):
        plong[boolmask] = plong[boolmask] - 2 * np.pi
    else:
        if boolmask:
            plong = plong - 2 * np.pi

    dm = np.rad2deg(a95 * (old_div(np.sin(p), np.cos(dip))))
    dp = np.rad2deg(a95 * (old_div((1 + 3 * (np.cos(p)**2)), 2)))
    plat = np.rad2deg(plat)
    plong = np.rad2deg(plong)
    return plong.tolist(), plat.tolist(), dp.tolist(), dm.tolist()
Converts directional data (declination, inclination, alpha95) at a given
    location (Site latitude, Site longitude) to pole position (pole longitude,
    pole latitude, dp, dm)

    Parameters
    ----------
    Takes input as (Dec, Inc, a95, Site latitude, Site longitude)
    Input can be as individual values (5 parameters)
    or
    as a list of lists: [[Dec, Inc, a95, lat, lon],[Dec, Inc, a95, lat, lon]]

    Returns
    ----------
    if input is individual values for one pole the return is:
    pole longitude, pole latitude, dp, dm

    if input is list of lists the return is:
    list of pole longitudes, list of pole latitude, list of dp, list of dm