Simple Implementation of 2D Ising model. Initialize system with edge length and temperature, sweep until equilibrated. Plot to taste.

Ising.cpp

#include "Ising.h"
#include <cmath>
#ifndef _CSTDLIB
#define _CSTDLIB
#include <cstdlib>
#endif

Ising::Ising(unsigned int L_, double temp_){
    L = L_;
    Lsq = L*L;
    temp = temp_;
    p2 = exp(-4.0/temp);
    p4 = exp(-8.0/temp);
    grid = new int[L*L];
    unsigned int i;
    for (i = 0; i < Lsq; i++) grid[i] = 1;
}

Ising::~Ising(){
    delete[] grid;
};

// Mersenne Twister for better random numbers
double random(){
#define N     624
#define M     397
#define HI    0x80000000
#define LO    0x7fffffff
#define MAX   0xFFFFFFFF
	static const uint32_t seed = 5489UL;
	static const uint32_t A[2] = { 0, 0x9908b0df };
	static uint32_t y[N];
	static int i = N + 1;

	if (i > N) {
		/* Initialize with pseudo random number */
		y[0] = seed;

		for (i = 1; i < N; ++i) {
			y[i] = (1812433253UL * (y[i - 1] ^ (y[i - 1] >> 30)) + i);
		}
	}

	if (i == N) {
		/* calculate new state vector */
		uint32_t h;

		for (i = 0; i < N - M; ++i) {
			h = (y[i] & HI) | (y[i + 1] & LO);
			y[i] = y[i + M] ^ (h >> 1) ^ A[h & 1];
		}
		for (; i < N - 1; ++i) {
			h = (y[i] & HI) | (y[i + 1] & LO);
			y[i] = y[i + (M - N)] ^ (h >> 1) ^ A[h & 1];
		}

		h = (y[N - 1] & HI) | (y[0] & LO);
		y[N - 1] = y[M - 1] ^ (h >> 1) ^ A[h & 1];
		i = 0;
	}

	uint32_t e = y[i++];
	/* tempering */
	e ^= (e >> 11);
	e ^= (e << 7) & 0x9d2c5680;
	e ^= (e << 15) & 0xefc60000;
	e ^= (e >> 18);

	return (double)e / MAX;
#undef N
#undef M
#undef HI
#undef LO
}

void Ising::sweep(){
    unsigned int i, xy, x, y;
    int deltaE;
    for(i = 0; i < Lsq; i++){
        xy = Lsq*random();
        y = xy%L;
        x = xy/L;
        deltaE = grid[xy]*(grid[L*((x-1)%L)+y]+grid[L*((x+1)%L)+y]+grid[L*x+(y-1)%L]+grid[L*x+(y+1)%L]); // calculates half of deltaE
        switch(deltaE){
        case 2:
            if(random() < p2) grid[xy] *= -1;
            break;
        case 4:
            if(random() < p4) grid[xy] *= -1;
            break;
        default: grid[xy] *= -1;
        }
    }
}

Ising.h

#ifndef ISING_H
#define ISING_H

class Ising{
public:
    Ising(unsigned int, double);
    ~Ising();
    void sweep();

    int* grid;


private:
    double temp, p2, p4;
    unsigned int L;
    unsigned int Lsq;
};

#endif //ISING_H