docs: def plot_ell(fignum, pars, col, lower, plot):
      """
      function to calcualte/plot points on an ellipse about Pdec,Pdip with angle beta,gamma
      Parameters
      _________
      fignum : matplotlib figure number
      pars : list of [Pdec, Pinc, beta, Bdec, Binc, gamma, Gdec, Ginc ]
           where P is direction, Bdec,Binc are beta direction, and Gdec,Ginc are gamma direction
      col : color for ellipse
      lower : boolean, if True, lower hemisphere projection
      plot : boolean, if False, return the points, if False, make the plot
      """
      plt.figure(num=fignum)
      rad = old_div(np.pi, 180.)
      Pdec, Pinc, beta, Bdec, Binc, gamma, Gdec, Ginc = pars[0], pars[
          1], pars[2], pars[3], pars[4], pars[5], pars[6], pars[7]
      if beta > 90. or gamma > 90:
          beta = 180. - beta
          gamma = 180. - gamma
          Pdec = Pdec - 180.
          Pinc = -Pinc
      beta, gamma = beta * rad, gamma * rad  # convert to radians
      X_ell, Y_ell, X_up, Y_up, PTS = [], [], [], [], []
      nums = 201
      xnum = old_div(float(nums - 1.), 2.)
  # set up t matrix
      t = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
      X = pmag.dir2cart((Pdec, Pinc, 1.0))  # convert to cartesian coordintes
      if lower == 1 and X[2] < 0:
          for i in range(3):
              X[i] = -X[i]
  # set up rotation matrix t
      t[0][2] = X[0]
      t[1][2] = X[1]
      t[2][2] = X[2]
      X = pmag.dir2cart((Bdec, Binc, 1.0))
      if lower == 1 and X[2] < 0:
          for i in range(3):
              X[i] = -X[i]
      t[0][0] = X[0]
      t[1][0] = X[1]
      t[2][0] = X[2]
      X = pmag.dir2cart((Gdec, Ginc, 1.0))
      if lower == 1 and X[2] < 0:
          for i in range(3):
              X[i] = -X[i]
      t[0][1] = X[0]
      t[1][1] = X[1]
      t[2][1] = X[2]
  # set up v matrix
      v = [0, 0, 0]
      for i in range(nums):  # incremental point along ellipse
          psi = float(i) * np.pi / xnum
          v[0] = np.sin(beta) * np.cos(psi)
          v[1] = np.sin(gamma) * np.sin(psi)
          v[2] = np.sqrt(1. - v[0]**2 - v[1]**2)
          elli = [0, 0, 0]
  # calculate points on the ellipse
          for j in range(3):
              for k in range(3):
                  # cartesian coordinate j of ellipse
                  elli[j] = elli[j] + t[j][k] * v[k]
          pts = pmag.cart2dir(elli)
          PTS.append([pts[0], pts[1]])
          # put on an equal area projection
          R = old_div(np.sqrt(
              1. - abs(elli[2])), (np.sqrt(elli[0]**2 + elli[1]**2)))
          if elli[2] <= 0:
              #            for i in range(3): elli[i]=-elli[i]
              X_up.append(elli[1] * R)
              Y_up.append(elli[0] * R)
          else:
              X_ell.append(elli[1] * R)
              Y_ell.append(elli[0] * R)
      if plot == 1:
          col = col[0]+'.'
          if X_ell != []:
              plt.plot(X_ell, Y_ell, col, markersize=3)
          if X_up != []:
              plt.plot(X_up, Y_up, col, markersize=3)
      else:
          return PTS