code: def watsonsV(Dir1, Dir2):
      
      counter, NumSims = 0, 500
  #
  # first calculate the fisher means and cartesian coordinates of each set of Directions
  #
      pars_1 = fisher_mean(Dir1)
      pars_2 = fisher_mean(Dir2)
  #
  # get V statistic for these
  #
      V = vfunc(pars_1, pars_2)
  #
  # do monte carlo simulation of datasets with same kappas, but common mean
  #
      Vp = []  # set of Vs from simulations
      print("Doing ", NumSims, " simulations")
      for k in range(NumSims):
          counter += 1
          if counter == 50:
              print(k + 1)
              counter = 0
          Dirp = []
  # get a set of N1 fisher distributed vectors with k1, calculate fisher stats
          for i in range(pars_1["n"]):
              Dirp.append(fshdev(pars_1["k"]))
          pars_p1 = fisher_mean(Dirp)
  # get a set of N2 fisher distributed vectors with k2, calculate fisher stats
          Dirp = []
          for i in range(pars_2["n"]):
              Dirp.append(fshdev(pars_2["k"]))
          pars_p2 = fisher_mean(Dirp)
  # get the V for these
          Vk = vfunc(pars_p1, pars_p2)
          Vp.append(Vk)
  #
  # sort the Vs, get Vcrit (95th one)
  #
      Vp.sort()
      k = int(.95 * NumSims)
      return V, Vp[k]
docs: calculates Watson's V statistic for two sets of directions