code: def calculate_best_fit_vectors(L, E, V, n_planes):
      
  
      U, XV = E[:], []  # make a copy of E to prevent mutation
      for pole in L:
          XV.append(vclose(pole, V))  # get some points on the great circle
          for c in range(3):
              U[c] = U[c] + XV[-1][c]
  # iterate to find best agreement
      angle_tol = 1.
      while angle_tol > 0.1:
          angles = []
          for k in range(n_planes):
              for c in range(3):
                  U[c] = U[c] - XV[k][c]
              R = np.sqrt(U[0]**2 + U[1]**2 + U[2]**2)
              for c in range(3):
                  V[c] = old_div(U[c], R)
              XX = vclose(L[k], V)
              ang = XX[0] * XV[k][0] + XX[1] * XV[k][1] + XX[2] * XV[k][2]
              angles.append(np.arccos(ang) * 180. / np.pi)
              for c in range(3):
                  XV[k][c] = XX[c]
                  U[c] = U[c] + XX[c]
              amax = -1
              for ang in angles:
                  if ang > amax:
                      amax = ang
              angle_tol = amax
  
      return XV
docs: Calculates the best fit vectors for a set of plane interpretations used in fisher mean calculations
      @param: L - a list of the "EL, EM, EN" array of MM88 or the cartisian form of dec and inc of the plane interpretation
      @param: E - the sum of the cartisian coordinates of all the line fits to be used in the mean
      @param: V - inital direction to start iterating from to get plane best fits
      @returns: nested list of n_plane by 3 dimension where the 3 are the cartisian dimension of the best fit vector