  We introduce a constrained Monte Carlo method which allows us to traverse the
phase space of a classical spin system while fixing the magnetization
direction. Subsequently we show the method's capability to model the
temperature dependence of magnetic anisotropy, and for bulk uniaxial and cubic
anisotropies we recover the low-temperature Callen-Callen power laws in M. We
also calculate the temperature scaling of the 2-ion anisotropy in L10 FePt, and
recover the experimentally observed M^2.1 scaling. The method is newly applied
to evaluate the temperature dependent effective anisotropy in the presence of
the N'eel surface anisotropy in thin films with different easy axis
configurations. In systems having different surface and bulk easy axes, we show
the capability to model the temperature-induced reorientation transition. The
intrinsic surface anisotropy is found to follow a linear temperature behavior
in a large range of temperatures.
