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))

