  We have computed models of rotating relativistic stars with a toroidal
magnetic field and investigated the combined effects of magnetic field and
rotation on the apparent shape (i.e. the surface deformation), which could be
relevant for the electromagnetic emission, and on the internal matter
distribution (i.e. the quadrupole distortion), which could be relevant for the
emission of gravitational waves. Using a sample of eight different cold nuclear
physics equations of state, we have computed models of maximum field strength,
as well as the distortion coefficients for the surface and the quadrupolar
deformations. Surprisingly, we find that non-rotating models admit arbitrary
levels of magnetization, accompanied by a growth of size and quadrupole
distortion to which we could not find a limit. Rotating models, on the other
hand, are subject to a mass-shedding limit at frequencies well below the
corresponding ones for unmagnetized stars. Overall, the space of solutions can
be split into three distinct classes for which the surface deformation and the
quadrupole distortion are either prolate and prolate, oblate and prolate, or
oblate and oblate, respectively. We also derive a simple formula expressing the
relativistic distortion coefficients, which allows one to compute the surface
deformation and the quadrupole distortion up to significant levels of rotation
and magnetization, essentially covering all known magnetars. Such a formula
replaces Newtonian equivalent expressions that overestimate the magnetic
quadrupole distortion by about a factor of 6 and are inadequate for strongly
relativistic objects like neutron stars.
