docs: 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()