2D, 3D Ising model; XY model C code. Want to simulate a ferromagnet?

The Ising model attempts to describe the low-temperature behaviour of ferromagnetic materials, which show a phase transition (means: from simple metal to magnet) at a certain temperature near 0°K. These ferromagnetic materials essentially become magnets just bringing their temperature below a critical value.
The reaching of this critical temperature manifests itself with sharp peaks or drops (jumps) among the values of various "observables", or macroscopic thermodynamical quantities, which are directly measurable given a computer simulation of the Ising model; this computer simulation is based on the Monte Carlo technique, and thus relies on random number generation in order to simulate the dynamics of spin flipping in a sample ferromagnet.
Such code decides where and when to flip spins in an initially random spin matrix (can be a 2D or 3D array) following two rules in logical OR:

• Energy of the whole system is lower if this spin is flipped
• A random probability associated with this spin flip is < or = than the negative exponential function of the change in energy and the inverse of the temperature

In my code, I used the Metropolis algorithm for the acceptance of the spin flip; basically it assumes that a spin gets always flipped if the resulting energy is lower than the starting one.

Ising 2D and 3D models are not that different: 3D model is just an extension of the 2D model code. The XY model is  based on the 2D Ising code too, but it relies on a different Hamiltonian for the system.

My code does the following: it is divided into 2 parts, Simulation and Data Analysis.

• Loops for N_beta values of the temperature. At every loop it first initializes the spin system to a fixed configuration, then it runs a MonteCarlo simulation of the system at that temperature. It writes the values of the total Magnetization per site and of the total energy per site into a measurement file.
• A data analysis program takes the measurement file and computes thermodynamic observable values, writing them to a final data file.

DIRECTIONS for the CODE:
Repeat for each Simulation (run first) and DataAnalysis:

• Set the SIDE of the spin array (I use around 14)
• Set the number of sweeps (MonteCarlo iterations, about 40000)
• Set beta_minimum, beta_maximum, N_beta (I use .25, .9,40)
• Compile and Run 

DOWNLOAD:
 https://drive.google.com/file/d/0B36YNo_GwIhxNnZPSVlfTXl2OXc/view?usp=sharing