starfield simulation, probably should use direct3d/opengl for better realism and faster frame rate. stars will not just be points and if they come close enough they can take up more than one pixel.
Take into account the 7 types of stars (http://www.enchantedlearning.com/subjects/astronomy/stars/startypes.shtml) for color, size, brightness and frequencies of the respective types?
or base the probabilities on this: http://mirror-uk-rb1.gallery.hd.org/_exhibits/natural-science/_more2000/_more10/stars-in-night-sky-over-London-England-tweaked-2-DHD.jpg
but they would be different in different parts of a galaxy or outside a galaxy.

if not direct3d, use directdraw or direct3d with a 2d surface, and the algorithm for calculating positions can be extremely simple:
each star has an absolute location in 3d space with x=0, y=0 being directly in front of you (and z being distance)
the monitor is simply a window to the stars in between them and your eye.
therefore, to calculate the the x and y positions, calculate the point at which the line between the star and your eye intersects the screen.
that will be the star's x and y scaled according to its z. for x and y just divide the star's x and y by its z and multiply by the supposed distance between your eye and the monitor (like 1.5 to 2 feet, but use the same units as your star distances) and multiply by the monitor's resolution and divide by the size of the monitor. your distance * resolution / size of monitor can be calculated once and kept in a variable (say "K"). so all you have to do is calculate x/z*K, y/z*K, or even better, scale your x's and y's or your z's to cancel out K when they're generated.
since the star's x and y are centered at 0 you'll have to add width/2 and height/2 to your resultant x and y values.
the star's x's and y's won't change, but each iteration its z will change. as you move forward in space they'll relatively move closer, decrementing their z. but if you make the motion through space and the screen updating asynchronous, don't go through and change each star's z for each motion update. just decrement a universal z offset and add that to the star's z when you draw to screen.
when a star moves out of view (when your calculation yields an x or a y out of range), delete it from your array and replace it with a newly randomized star.
to have stars bigger than one pixel, use a circle and blit it to a size proportional to 1/(x^2+y^+z^2). for results smaller than 1 pixel, dim the pixels proportionately.
so it looks like this, simplified:

K = 1.5 feet * monitor's DPI # adjust feet and dpi to the units used for stars x/y/z
uz = 0
wo2 = width/2
ho2 = height/2
def randomstar():
  x = random positive or negative value * K
  y = random positive or negative value * K
  z = random value + uz
  return x, y, z
stars = []
for star in xrange(numberofstars):
  stars.append(randomstar())
while 1:
  screen.clear()
  uz += distance you travel forward
  for starindex in xrange(numberofstars):
    starx, stary, starz = stars[starindex]
    x = starx/(starz-uz) # check to not divide by zero, although that could be exceedinly unlikely
    y = stary/(starz-uz)
    if 0 <= x < width and 0 <= y < height: # modify this a little becasue part of a star might show while the center of it is out of the picture
      draw star at x, y
    else:
      star[starindex] = randomstar()
  screen.blit()