  We construct two-fluid equilibrium configurations for neutron stars with
magnetic fields, using a self-consistent and nonlinear numerical approach. The
two-fluid approach - likely to be valid for large regions of all but the
youngest NSs - provides us with a straightforward way to introduce
stratification and allows for more realistic models than the ubiquitous
barotropic assumption. In all our models the neutrons are modelled as a
superfluid, whilst for the protons we consider two cases: one where they are a
normal fluid and another where they form a type-II superconductor. We consider
a variety of field configurations in the normal-proton case and purely toroidal
fields in the superconducting case. We find that stratification allows for a
stronger toroidal component in mixed-field configurations, though the poloidal
component remains the largest in all our models. We provide quantitative
results for magnetic ellipticities of NSs, both in the normal- and
superconducting-proton cases.
