#ifndef _INCLUDED_L3_ISING_H
#define _INCLUDED_L3_ISING_H

// *Ising Lattice Of Spins +0.5/-0.5* Copyright (C) Krzysztof Bosak, 1999-10-02...1999-10-04
// All rights reserved.
// kbosak@box43.pl
// http://www.kbosak.prv.pl

#include <stdlib.h>
#include "l3_rand.h"
#include "l3_spin.h"
#include "l3_array.h"

class IsingLattice
{
private:
	int _size;
	basearray<int> _xindexraw;
	basearray<int> _yindexraw;
	basearray<int> _exp_cache;
	SpinEnsemble _spins;
	int * _xindex;
	int * _yindex;
	double _temperature;

public:
	inline IsingLattice(int n, double t);
	inline explicit IsingLattice(const IsingLattice &il);
	inline IsingLattice& operator=(const IsingLattice &il);// EXCEPTION NEUTRAL
	inline void Init(int m);
	inline void InitSpinsUp(int m);
	inline void SetTemperature(double t);
	inline void MCS();
	inline void MCStotal();
	inline void Show() const;
	inline char Spin(int x, int y) const;
	inline int Width() const;
	inline int SpinsTotal() const;
	inline int SpinsUp() const;
	inline int SpinsDown() const;
	inline double Magnetisation() const;
};
void IsingLattice::SetTemperature(double t)
{// Sets temperature and updates caching of probability
	assert(t>0);
	_temperature=t;
	_exp_cache[2]=static_cast<int>(RAND.irandommax());
	_exp_cache[3]=static_cast<int>(RAND.irandommax()*exp(-1/_temperature));
	_exp_cache[4]=static_cast<int>(RAND.irandommax()*exp(-2/_temperature));
	_exp_cache[1]=_exp_cache[3];
	_exp_cache[0]=_exp_cache[4];
}
char IsingLattice::Spin(int x, int y) const
{// Returns true if spin at (x,y) is spin-up, periodic boundaries included in arrays...
	assert(x>=-1);
	assert(x<=_size);
	assert(y>=-1);
	assert(y<=_size);
	return _spins[_xindex[x]+_yindex[y]];
}
IsingLattice::IsingLattice(int n, double t)
	: _size(n),
	_xindexraw(n*3),
	_yindexraw(n*3),
	_exp_cache(5),
	_spins(n*n)
{// Ising Lattice constructor of nxn dimension
	_xindex=&_xindexraw[n];
	_yindex=&_yindexraw[n];
	assert(n>=5);// due to MCS evaluation algorithm
	const int sizex2=_size<<1;
	int i;
	int temp=0;
	for(i=0; i<_size; i++, temp+=_size)
	{
		_xindexraw[i]=i;
		_xindexraw[i+_size]=i;
		_xindexraw[i+sizex2]=i;
		_yindexraw[i]=temp;
		_yindexraw[i+_size]=temp;
		_yindexraw[i+sizex2]=temp;
	}
	SetTemperature(t);
}
inline IsingLattice::IsingLattice(const IsingLattice &il)
	: _size(il._size),
	_xindexraw(il._xindexraw),
	_yindexraw(il._yindexraw),
	_exp_cache(il._exp_cache),
	_spins(il._spins)
{
	_xindex=&_xindexraw[il._size];
	_yindex=&_yindexraw[il._size];
}
inline IsingLattice& IsingLattice::operator=(const IsingLattice &il)// EXCEPTION NEUTRAL
{
	if(&il==this)
	{
		return *this;
	}
	_size=il._size;
	_xindexraw=il._xindexraw;
	_xindex=&_xindexraw[il._size];
	_yindexraw=il._yindexraw;
	_yindex=&_yindexraw[il._size];
	_exp_cache=il._exp_cache;
	_spins=il._spins;
	return *this;
}
void IsingLattice::Init(int m)
{// Set spins: spins-up - spins-down == m
	const int size2=_spins.SpinsTotal()>>1;
	assert(abs(m)<=size2);
	_spins.Fill(m+size2);
}
void IsingLattice::InitSpinsUp(int spup)
{// Set spup spins to spin-up
	assert(spup>=0);
	assert(spup<=_spins.SpinsTotal());
	_spins.Fill(spup);
}
void IsingLattice::MCS()
{// Monte-Carlo Step using Metropolis Algorithm
	const int position=RAND.irandom(_spins.SpinsTotal());
	const int y=position/_size;
	const int x=position-y*_size;
	int sum=Spin(x-1,y)+Spin(x+1,y)+Spin(x,y-1)+Spin(x,y+1);
	assert(sum>=0 && sum<=4);
	char * const ptr=&_spins[position];
	if((*ptr)==1)
	{
		if(sum<=2 || RAND.irandom()<_exp_cache[sum])
			(*ptr)=0;
	}
	else
	{
		if(sum>=2 || RAND.irandom()<_exp_cache[sum])
			(*ptr)=1;
	}
}
void IsingLattice::MCStotal()
{// Monte-Carlo Step using Metropolis Algorithm for all spins
	int x;
	char *ptr=&_spins[0];
	for(x=0; x<_size; x++, ptr++)
	{// upper row
		const int sum=Spin(x-1,0)
			+Spin(x+1,0)
			+Spin(x,-1)
			+ptr[_size];
		assert(sum>=0 && sum<=4);
		if(*ptr==1)
		{
			if(sum<=2 || RAND.irandom()<_exp_cache[sum])
			{
				*ptr=0;
			}
		}
		else
		{
			if(sum>=2 || RAND.irandom()<_exp_cache[sum])
				(*ptr)=1;
		}
	}
	const int wm1=_size-1;
	ptr=&_spins[_yindex[wm1]];
	for(x=0; x<_size; x++, ptr++)
	{// lower row
		const int sum=Spin(x-1,wm1)
			+Spin(x+1,wm1)
			+ptr[-_size]
			+Spin(x,_size);
		assert(sum>=0 && sum<=4);
		if(*ptr==1)
		{
			if(sum<=2 || RAND.irandom()<_exp_cache[sum])
			{
				*ptr=0;
			}
		}
		else
		{
			if(sum>=2 || RAND.irandom()<_exp_cache[sum])
			{
				*ptr=1;
			}
		}
	}
	int y;
	ptr=&_spins[_size];
	for(y=1; y<wm1; y++, ptr+=_size)
	{// left column without edges
		const int sum=Spin(-1,y)+
			ptr[1]+
			ptr[-_size]+
			ptr[_size];
		assert(sum>=0 && sum<=4);
		if(*ptr==1)
		{
			if(sum<=2 || RAND.irandom()<_exp_cache[sum])
			{
				(*ptr)=0;
			}
		}
		else
		{
			if(sum>=2 || RAND.irandom()<_exp_cache[sum])
			{
				*ptr=1;
			}
		}
	}
	ptr=&_spins[_size+wm1];
	for(y=1; y<wm1; y++, ptr+=_size)
	{// right column without edges
		const int sum=ptr[-1]
			+Spin(wm1,y)
			+ptr[-_size]
			+ptr[_size];
		assert(sum>=0 && sum<=4);
		if(*ptr==1)
		{
			if(sum<=2 || RAND.irandom()<_exp_cache[sum])
			{
				*ptr=0;
			}
		}
		else
		{
			if(sum>=2 || RAND.irandom()<_exp_cache[sum])
			{
				*ptr=1;
			}
		}
	}
	char *ptry=&_spins[_size];
	for(y=1; y<wm1; y++, ptry+=_size)
	{// the rest - interior of the lattice
		for(x=1; x<wm1; x++)
		{
			char * const ptrinter=&ptry[x];
			const int sum=ptrinter[-1]+ptrinter[1]+ptrinter[-_size]+ptrinter[_size];
			assert(sum>=0 && sum<=4);
			if(*ptrinter==1)
			{
				if(sum<=2 || RAND.irandom()<_exp_cache[sum])
				{
					*ptrinter=0;
				}
			}
			else
			{
				if(sum>=2 || RAND.irandom()<_exp_cache[sum])
				{
					*ptrinter=1;
				}
			}
		}
	}
}
void IsingLattice::Show() const
{// Shows text output of lattice
	for(int x=0; x<_size; x++)
	{
		for(int y=0; y<_size; y++)
		{
			if(Spin(x,y))
			{
				cout<<'#';
			}
			else
			{
				cout<<' ';
			}
		}
		cout<<endl;
	}
}
int IsingLattice::Width() const
{// Returns lattice width
	return _size;
}
int IsingLattice::SpinsTotal() const
{// Returns number of all spins in lattice
	return _spins.SpinsTotal();
}
int IsingLattice::SpinsUp() const
{// Actual number of spins-up
	return _spins.SpinsUp();
}
int IsingLattice::SpinsDown() const
{// Actual number of spins-down
	return _spins.SpinsDown();
}
double IsingLattice::Magnetisation() const
{// Returns abs value of ( spins-up - spins-down ) / all_spins
	return _spins.Magnetisation();
}

#endif //_INCLUDED_L3_ISING_H