The python package for simulating Gaia data (using the formula from the Gaia science performance web pages) can be obtained from:

https://pypi.python.org/pypi/PyGaia

This code allows you to transform your phase space data (x,y,z,vx,vy,vz) into astrometry and radial velocities. In addition you can calculate the predicted Gaia errors as a function of magnitude and colour of your simulated stars.

For question/comments/suggestions/corrections for this code contact:

Anthony Brown <brown _at_ strw.leidenuniv.nl>

Assuming you've installed everything successfully, here's an example of simulating Gaia catalogue data.

from numpy.random import randn
from pygaia.astrometry.vectorastrometry import phaseSpaceToAstrometry
from pygaia.errors.astrometric import parallaxErrorSkyAvg
from pygaia.photometry.utils import vminiFromSpt
from pygaia.photometry.transformations import gminvFromVmini

# Convert your simulated phase space coordinates to Gaia observables.
# The coordinates are assumed to be barycentric (heliocentric). So make sure
# you agree on what position and velocity for the Sun to use.

l, b, parallax, mu_l, mu_b, vrad = phaseSpaceToAstrometry(x, y, z, vx, vy, vz)

# Obtain the errors on parallax for your set of stars.
# You have for each star V and (V-I) (possibly including extinction and reddening).

gmag = vmag + gminvFromVmini(vmini)
sigma_parallax = parallaxErrorSkyAvg(gmag, vmini)

# Simulate the observed parallaxes

parallax_observed = parallax + randn(parallax.size)*sigma_parallax

The code above calculates sky averaged parallaxes. If you want to account for the variations of the errors across the sky here's an examples of how to do that.

from pygaia.astrometry.coordinates import Transformations
from pygaia.astrometry.coordinates import CoordinateTransformation
from pygaia.errors.astrometric import parallaxError

# First obtain ecliptic sky-coordinates for your stars
coordinate_transformation = CoordinateTransformation(Transformation.ICRS2ECL)
lambda, beta = coordinate_transformation.transformSkyCoordinates(ra, dec)

# Calculate parallax errors as a function also of ecliptic latitude
sigma_parallax = parallaxError(gmag, vmini, beta)

For errors on the sky positions, proper motions, and radial velocities you can write code along the same lines.

Note that in it's current state the code happily returns errors for stars beyond Gaia's magnitude range G=20 (Grvs=17 for radial velocities). It's you responsibility to make sure stars fainter than these limits are not included in the simulated data.