docs: def bingham_pdf(fit):
      """
      From the *Encyclopedia of Paleomagnetism*
  
      From Onstott, 1980:
      Vector resultant: R is analogous to eigenvectors
      of T.
      Eigenvalues are analogous to |R|/N.
      """
      # Uses eigenvectors of the covariance matrix
      e = fit.hyperbolic_axes #singular_values
      #e = sampling_covariance(fit) # not sure
      e = e[2]**2/e
  
      kappa = (e-e[2])[:-1]
      kappa /= kappa[-1]
      F = N.sqrt(N.pi)*confluent_hypergeometric_function(*kappa)
  
  
      ax = fit.axes
  
      Z = 1/e
      M = ax
      F = 1/hyp1f1(*1/Z)
  
      def pdf(coords):
          lon,lat = coords
  
          I = lat
          D = lon# + N.pi/2
          #D,I = _rotate(N.degrees(D),N.degrees(I),90)
  
          # Bingham is given in spherical coordinates of inclination
          # and declination in radians
  
          # From USGS bingham statistics reference
  
          xhat = N.array(sph2cart(lon,lat)).T
  
          #return F*expm(dot(xhat.T, M, N.diag(Z), M.T, xhat))
  
          return 1/(F*N.exp(dot(xhat.T, M, N.diag(Z), M.T, xhat)))
  
      return pdf