#ifndef _INCLUDED_L3_SPIN_H
#define _INCLUDED_L3_SPIN_H

// *Spin Up-Down Ensemble* Copyright (C) Krzysztof Bosak, 1999-10-02...1999-10-03
// All rights reserved.
// kbosak@box43.pl
// http://www.kbosak.prv.pl

#include "l3_array.h"
#include "l3_rand.h"

class SpinEnsemble
{
private:
	int _size;
	basearray<char> _spin;

public:
	inline explicit SpinEnsemble(int n);
	inline int SpinsTotal() const;
	inline int SpinsUp() const;
	inline int SpinsDown() const;
	inline double Magnetisation() const;
	inline char operator[](int pos) const;
	inline char& operator[](int pos);
	inline void SetUp();
	inline void SetDown();
	inline void ReverseUp(int n);
	inline void ReverseDown(int n);
	inline void Fill(int spinup);
};

SpinEnsemble::SpinEnsemble(int size)
	: _size(size),
	_spin(size)
{// creates ensemble of 'size' spins
}
int SpinEnsemble::SpinsTotal() const
{// returns number of spins in set
	return _size;
}
int SpinEnsemble::SpinsUp() const
{// counts number of spins-up
	int up=0;
	for(int i=0; i<_size; i++)
	{
		if(_spin[i]==1)
		{
			up++;
		}
	}
	return up;
}
int SpinEnsemble::SpinsDown() const
{// counts number of spins-down
	return _size-SpinsUp();
}
char SpinEnsemble::operator[](int pos) const
{// lookup at spins without modyfying them
	assert(pos>=0);
	assert(pos<_size);
	return _spin[pos];
}
char& SpinEnsemble::operator[](int pos)
{// lookup at spins with modifications allowed
	assert(pos>=0);
	assert(pos<_size);
	return _spin[pos];
}
void SpinEnsemble::SetUp()
{// sets all spins to spin-up
	_spin.fill(1);
}
void SpinEnsemble::SetDown()
{// sets all spins to spin-down
	_spin.erase();
}
void SpinEnsemble::ReverseUp(int n)
{// reverses n spins-down to spin-up
 // warning: could be impossible! (f.ex if all spins are already spin-up)
	assert(n>=0);
	int inversions=0;
	while(inversions<n)
	{
		const int rand_index=RAND.irandom(_size);
		if(_spin[rand_index]==0)
		{
			_spin[rand_index]=1;
			inversions++;
		}
	}
}
void SpinEnsemble::ReverseDown(int n)
{// reverses n spins-up to spin-down
 // warning: could be impossible! (f.ex if all spins are already spin-down)
	assert(n>=0);
	int inversions=0;
	while(inversions<n)
	{
		const int rand_index=RAND.irandom(_size);
		if(_spin[rand_index]==1)
		{
			_spin[rand_index]=0;
			inversions++;
		}
	}
}
void SpinEnsemble::Fill(int spinup)
{// sets n spins to spin-up and the rest to spin-down
	assert(spinup>=0);
	assert(spinup<=_size);
	if(spinup < _size/2)
	{
		SetDown();
		ReverseUp(spinup);
	}
	else
	{
		SetUp();
		ReverseDown(_size-spinup);
	}
}
double SpinEnsemble::Magnetisation() const
{// returns abs value of ( spins-up - spins-down ) / all_spins
	const double mag=2.*SpinsUp()/SpinsTotal()-1;
	if(mag<0)
	{
		return -mag;
	}
	return mag;
}

#endif //_INCLUDED_L3_SPIN_H