<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def watsonsV(Dir1, Dir2):
 """ calculates Watson's V statistic for two sets of directions """

    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]