  Thermal surface radiation has been detected with X-ray instruments for
several neutron stars with high spectral, spatial and timing resolution. These
observations allow for direct study of fundamental properties of the source,
but require model atmospheres and spectra for careful interpretation. We
describe a robust and extensible implementation of complete linearization for
computing the spectra of isolated cooling neutron stars for a broad range of
temperature and magnetic field. Self-consistent spectra are derived for
arbitrary magnetic field geometries at $B \leq 10^{14}$ G.
