Distribute points on a half sphere randomly without clustering in the center

This is a follow up question to this one: pyephem problem with alt/az to ra/dec and back

I want to distribute stars randomly above an observer. I was able to to do this in the last question.

The issue I now get is that the stars seems clustered towards the center.

How can I get rid of this?

import matplotlib.pyplot as plt
import numpy as np
import  math
import ephem
from datetime import datetime
import random

home = ephem.Observer()

home.lon = '-70.4'  # +E
home.lat = '-24.62'  # +N
home.elevation = 0  # meters
home.date = datetime.utcnow()

theta_plot = []
r_plot=[]

stars_pre_selected = []

for n in range(4000):
    ra, dec = home.radec_of(az=random.uniform(0, 2 * math.pi),
                            alt=random.uniform(0, math.pi))

    body = ephem.FixedBody()
    body._epoch = ephem.J2000
    body._ra = ephem.degrees(ra)
    body._dec = ephem.degrees(dec)

    body.compute(home)


    if math.degrees(body.alt) > 0:
        stars_pre_selected.append(body)


    else:
        print(math.degrees(body.alt))


for body in stars_pre_selected:
    body.compute(home)
    theta_plot.append(body.az)
    r_plot.append(math.cos(body.alt))
    print(body.alt,body.az)

ax = plt.subplot(111, polar=True)
#
ax.set_theta_direction(-1)
ax.set_theta_offset(np.pi / 2)
ax.grid(True)
ax.scatter(theta_plot, r_plot, c='blue', s=2)
ax.set_rmax(1.0)

plt.pause(0.1)

plt.show()

UPDATE:

    ra, dec = home.radec_of(az=random.uniform(0, 2 * math.pi),
    alt=math.asin(2*random.random()-1))
    

this will work too, but then the code hast to be changed to math.sin instead of math.cos :

    ra, dec = home.radec_of(az=random.uniform(0, 2 * math.pi),
    alt=math.acos(2*random.random()-1))

ra, dec = home.radec_of(az=random.uniform(0, 2 * math.pi),
alt=math.asin(2*random.random()-1))