Sampling angles on a sphere: From a general covariance matrix to concentration parameter

Question

I am trying to extract points on a sphere in python.

I have to localise event in the sky and produce a map using healpy. During the test session I used the von Mises-Fisher since I assumed that variance for $'\theta'$ is the same as the one for $'\phi'$. All worked well and I was able to obtain a concentration parameter $'k'$ by using $'k=1/\sigma^2'$.

My function to evaluate the probability in a pixel was

def Mises_Fisher(theta,phi,DS_theta,DS_phi,conc):
    meanvec=hp.ang2vec(DS_theta,DS_phi)
    meanvec=np.asarray(meanvec,dtype=np.float128)
    norm=np.sqrt(np.dot(meanvec,meanvec))
    meanvec=meanvec/norm
    
    var=hp.ang2vec(theta,phi)
    var=np.asarray(var,dtype=np.float128)
    norm=np.sqrt(np.dot(var,var))
    var=var/norm
    
    factor=np.dot(conc*var,meanvec)
    factor=np.float128(factor)
    #Normalization is futile, we will devide by the sum
    #fullnorm=conc/(2*np.pi*(np.exp(conc)-np.exp(-conc)))
    ret=np.float128(np.exp(factor))#/fullnorm
    #ret=factor
    return ret

and then for each pixel in the healpy map, I can assign a probability.

Now I have a covariance matrix that is no more isotropic, how can I generate points on a sphere in this case? I know the existence of the Kent distribution, but I don't have the required quantities, like $'\gamma'$,$'\beta'$ and $'k'$. I don't know if there is a way to obtain those given a covariance matrix.

To further complicate the problem, the matrix is in 11 dimension but I need only the projection to the (distance,theta,phi)-space to then produce the map.

For now I tried to extract from a multivariate Gaussian but since the Gaussian is on the tangent plane, for some variances the distribution produce angles that are outside the domain for $'\theta'$ and $'\phi'$.

Here is the code I'm using to extract the points

num_samples = 10**7
samples = np.random.multivariate_normal(perm_mean, perm_cov, num_samples)

phi = samples[:, 2]
theta = samples[:, 1]
print(np.min(theta),np.max(theta))
print(np.min(phi),np.max(phi))
print('starting mean values')
print('theta={}, phi={}'.format(perm_mean[1],perm_mean[2]))

The point is that for some matrices, the variance is such that angles can go outside the domain. I can't act on the matrices, those are given to me by some other codes and should be treated as fixed data.

Thank you all in advance!