Generate random sampling on a sphere surface with probability function in latitude - python

Question

I'm trying to generate a random events in the sky given the right ascension (longitude) and declination (latitude). But i don't want to be uniform distributed. The events needs to be distributed according to a probability function that i know how to calculate in each point. That function depends only in the declination (latitude). The probability function (unormalized) can be viewed below: Probability function unormalized

I tried to generate the random sampling using two methods: The first one is generating a random declination (latitude) from the CDF, where i generate a random number in the interval [0,1] and, based on the CDF, generate a random declination. Then, for the longitude i just generate a random number in the interval [0,2*pi]: Sampling using the first method The problem of this method is that, normaly, when using an uniform distribution to generate the samples in the surface of a sphere, you need to account for the fact that just generating a random latitude and longitude, near the poles, you will have more points than near the equator, and i don't think that this is being taking in account in this method.

The second method is using the Metropolis Algorithm (MCMC). I generate a random point in the sphere, and then, i walk in some direction by 45°: Sampling using the second method The problem is that, the "size" of the walk (45°) appear to influence in the result, and, when i generate less points, it looks like even worse.

The first method is correct or there is a way to correct it? There is a better way to solve this problem?

Below is some of the code:

def ExposicaoRelativa(dec):
if (np.cos(dec) == 0):
    if ((np.cos(zenitalMax)-np.sin(latitudeDetector)*np.sin(dec)) > 0):
        alpha = 0.0
    elif (np.cos(zenitalMax)-np.sin(latitudeDetector)*np.sin(dec) < 0):
        alpha = np.pi
else:   
    sup = (np.cos(zenitalMax)-np.sin(latitudeDetector)*np.sin(dec))/(np.cos(latitudeDetector)*np.cos(dec))
    if (sup > 1):
        alpha = 0
    elif (sup < -1):
        alpha = np.pi
    else:
        alpha = np.arccos(sup)
return (np.sin(alpha)*np.cos(latitudeDetector)*np.cos(dec)+alpha*np.sin(latitudeDetector)*np.sin(dec))/np.pi

where 'dec' is the latitude, zenitalMax and latitudeDectector are variables given in another part of the code.