  A new numerical code for computing stationary axisymmetric rapidly rotating
stars in general relativity is presented. The formulation is based on a fully
constrained-evolution scheme for 3+1 numerical relativity using the Dirac gauge
and maximal slicing. We use both the polytropic and MIT bag model equations of
state to demonstrate that the code can construct rapidly rotating neutron star
and strange star models. We compare numerical models obtained by our code and a
well-established code, which uses a different gauge condition, and show that
the two codes agree to high accuracy.
