How to depict small charges on a spherical object using vpython library?

I am working on a project related to charge distribution on the sphere and I decided to simulate the problem using vpython and Coulomb's law. I ran into an issue when I created a sphere because I am trying to evenly place out like 1000 points (charges) on the sphere and I can't seem to succeed, I have tried several ways but can't seem to make the points be on the sphere.

I defined an arbitrary value SOYDNR as a number to divide the diameter of the sphere into smaller segments. This would allow me to create a smaller rings of charges and fill out the surface of the spahre with charges. Then I make a list with 4 values that represent different parts of the radius to create the charge rings on the surface. Then I run into a problem and I am not sure how to deal with it. I have tried calculating the radius at those specific heights but yet I couldn't implement it. This is how it looks visually:![Sphere with charges on the surface].(https://i.sstatic.net/3N4x6.png) If anyone has any suggestions, I would be greatful, thanks!

 SOYDNR = 10 #can be changed
  SOYD = 2*radi/SOYDNR # strips of Y direction; initial height + SOYD until = 2*radi
  theta = 0
  dtheta1 = 2*pi/NCOS
  y_list = np.arange(height - radi + SOYD, height, SOYD).tolist()
  print(y_list)
  for i in y_list:
   while Nr<NCOS and theta<2*pi:
    position = radi*vector(cos(theta),i*height/radi,sin(theta))
    points_on_sphere = points_on_sphere + [sphere(pos=position, radius=radi/50, color=vector(1, 0, 0))]
    Nr = Nr + 1
    theta = theta + dtheta1
   Nr = 0
   theta = 0

Solution

I found a great way to do it, it creates a bunch of spheres in the area that is described by an if statement this is the code I am using for my simulation that creates the sphere with points on it.

def SOSE (radi, number_of_charges, height):
  Charged_Sphere = sphere(pos=vector(0,height,0), radius=radi, color=vector(3.5, 3.5, 3.5), opacity=(0.2))
  points_on_sphere = []
  NCOS = number_of_charges
  theta = 0
  dtheta = 2*pi/NCOS
  dr = radi/60
  direcVector = vector(0, height, 0)
  while theta<2*pi:
    posvec1 = radi*vector(1-radi*random(),1-radi*random()/radi,1-radi*random())
    posvec2 = radi*vector(1-radi*random(),-1+radi*random()/radi,1-radi*random())
    if mag(posvec1)<radi and mag(posvec1)>(radi-dr):
        posvec1 = posvec1+direcVector
        points_on_sphere=points_on_sphere+[sphere(pos=posvec1,radius=radi/60,color=vector(1, 0, 0))] 
        theta=theta + dtheta
    if mag(posvec2)<radi and mag(posvec2)>(radi-dr):
        posvec2 = posvec2+direcVector
        points_on_sphere=points_on_sphere+[sphere(pos=posvec2,radius=radi/60,color=vector(1, 0, 0))]
        theta=theta + dtheta

This code can be edited to add more points and I have two if statements because I want to change the height at which the sphere is present, and if I have just one statement I only see half of the sphere. :)