  The data from all observations of RBS 1223 (1RXS J130848.6+212708) conducted
by XMM-Newton EPIC pn with the same instrumental setup in 2003-2007 were
combined to form spin-phase resolved spectra. A number of complex models of
neutron stars with strongly magnetized ($B_{pole} 10^{13}-10^{14}$ G) surface,
with temperature and magnetic field distributions around magnetic poles, and
partially ionized hydrogen thin atmosphere above it have been implemented into
the X-ray spectral fitting package XSPEC for simultaneous fitting of
phase-resolved spectra. A Markov-Chain-Monte-Carlo (MCMC) approach is also
applied to verify results of fitting and estimating in multi parameter models.
The spectra in different rotational phase intervals and light curves in
different energy bands with high S/N ratio show a high complexity. The spectra
can be parameterized with a Gaussian absorption line superimposed on a
blackbody spectrum, while the light curves with double-humped shape show strong
dependence of pulsed fraction upon the energy band (13%-42%), which indicates
that radiation emerges from at least two emitting areas. A model with condensed
iron surface and partially ionized hydrogen thin atmosphere above it allows us
to fit simultaneously the observed general spectral shape and the broad
absorption feature observed at 0.3 keV in different spin phases of RBS 1223. It
allowed to constrain some physical properties of X-ray emitting areas, i.e. the
temperatures ($T_{p1} ~ 105$ eV, $T_{p2} ~ 99$ eV), magnetic field strengths
($B_{p1}= B_{p2} ~ 8.6x10^{13}$G) at the poles, and their distributions
parameters ($a_{1} 0.61, a_{2} 0.29$, indicating an absence of strong toroidal
magnetic field component). In addition, it puts some constraints on the
geometry of the emerging X-ray emission and gravitational redshift
($z=0.16^{+0.03}_{-0.01}$) of RBS 1223.
