code: def sky(lon=None,lat=None,size=1):
      
      if lon is None:
          umin,umax = 0,1
      else:
          lon = np.asarray(lon)
          lon = np.radians(lon + 360.*(lon<0))
          if   lon.size==1: umin=umax=lon/(2*np.pi)
          elif lon.size==2: umin,umax=lon/(2*np.pi)
          else: raise Exception('...')
          
      if lat is None:
          vmin,vmax = -1,1
      else:
          lat = np.asarray(lat)
          lat = np.radians(90 - lat)
          if   lat.size==1: vmin=vmax=np.cos(lat)
          elif lat.size==2: vmin,vmax=np.cos(lat)
          else: raise Exception('...')
  
      phi = 2*np.pi*np.random.uniform(umin,umax,size=size)
      theta = np.arcsin(np.random.uniform(vmin,vmax,size=size))
      return np.degrees(phi),np.degrees(theta)
docs: Outputs uniform points on sphere from:
      [0 < lon < 360] & [-90 < lat < 90]