/* Steven Andrews, 12/02, 1/03 */
/* See documentation called Ising doc */
/* Copyright 2003 by Steven Andrews.  Permission is granted
   for non-commercial use of and modifications to the code. */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "Ising.h"

#define RGTS(s) (s)
#define DWNS(s) (s&~3|s+1&3)
#define LFTS(s) (s&~3|s+2&3)
#define UPPS(s) (s&~3|s+3&3)
#define RGTH(s) (s)
#define RDWH(s) (s&~7|((s&7)+1)%6)
#define LDWH(s) (s&~7|((s&7)+2)%6)
#define LFTH(s) (s&~7|((s&7)+3)%6)
#define LUPH(s) (s&~7|((s&7)+4)%6)
#define RUPH(s) (s&~7|((s&7)+5)%6)
#define RGTT(s) (s)
#define LDWT(s) (s&~3|((s&3)+1)%3)
#define LUPT(s) (s&~3|((s&3)+2)%3)

int clsizeS(int *g,int nx,int ny,int ix,int iy);
int clsizeH(int *g,int nx,int ny,int ix,int iy);
int clsizeT(int *g,int nx,int ny,int ix,int iy);


Lattice SetUpIsing(char shape,int nx,int ny,int spc,int *nunit,int rot) {
	int i,ix,iy,s,*g,ntot,sm,max;
	Lattice lt;

	if(shape!='S'&&shape!='H'&&shape!='T') return NULL;
	if(nx<1||nx>RAND_MAX+1||ny<1||ny>RAND_MAX+1) return NULL;
	if(spc<1) return NULL;
	if(nunit) for(ntot=0,s=0;s<spc;s++) ntot+=nunit[s];
	else ntot=0;
	if(nx*ny<ntot) return NULL;
	if(rot<0||rot>1) return NULL;

	lt=(Lattice)malloc(sizeof(struct LatticeType));
	if(!lt) return NULL;
	lt->g=NULL;
	lt->shape=shape;
	lt->nx=nx;
	lt->ny=ny;
	lt->spc=spc;
	lt->rot=rot;
	lt->g=g=(int*)calloc(nx*ny,sizeof(int));
	if(!g) {free(lt);return NULL;}
	if(!nunit) {
		for(i=0;i<nx*ny;i++) g[i]=0;
		return lt; }
	sm=-1;
	max=nx*ny-ntot;
	for(s=0;s<spc;s++)
		if(nunit[s]>max) {
			max=nunit[s];
			sm=s; }
	for(i=0;i<nx*ny;i++) g[i]=sm+1;
	for(s=-1;s<spc;s++)
		if(s!=sm)
			for(i=0;i<(s==-1?nx*ny-ntot:nunit[s]);)
				if(g[(iy=rand()%ny)*nx+(ix=rand()%nx)]==sm+1) {
					g[iy*nx+ix]=s+1;
					i++; }
	if(rot) {
		if(shape=='S') rot=4;
		else if(shape=='H') rot=6;
		else rot=3;
		for(i=0;i<nx*ny;i++) {
			g[i]*=shape=='H'?8:4;
			if(g[i]) g[i]+=rand()%rot; }}
	return lt; }


void FreeIsing(Lattice lt) {
	if(!lt) return;
	if(lt->g) free(lt->g);
	free(lt);
	return; }


void DisplayIsing(Lattice lt) {
	int ix,iy,nx,ny,*g,spc,rot;
	char c,shape;

	if(!lt||!lt->g) {
		printf("NULL lattice.\n");
		return; }
	g=lt->g;
	nx=lt->nx;
	ny=lt->ny;
	shape=lt->shape;
	spc=lt->spc;
	if(lt->rot) rot=shape=='H'?8:4;
	else rot=0;
	if(spc==1) c='x'-1;
	else c='a'-1;
	for(ix=0;ix<(rot?3*nx:2*nx);ix++) putchar('-');
	putchar('\n');
	for(iy=0;iy<ny;iy++) {
		if(shape=='H'&&iy&1||shape=='T'&&iy&2) putchar(' ');
		for(ix=0;ix<nx;ix++)
			if(g[iy*nx+ix]&&rot) printf("%c%i ",c+g[iy*nx+ix]/rot,g[iy*nx+ix]%rot);
			else if(rot) printf("   ");
			else if(g[iy*nx+ix]) printf("%c ",c+g[iy*nx+ix]);
			else printf("  ");
		putchar('\n'); }
	for(ix=0;ix<(rot?3*nx:2*nx);ix++) putchar('-');
	putchar('\n');
	return; }


void MetropolisIsingS(Lattice lt,int itmax,float *eps) {
	static int hrz[8]={1,1,0,-1,-1,-1,0,1};
	static int vrt[8]={0,1,1,1,0,-1,-1,-1};
	static int hrzo[16]={1,1,1,0,0,-1,-1,-1,-1,0,0,1,0,0,0,0};
	static int vrto[16]={0,0,1,1,1,1,0,0,-1,-1,-1,-1,0,0,0,0};
	int rnd,ix,iy,jx,jy,b,nx,ny,*g,small;
	int is,js,is2,js2,isspc,jsspc,spc,rot;
	int kir,kid,kil,kiu,kjr,kjd,kjl,kju;
	float de;
	int prob[9],db;

	g=lt->g;
	nx=lt->nx;
	ny=lt->ny;
	if((nx*ny&-nx*ny)!=nx*ny) return;
	spc=lt->spc;
	rot=lt->rot;
	small=nx*ny<=RAND_MAX+1?nx*ny-1:0;
	nx--;
	ny--;
	for(b=0;nx>>b&1;b+=1);

	if(rot) {
		spc=(spc+1)*4;
		while(itmax--) {
			if(small) {
				while(!g[ix=rand()&small]);
				iy=ix>>b;
				ix&=nx; }
			else while(!g[(iy=rand()&ny)<<b|(ix=(rand()&nx))]);
			is=g[iy<<b|ix];
			kir=g[iy<<b|(ix+1)&nx];
			kid=g[((iy+1)&ny)<<b|ix];
			kil=g[iy<<b|(ix-1)&nx];
			kiu=g[((iy-1)&ny)<<b|ix];
			rnd=rand();
			if((rnd&12)==12) {													// S, rot, no translation
				de=-eps[RGTS(is)*spc+LFTS(kir)];
				de-=eps[DWNS(is)*spc+UPPS(kid)];
				de-=eps[LFTS(is)*spc+RGTS(kil)];
				de-=eps[UPPS(is)*spc+DWNS(kiu)];
				g[iy<<b|ix]=is2=is&~3|rnd>>4&3;
				de+=eps[RGTS(is2)*spc+LFTS(kir)];
				de+=eps[DWNS(is2)*spc+UPPS(kid)];
				de+=eps[LFTS(is2)*spc+RGTS(kil)];
				de+=eps[UPPS(is2)*spc+DWNS(kiu)];
				if(de>0&&rand()>RAND_MAX*exp(-de))
					g[iy<<b|ix]=is; }
			else {																			// S, rot, with translation
				jx=(ix+hrzo[rnd&15])&nx;
				jy=(iy+vrto[rnd&15])&ny;
				js=g[jy<<b|jx];
				kjr=g[jy<<b|(jx+1)&nx];
				kjd=g[((jy+1)&ny)<<b|jx];
				kjl=g[jy<<b|(jx-1)&nx];
				kju=g[((jy-1)&ny)<<b|jx];
				de=-eps[RGTS(is)*spc+LFTS(kir)];
				de-=eps[DWNS(is)*spc+UPPS(kid)];
				de-=eps[LFTS(is)*spc+RGTS(kil)];
				de-=eps[UPPS(is)*spc+DWNS(kiu)];
				de-=eps[RGTS(js)*spc+LFTS(kjr)];
				de-=eps[DWNS(js)*spc+UPPS(kjd)];
				de-=eps[LFTS(js)*spc+RGTS(kjl)];
				de-=eps[UPPS(js)*spc+DWNS(kju)];
				g[iy<<b|ix]=is2=js?(js&~3|rnd>>4&3):0;
				g[jy<<b|jx]=js2=is&~3|rnd>>6&3;
				kir=g[iy<<b|(ix+1)&nx];
				kid=g[((iy+1)&ny)<<b|ix];
				kil=g[iy<<b|(ix-1)&nx];
				kiu=g[((iy-1)&ny)<<b|ix];
				kjr=g[jy<<b|(jx+1)&nx];
				kjd=g[((jy+1)&ny)<<b|jx];
				kjl=g[jy<<b|(jx-1)&nx];
				kju=g[((jy-1)&ny)<<b|jx];
				de+=eps[RGTS(is2)*spc+LFTS(kir)];
				de+=eps[DWNS(is2)*spc+UPPS(kid)];
				de+=eps[LFTS(is2)*spc+RGTS(kil)];
				de+=eps[UPPS(is2)*spc+DWNS(kiu)];
				de+=eps[RGTS(js2)*spc+LFTS(kjr)];
				de+=eps[DWNS(js2)*spc+UPPS(kjd)];
				de+=eps[LFTS(js2)*spc+RGTS(kjl)];
				de+=eps[UPPS(js2)*spc+DWNS(kju)];
				if(de>0&&rand()>RAND_MAX*exp(-de)) {
					g[iy<<b|ix]=is;
					g[jy<<b|jx]=js; }}}}
	else if(spc>1) {																// S, no rotation
		spc++;
		while(itmax--) {
			if(small) {
				while(!g[ix=rand()&small]);
				iy=ix>>b;
				ix&=nx; }
			else while(!g[(iy=rand()&ny)<<b|(ix=(rand()&nx))]);
			is=g[iy<<b|ix];
			rnd=rand()&7;
			jx=(ix+hrz[rnd])&nx;
			jy=(iy+vrt[rnd])&ny;
			js=g[jy<<b|jx];
			isspc=is*spc;
			jsspc=js*spc;
			de=-eps[isspc+g[iy<<b|(ix+1)&nx]];
			de-=eps[isspc+g[((iy+1)&ny)<<b|ix]];
			de-=eps[isspc+g[iy<<b|(ix-1)&nx]];
			de-=eps[isspc+g[((iy-1)&ny)<<b|ix]];
			de-=eps[jsspc+g[jy<<b|(jx+1)&nx]];
			de-=eps[jsspc+g[((jy+1)&ny)<<b|jx]];
			de-=eps[jsspc+g[jy<<b|(jx-1)&nx]];
			de-=eps[jsspc+g[((jy-1)&ny)<<b|jx]];
			g[iy<<b|ix]=js;
			g[jy<<b|jx]=is;
			de+=eps[isspc+g[jy<<b|(jx+1)&nx]];
			de+=eps[isspc+g[((jy+1)&ny)<<b|jx]];
			de+=eps[isspc+g[jy<<b|(jx-1)&nx]];
			de+=eps[isspc+g[((jy-1)&ny)<<b|jx]];
			de+=eps[jsspc+g[iy<<b|(ix+1)&nx]];
			de+=eps[jsspc+g[((iy+1)&ny)<<b|ix]];
			de+=eps[jsspc+g[iy<<b|(ix-1)&nx]];
			de+=eps[jsspc+g[((iy-1)&ny)<<b|ix]];
			if(de>0&&rand()>RAND_MAX*exp(-de)) {
				g[iy<<b|ix]=is;
				g[jy<<b|jx]=js; }}}
	else {																					// S, no rotation, no species
		for(db=-4;db<=4;db++) prob[db+4]=exp(-db**eps)<1?RAND_MAX*exp(-db**eps):RAND_MAX;
		while(itmax--) {
			if(small) {
				while(!g[ix=rand()&small]);
				iy=ix>>b;
				ix&=nx; }
			else while(!g[(iy=rand()&ny)<<b|(ix=(rand()&nx))]);
			rnd=rand()&7;
			jx=(ix+hrz[rnd])&nx;
			jy=(iy+vrt[rnd])&ny;
			if(!g[jy<<b|jx]) {
				g[iy<<b|ix]=0;
				db=-g[iy<<b|(ix+1)&nx];
				db-=g[((iy+1)&ny)<<b|ix];
				db-=g[iy<<b|(ix-1)&nx];
				db-=g[((iy-1)&ny)<<b|ix];
				db+=g[jy<<b|(jx+1)&nx];
				db+=g[((jy+1)&ny)<<b|jx];
				db+=g[jy<<b|(jx-1)&nx];
				db+=g[((jy-1)&ny)<<b|jx];
				if(prob[db+4]==RAND_MAX||rand()<=prob[db+4]) g[jy<<b|jx]=1;
				else g[iy<<b|ix]=1; }}}
	return; }


void MetropolisIsingH(Lattice lt,int itmax,float *eps) {
	static int hrz[16]={0,0, 0,0, 1,1, 0,1, -1,0, -1,-1, -1,0, 0,1};
	static int vrt[16]={0,0, 0,0, 0,0, 1,1, 1,1, 0,0, -1,-1, -1,-1};
	int rnd,ix,iy,jx,jy,nx,ny,b,*g,small;
	int is,is2,js,js2,spc,isspc,jsspc,rot;
	int ki0,ki1,ki2,ki3,ki4,ki5,kj0,kj1,kj2,kj3,kj4,kj5;
	float de;
	int prob[13],db;

	g=lt->g;
	nx=lt->nx;
	ny=lt->ny;
	if((nx*ny&-nx*ny)!=nx*ny) return;
	spc=lt->spc;
	rot=lt->rot;
	small=nx*ny<=RAND_MAX+1?nx*ny-1:0;
	nx--;
	ny--;
	for(b=0;nx>>b&1;b+=1);

	if(rot) {
		spc=(spc+1)*8;
		while(itmax--) {
			if(small) {
				while(!g[ix=rand()&small]);
				iy=ix>>b;
				ix&=nx; }
			else while(!g[(iy=rand()&ny)<<b|(ix=(rand()&nx))]);
			is=g[iy<<b|ix];
			if(iy&1) {
				ki0=g[iy<<b|(ix+1)&nx];
				ki1=g[((iy+1)&ny)<<b|(ix+1)&nx];
				ki2=g[((iy+1)&ny)<<b|ix];
				ki3=g[iy<<b|(ix-1)&nx];
				ki4=g[((iy-1)&ny)<<b|ix];
				ki5=g[((iy-1)&ny)<<b|(ix+1)&nx];
				rnd=rand()|1; }
			else {
				ki0=g[iy<<b|(ix+1)&nx];
				ki1=g[((iy+1)&ny)<<b|ix];
				ki2=g[((iy+1)&ny)<<b|(ix-1)&nx];
				ki3=g[iy<<b|(ix-1)&nx];
				ki4=g[((iy-1)&ny)<<b|(ix-1)&nx];
				ki5=g[((iy-1)&ny)<<b|ix];
				rnd=rand()&~1; }
			if(!(rnd&12)) {															// H, rot, no translation
				de=-eps[RGTH(is)*spc+LFTH(ki0)];
				de-=eps[RDWH(is)*spc+LUPH(ki1)];
				de-=eps[LDWH(is)*spc+RUPH(ki2)];
				de-=eps[LFTH(is)*spc+RGTH(ki3)];
				de-=eps[LUPH(is)*spc+RDWH(ki4)];
				de-=eps[RUPH(is)*spc+LDWH(ki5)];
				g[iy<<b|ix]=is2=is&~7|(rnd>>4)%6;
				de+=eps[RGTH(is2)*spc+LFTH(ki0)];
				de+=eps[RDWH(is2)*spc+LUPH(ki1)];
				de+=eps[LDWH(is2)*spc+RUPH(ki2)];
				de+=eps[LFTH(is2)*spc+RGTH(ki3)];
				de+=eps[LUPH(is2)*spc+RDWH(ki4)];
				de+=eps[RUPH(is2)*spc+LDWH(ki5)];
				if(de>0&&rand()>RAND_MAX*exp(-de))
					g[iy<<b|ix]=is; }
			else {																			// H, rot, with translation
				jx=(ix+hrz[rnd&15])&nx;
				jy=(iy+vrt[rnd&15])&ny;
				js=g[jy<<b|jx];
				if(jy&1) {
					kj0=g[jy<<b|(jx+1)&nx];
					kj1=g[((jy+1)&ny)<<b|(jx+1)&nx];
					kj2=g[((jy+1)&ny)<<b|jx];
					kj3=g[jy<<b|(jx-1)&nx];
					kj4=g[((jy-1)&ny)<<b|jx];
					kj5=g[((jy-1)&ny)<<b|(jx+1)&nx]; }
				else {
					kj0=g[jy<<b|(jx+1)&nx];
					kj1=g[((jy+1)&ny)<<b|jx];
					kj2=g[((jy+1)&ny)<<b|(jx-1)&nx];
					kj3=g[jy<<b|(jx-1)&nx];
					kj4=g[((jy-1)&ny)<<b|(jx-1)&nx];
					kj5=g[((jy-1)&ny)<<b|jx]; }
				de=-eps[RGTH(is)*spc+LFTH(ki0)];
				de-=eps[RDWH(is)*spc+LUPH(ki1)];
				de-=eps[LDWH(is)*spc+RUPH(ki2)];
				de-=eps[LFTH(is)*spc+RGTH(ki3)];
				de-=eps[LUPH(is)*spc+RDWH(ki4)];
				de-=eps[RUPH(is)*spc+LDWH(ki5)];
				de-=eps[RGTH(js)*spc+LFTH(kj0)];
				de-=eps[RDWH(js)*spc+LUPH(kj1)];
				de-=eps[LDWH(js)*spc+RUPH(kj2)];
				de-=eps[LFTH(js)*spc+RGTH(kj3)];
				de-=eps[LUPH(js)*spc+RDWH(kj4)];
				de-=eps[RUPH(js)*spc+LDWH(kj5)];
				g[iy<<b|ix]=is2=js?js&~7|(rnd>>4)%6:0;
				g[jy<<b|jx]=js2=is&~7|(rnd>>6)%6;
				if(iy&1) {
					ki0=g[iy<<b|(ix+1)&nx];
					ki1=g[((iy+1)&ny)<<b|(ix+1)&nx];
					ki2=g[((iy+1)&ny)<<b|ix];
					ki3=g[iy<<b|(ix-1)&nx];
					ki4=g[((iy-1)&ny)<<b|ix];
					ki5=g[((iy-1)&ny)<<b|(ix+1)&nx]; }
				else {
					ki0=g[iy<<b|(ix+1)&nx];
					ki1=g[((iy+1)&ny)<<b|ix];
					ki2=g[((iy+1)&ny)<<b|(ix-1)&nx];
					ki3=g[iy<<b|(ix-1)&nx];
					ki4=g[((iy-1)&ny)<<b|(ix-1)&nx];
					ki5=g[((iy-1)&ny)<<b|ix]; }
				if(jy&1) {
					kj0=g[jy<<b|(jx+1)&nx];
					kj1=g[((jy+1)&ny)<<b|(jx+1)&nx];
					kj2=g[((jy+1)&ny)<<b|jx];
					kj3=g[jy<<b|(jx-1)&nx];
					kj4=g[((jy-1)&ny)<<b|jx];
					kj5=g[((jy-1)&ny)<<b|(jx+1)&nx]; }
				else {
					kj0=g[jy<<b|(jx+1)&nx];
					kj1=g[((jy+1)&ny)<<b|jx];
					kj2=g[((jy+1)&ny)<<b|(jx-1)&nx];
					kj3=g[jy<<b|(jx-1)&nx];
					kj4=g[((jy-1)&ny)<<b|(jx-1)&nx];
					kj5=g[((jy-1)&ny)<<b|jx]; }
				de+=eps[RGTH(is2)*spc+LFTH(ki0)];
				de+=eps[RDWH(is2)*spc+LUPH(ki1)];
				de+=eps[LDWH(is2)*spc+RUPH(ki2)];
				de+=eps[LFTH(is2)*spc+RGTH(ki3)];
				de+=eps[LUPH(is2)*spc+RDWH(ki4)];
				de+=eps[RUPH(is2)*spc+LDWH(ki5)];
				de+=eps[RGTH(js2)*spc+LFTH(kj0)];
				de+=eps[RDWH(js2)*spc+LUPH(kj1)];
				de+=eps[LDWH(js2)*spc+RUPH(kj2)];
				de+=eps[LFTH(js2)*spc+RGTH(kj3)];
				de+=eps[LUPH(js2)*spc+RDWH(kj4)];
				de+=eps[RUPH(js2)*spc+LDWH(kj5)];
				if(de>0&&rand()>RAND_MAX*exp(-de)) {
					g[iy<<b|ix]=is;
					g[jy<<b|jx]=js; }}}}
	else if(spc>1) {																// H, no rotation
		spc++;
		while(itmax--) {
			if(small) {
				while(!g[ix=rand()&small]);
				iy=ix>>b;
				ix&=nx; }
			else while(!g[(iy=rand()&ny)<<b|(ix=(rand()&nx))]);
			is=g[iy<<b|ix];
			rnd=rand()&14|iy&1;
			jx=(ix+hrz[rnd])&nx;
			jy=(iy+vrt[rnd])&ny;
			js=g[jy<<b|jx];
			isspc=is*spc;
			jsspc=js*spc;
			de=-eps[isspc+g[iy<<b|(ix+1)&nx]];
			de-=eps[isspc+g[iy<<b|(ix-1)&nx]];
			de-=eps[isspc+g[((iy+1)&ny)<<b|ix]];
			de-=eps[isspc+g[((iy-1)&ny)<<b|ix]];
			de-=eps[isspc+g[((iy-1)&ny)<<b|(ix+(iy&1?1:-1))&nx]];
			de-=eps[isspc+g[((iy+1)&ny)<<b|(ix+(iy&1?1:-1))&nx]];
			de-=eps[jsspc+g[jy<<b|(jx+1)&nx]];
			de-=eps[jsspc+g[jy<<b|(jx-1)&nx]];
			de-=eps[jsspc+g[((jy+1)&ny)<<b|jx]];
			de-=eps[jsspc+g[((jy-1)&ny)<<b|jx]];
			de-=eps[jsspc+g[((jy-1)&ny)<<b|(jx+(jy&1?1:-1))&nx]];
			de-=eps[jsspc+g[((jy+1)&ny)<<b|(jx+(jy&1?1:-1))&nx]];
			g[iy<<b|ix]=js;
			g[jy<<b|jx]=is;
			de+=eps[isspc+g[jy<<b|(jx+1)&nx]];
			de+=eps[isspc+g[jy<<b|(jx-1)&nx]];
			de+=eps[isspc+g[((jy+1)&ny)<<b|jx]];
			de+=eps[isspc+g[((jy-1)&ny)<<b|jx]];
			de+=eps[isspc+g[((jy-1)&ny)<<b|(jx+(jy&1?1:-1))&nx]];
			de+=eps[isspc+g[((jy+1)&ny)<<b|(jx+(jy&1?1:-1))&nx]];
			de+=eps[jsspc+g[iy<<b|(ix+1)&nx]];
			de+=eps[jsspc+g[iy<<b|(ix-1)&nx]];
			de+=eps[jsspc+g[((iy+1)&ny)<<b|ix]];
			de+=eps[jsspc+g[((iy-1)&ny)<<b|ix]];
			de+=eps[jsspc+g[((iy-1)&ny)<<b|(ix+(iy&1?1:-1))&nx]];
			de+=eps[jsspc+g[((iy+1)&ny)<<b|(ix+(iy&1?1:-1))&nx]];
			if(de>0&&rand()>RAND_MAX*exp(-de)) {
				g[iy<<b|ix]=is;
				g[jy<<b|jx]=js; }}}
	else {																					// H, no rotation, no species
		for(db=-6;db<=6;db++) prob[db+6]=exp(-db**eps)<1?RAND_MAX*exp(-db**eps):RAND_MAX;
		while(itmax--) {
			if(small) {
				while(!g[ix=rand()&small]);
				iy=ix>>b;
				ix&=nx; }
			else while(!g[(iy=rand()&ny)<<b|(ix=(rand()&nx))]);
			rnd=rand()&14|iy&1;
			jx=(ix+hrz[rnd])&nx;
			jy=(iy+vrt[rnd])&ny;
			if(!g[jy<<b|jx]) {
				g[iy<<b|ix]=0;
				db =g[jy<<b|(jx+1)&nx];
				db-=g[iy<<b|(ix+1)&nx];
				db+=g[jy<<b|(jx-1)&nx];
				db-=g[iy<<b|(ix-1)&nx];
				db+=g[((jy+1)&ny)<<b|jx];
				db-=g[((iy+1)&ny)<<b|ix];
				db+=g[((jy-1)&ny)<<b|jx];
				db-=g[((iy-1)&ny)<<b|ix];
				db+=g[((jy+1)&ny)<<b|(jx+(jy&1?1:-1))&nx];
				db-=g[((iy+1)&ny)<<b|(ix+(iy&1?1:-1))&nx];
				db+=g[((jy-1)&ny)<<b|(jx+(jy&1?1:-1))&nx];
				db-=g[((iy-1)&ny)<<b|(ix+(iy&1?1:-1))&nx];
				if(prob[db+6]==RAND_MAX||rand()<=prob[db+6]) g[jy<<b|jx]=1;
				else g[iy<<b|ix]=1; }}}
	return; }


void MetropolisIsingT(Lattice lt,int itmax,float *eps) {
	static int hrz[16]={0,0,0,0, 0,0,0,0, 0,0,0,0, -1,-1,1,1};
	static int vrt[16]={0,0,0,0, -1,-1,-1,-1, 1,1,1,1, -1,1,-1,1};
	int rnd,ix,iy,jx,jy,b,nx,ny,*g,small;
	int is,js,is2,js2,spc,isspc,jsspc,rot;
	int ki0,ki1,ki2,kj0,kj1,kj2;
	float de;
	int prob[7],db;

	g=lt->g;
	nx=lt->nx;
	ny=lt->ny;
	if((nx*ny&-nx*ny)!=nx*ny) return;
	spc=lt->spc;
	rot=lt->rot;
	small=nx*ny<=RAND_MAX+1?nx*ny-1:0;
	nx--;
	ny--;
	for(b=0;nx>>b&1;b+=1);

	if(rot) {
		spc=(spc+1)*4;
		while(itmax--) {
			if(small) {
				while(!g[ix=rand()&small]);
				iy=ix>>b;
				ix&=nx; }
			else while(!g[(iy=rand()&ny)<<b|(ix=(rand()&nx))]);
			is=g[iy<<b|ix];
			if((iy&3)==0) {
				ki0=g[((iy-1)&ny)<<b|ix];
				ki1=g[((iy+1)&ny)<<b|ix];
				ki2=g[((iy-1)&ny)<<b|(ix-1)&nx];
				rnd=rand()&~3; }
			else if((iy&3)==1) {
				ki0=g[((iy+1)&ny)<<b|ix];
				ki1=g[((iy+1)&ny)<<b|(ix-1)&nx];
				ki2=g[((iy-1)&ny)<<b|ix];
				rnd=rand()&~3|1; }
			else if((iy&3)==2) {
				ki0=g[((iy-1)&ny)<<b|(ix+1)&nx];
				ki1=g[((iy+1)&ny)<<b|ix];
				ki2=g[((iy-1)&ny)<<b|ix];
				rnd=rand()&~3|2; }
			else {
				ki0=g[((iy+1)&ny)<<b|(ix+1)&nx];
				ki1=g[((iy+1)&ny)<<b|ix];
				ki2=g[((iy-1)&ny)<<b|ix];
				rnd=rand()|3; }
			if(!(rnd&12)) {															// T, rot, no translation
				if(iy&1) {
					de=-eps[RGTT(is)*spc+LUPT(ki0)];
					de-=eps[LDWT(is)*spc+RGTT(ki1)];
					de-=eps[LUPT(is)*spc+LDWT(ki2)]; }
				else {
					de=-eps[RGTT(is)*spc+LDWT(ki0)];
					de-=eps[LDWT(is)*spc+LUPT(ki1)];
					de-=eps[LUPT(is)*spc+RGTT(ki2)]; }
				g[iy<<b|ix]=is2=is&~3|(rnd>>4)%3;
				if(iy&1) {
					de+=eps[RGTT(is2)*spc+LUPT(ki0)];
					de+=eps[LDWT(is2)*spc+RGTT(ki1)];
					de+=eps[LUPT(is2)*spc+LDWT(ki2)]; }
				else {
					de+=eps[RGTT(is2)*spc+LDWT(ki0)];
					de+=eps[LDWT(is2)*spc+LUPT(ki1)];
					de+=eps[LUPT(is2)*spc+RGTT(ki2)]; }
				if(de>0&&rand()>RAND_MAX*exp(-de))
					g[iy<<b|ix]=is; }
			else {																			// T, rot, with translation
				jx=(ix+hrz[rnd&15])&nx;
				jy=(iy+vrt[rnd&15])&ny;
				js=g[jy<<b|jx];
				if((jy&3)==0) {
					kj0=g[((jy-1)&ny)<<b|jx];
					kj1=g[((jy+1)&ny)<<b|jx];
					kj2=g[((jy-1)&ny)<<b|(jx-1)&nx]; }
				else if((jy&3)==1) {
					kj0=g[((jy+1)&ny)<<b|jx];
					kj1=g[((jy+1)&ny)<<b|(jx-1)&nx];
					kj2=g[((jy-1)&ny)<<b|jx]; }
				else if((jy&3)==2) {
					kj0=g[((jy-1)&ny)<<b|(jx+1)&nx];
					kj1=g[((jy+1)&ny)<<b|jx];
					kj2=g[((jy-1)&ny)<<b|jx]; }
				else {
					kj0=g[((jy+1)&ny)<<b|(jx+1)&nx];
					kj1=g[((jy+1)&ny)<<b|jx];
					kj2=g[((jy-1)&ny)<<b|jx]; }
				if(iy&1) {
					de=-eps[RGTT(is)*spc+LUPT(ki0)];
					de-=eps[LDWT(is)*spc+RGTT(ki1)];
					de-=eps[LUPT(is)*spc+LDWT(ki2)]; }
				else {
					de=-eps[RGTT(is)*spc+LDWT(ki0)];
					de-=eps[LDWT(is)*spc+LUPT(ki1)];
					de-=eps[LUPT(is)*spc+RGTT(ki2)]; }
				if(jy&1) {
					de-=eps[RGTT(js)*spc+LUPT(kj0)];
					de-=eps[LDWT(js)*spc+RGTT(kj1)];
					de-=eps[LUPT(js)*spc+LDWT(kj2)]; }
				else {
					de-=eps[RGTT(js)*spc+LDWT(kj0)];
					de-=eps[LDWT(js)*spc+LUPT(kj1)];
					de-=eps[LUPT(js)*spc+RGTT(kj2)]; }
				g[iy<<b|ix]=is2=js?js&~3|(rnd>>4)%3:0;
				g[jy<<b|jx]=js2=is&~3|(rnd>>6)%3;
				if((iy&3)==0) {
					ki0=g[((iy-1)&ny)<<b|ix];
					ki1=g[((iy+1)&ny)<<b|ix];
					ki2=g[((iy-1)&ny)<<b|(ix-1)&nx]; }
				else if((iy&3)==1) {
					ki0=g[((iy+1)&ny)<<b|ix];
					ki1=g[((iy+1)&ny)<<b|(ix-1)&nx];
					ki2=g[((iy-1)&ny)<<b|ix]; }
				else if((iy&3)==2) {
					ki0=g[((iy-1)&ny)<<b|(ix+1)&nx];
					ki1=g[((iy+1)&ny)<<b|ix];
					ki2=g[((iy-1)&ny)<<b|ix]; }
				else {
					ki0=g[((iy+1)&ny)<<b|(ix+1)&nx];
					ki1=g[((iy+1)&ny)<<b|ix];
					ki2=g[((iy-1)&ny)<<b|ix]; }
				if((jy&3)==0) {
					kj0=g[((jy-1)&ny)<<b|jx];
					kj1=g[((jy+1)&ny)<<b|jx];
					kj2=g[((jy-1)&ny)<<b|(jx-1)&nx]; }
				else if((jy&3)==1) {
					kj0=g[((jy+1)&ny)<<b|jx];
					kj1=g[((jy+1)&ny)<<b|(jx-1)&nx];
					kj2=g[((jy-1)&ny)<<b|jx]; }
				else if((jy&3)==2) {
					kj0=g[((jy-1)&ny)<<b|(jx+1)&nx];
					kj1=g[((jy+1)&ny)<<b|jx];
					kj2=g[((jy-1)&ny)<<b|jx]; }
				else {
					kj0=g[((jy+1)&ny)<<b|(jx+1)&nx];
					kj1=g[((jy+1)&ny)<<b|jx];
					kj2=g[((jy-1)&ny)<<b|jx]; }
				if(iy&1) {
					de+=eps[RGTT(is2)*spc+LUPT(ki0)];
					de+=eps[LDWT(is2)*spc+RGTT(ki1)];
					de+=eps[LUPT(is2)*spc+LDWT(ki2)]; }
				else {
					de+=eps[RGTT(is2)*spc+LDWT(ki0)];
					de+=eps[LDWT(is2)*spc+LUPT(ki1)];
					de+=eps[LUPT(is2)*spc+RGTT(ki2)]; }
				if(jy&1) {
					de+=eps[RGTT(js2)*spc+LUPT(kj0)];
					de+=eps[LDWT(js2)*spc+RGTT(kj1)];
					de+=eps[LUPT(js2)*spc+LDWT(kj2)]; }
				else {
					de+=eps[RGTT(js2)*spc+LDWT(kj0)];
					de+=eps[LDWT(js2)*spc+LUPT(kj1)];
					de+=eps[LUPT(js2)*spc+RGTT(kj2)]; }
				if(de>0&&rand()>RAND_MAX*exp(-de)) {
					g[iy<<b|ix]=is;
					g[jy<<b|jx]=js; }}}}
	else if(spc>1) {																// T, no rotation
		spc++;
		while(itmax--) {
			if(small) {
				while(!g[ix=rand()&small]);
				iy=ix>>b;
				ix&=nx; }
			else while(!g[(iy=rand()&ny)<<b|(ix=(rand()&nx))]);
			is=g[iy<<b|ix];
			rnd=rand()&12|iy&3;
			jx=(ix+hrz[rnd])&nx;
			jy=(iy+vrt[rnd])&ny;
			js=g[jy<<b|jx];
			isspc=is*spc;
			jsspc=js*spc;
			de=-eps[isspc+g[((iy+1)&ny)<<b|ix]];
			de-=eps[isspc+g[((iy-1)&ny)<<b|ix]];
			de-=eps[isspc+g[((iy+(iy&1?1:-1))&ny)<<b|(ix+(iy&2?1:-1))&nx]];
			de-=eps[jsspc+g[((jy+1)&ny)<<b|jx]];
			de-=eps[jsspc+g[((jy-1)&ny)<<b|jx]];
			de-=eps[jsspc+g[((jy+(jy&1?1:-1))&ny)<<b|(jx+(jy&2?1:-1))&nx]];
			g[iy<<b|ix]=js;
			g[jy<<b|jx]=is;
			de+=eps[isspc+g[((jy+1)&ny)<<b|jx]];
			de+=eps[isspc+g[((jy-1)&ny)<<b|jx]];
			de+=eps[isspc+g[((jy+(jy&1?1:-1))&ny)<<b|(jx+(jy&2?1:-1))&nx]];
			de+=eps[jsspc+g[((iy+1)&ny)<<b|ix]];
			de+=eps[jsspc+g[((iy-1)&ny)<<b|ix]];
			de+=eps[jsspc+g[((iy+(iy&1?1:-1))&ny)<<b|(ix+(iy&2?1:-1))&nx]];
			if(de>0&&rand()>RAND_MAX*exp(-de)) {
				g[iy<<b|ix]=is;
				g[jy<<b|jx]=js; }}}
	else {																					// T, no rotation, no species
		for(db=-3;db<=3;db++) prob[db+3]=exp(-db**eps)<1?RAND_MAX*exp(-db**eps):RAND_MAX;
		while(itmax--) {
			if(small) {
				while(!g[ix=rand()&small]);
				iy=ix>>b;
				ix&=nx; }
			else while(!g[(iy=rand()&ny)<<b|(ix=(rand()&nx))]);
			rnd=rand()&12|iy&3;
			jx=(ix+hrz[rnd])&nx;
			jy=(iy+vrt[rnd])&ny;
			if(!g[jy<<b|jx]) {
				g[iy<<b|ix]=0;
				db =g[((jy+1)&ny)<<b|jx];
				db-=g[((iy+1)&ny)<<b|ix];
				db+=g[((jy-1)&ny)<<b|jx];
				db-=g[((iy-1)&ny)<<b|ix];
				db+=g[((jy+(jy&1?1:-1))&ny)<<b|(jx+(jy&2?1:-1))&nx];
				db-=g[((iy+(iy&1?1:-1))&ny)<<b|(ix+(iy&2?1:-1))&nx];
				if(prob[db+3]==RAND_MAX||rand()<=prob[db+3]) g[jy<<b|jx]=1;
				else g[iy<<b|ix]=1; }}}
	return; }


void MetropolisIsing(Lattice lt,int itmax,float *eps) {
	if(lt->shape=='S') MetropolisIsingS(lt,itmax,eps);
	else if(lt->shape=='H') MetropolisIsingH(lt,itmax,eps);
	else MetropolisIsingT(lt,itmax,eps);
	return; }


float EnergyAtIsingS(Lattice lt,float *eps,int ix,int iy,int spc1,int rot1) {
	float e;
	int nx,ny,bond,is2,spc,rot,*g;
	int k0,k1,k2,k3;

	g=lt->g;
	nx=lt->nx;
	ny=lt->ny;
	spc=lt->spc;
	rot=lt->rot;

	if(rot) {
		spc=(spc+1)*4;
		is2=spc1<<2|rot1;
		k0=g[iy*nx+(ix+1)%nx];
		k1=g[(iy+1)%ny*nx+ix];
		k2=g[iy*nx+(ix-1+nx)%nx];
		k3=g[(iy-1+ny)%ny*nx+ix];
		e=eps[RGTS(is2)*spc+LFTS(k0)];
		e+=eps[DWNS(is2)*spc+UPPS(k1)];
		e+=eps[LFTS(is2)*spc+RGTS(k2)];
		e+=eps[UPPS(is2)*spc+DWNS(k3)]; }
	else if(spc>1) {
		spc++;
		is2=spc1*spc;
		e=eps[is2+g[iy*nx+(ix+1)%nx]];
		e+=eps[is2+g[(iy+1)%ny*nx+ix]];
		e+=eps[is2+g[iy*nx+(ix-1+nx)%nx]];
		e+=eps[is2+g[(iy-1+ny)%ny*nx+ix]]; }
	else if(spc1) {
		bond=g[iy*nx+(ix+1)%nx];
		bond+=g[(iy+1)%ny*nx+ix];
		bond+=g[iy*nx+(ix-1+nx)%nx];
		bond+=g[(iy-1+ny)%ny*nx+ix];
		e=*eps*bond; }
	else e=0;
	return e; }


float EnergyAtIsingH(Lattice lt,float *eps,int ix,int iy,int spc1,int rot1) {
	float e;
	int nx,ny,bond,is2,spc,rot,*g;
	int k0,k1,k2,k3,k4,k5;

	g=lt->g;
	nx=lt->nx;
	ny=lt->ny;
	spc=lt->spc;
	rot=lt->rot;

	if(rot) {
		spc=(spc+1)*8;
		is2=spc1<<3|rot1;
		if(iy&1) {
			k0=g[iy*nx+(ix+1)%nx];
			k1=g[(iy+1)%ny*nx+(ix+1)%nx];
			k2=g[(iy+1)%ny*nx+ix];
			k3=g[iy*nx+(ix-1+nx)%nx];
			k4=g[(iy-1+ny)%ny*nx+ix];
			k5=g[(iy-1+ny)%ny*nx+(ix+1)%nx]; }
		else {
			k0=g[iy*nx+(ix+1)%nx];
			k1=g[(iy+1)%ny*nx+ix];
			k2=g[(iy+1)%ny*nx+(ix-1+nx)%nx];
			k3=g[iy*nx+(ix-1+nx)%nx];
			k4=g[(iy-1+ny)%ny*nx+(ix-1+nx)%nx];
			k5=g[(iy-1+ny)%ny*nx+ix]; }
		e=eps[RGTH(is2)*spc+LFTH(k0)];
		e+=eps[RDWH(is2)*spc+LUPH(k1)];
		e+=eps[LDWH(is2)*spc+RUPH(k2)];
		e+=eps[LFTH(is2)*spc+RGTH(k3)];
		e+=eps[LUPH(is2)*spc+RDWH(k4)];
		e+=eps[RUPH(is2)*spc+LDWH(k5)]; }
	else if(spc>1) {
		spc++;
		is2=spc1*spc;
		e=eps[is2+g[iy*nx+(ix+1)%nx]];
		e+=eps[is2+g[(iy+1)%ny*nx+ix]];
		e+=eps[is2+g[iy*nx+(ix-1+nx)%nx]];
		e+=eps[is2+g[(iy-1+ny)%ny*nx+ix]];
		e+=eps[is2+g[(iy-1+ny)%ny*nx+(ix+(iy&1?1:-1)+nx)%nx]];
		e+=eps[is2+g[(iy+1)%ny*nx+(ix+(iy&1?1:-1)+nx)%nx]]; }
	else if(spc1) {
		bond=g[iy*nx+(ix+1)%nx];
		bond+=g[(iy+1)%ny*nx+ix];
		bond+=g[iy*nx+(ix-1+nx)%nx];
		bond+=g[(iy-1+ny)%ny*nx+ix];
		bond+=g[(iy-1+ny)%ny*nx+(ix+(iy&1?1:-1)+nx)%nx];
		bond+=g[(iy+1)%ny*nx+(ix+(iy&1?1:-1)+nx)%nx];
		e=*eps*bond; }
	else e=0;
	return e; }


float EnergyAtIsingT(Lattice lt,float *eps,int ix,int iy,int spc1,int rot1) {
	float e;
	int nx,ny,bond,is2,spc,rot,*g;
	int k0,k1,k2;

	g=lt->g;
	nx=lt->nx;
	ny=lt->ny;
	spc=lt->spc;
	rot=lt->rot;

	if(rot) {
		spc=(spc+1)*4;
		is2=spc1<<2|rot1;
		if((iy&3)==0) {
			k0=g[(iy-1+ny)%ny*nx+ix];
			k1=g[(iy+1)%ny*nx+ix];
			k2=g[(iy-1+ny)%ny*nx+(ix-1+nx)%nx]; }
		else if((iy&3)==1) {
			k0=g[(iy+1)%ny*nx+ix];
			k1=g[(iy+1)%ny*nx+(ix-1+nx)%nx];
			k2=g[(iy-1+ny)%ny*nx+ix]; }
		else if((iy&3)==2) {
			k0=g[(iy-1+ny)%ny*nx+(ix+1)%nx];
			k1=g[(iy+1)%ny*nx+ix];
			k2=g[(iy-1+ny)%ny*nx+ix]; }
		else {
			k0=g[(iy+1)%ny*nx+(ix+1)%nx];
			k1=g[(iy+1)%ny*nx+ix];
			k2=g[(iy-1+ny)%ny*nx+ix]; }
		if(iy&1) {
			e=eps[RGTT(is2)*spc+LUPT(k0)];
			e+=eps[LDWT(is2)*spc+RGTT(k1)];
			e+=eps[LUPT(is2)*spc+LDWT(k2)]; }
		else {
			e=eps[RGTT(is2)*spc+LDWT(k0)];
			e+=eps[LDWT(is2)*spc+LUPT(k1)];
			e+=eps[LUPT(is2)*spc+RGTT(k2)]; }}
	else if(spc>1) {
		spc++;
		is2=spc1*spc;
		e=eps[is2+g[(iy-1+ny)%ny*nx+ix]];
		e+=eps[is2+g[(iy+1)%ny*nx+ix]];
		e+=eps[is2+g[(iy+(iy&1?1:-1)+ny)%ny*nx+(ix+(iy&2?1:-1)+nx)%nx]]; }
	else if(spc1) {
		bond=g[(iy-1+ny)%ny*nx+ix];
		bond+=g[(iy+1)%ny*nx+ix];
		bond+=g[(iy+(iy&1?1:-1)+ny)%ny*nx+(ix+(iy&2?1:-1)+nx)%nx];
		e=*eps*bond; }
	else e=0;
	return e; }


float EnergyAtIsing(Lattice lt,float *eps,int ix,int iy,int spc1,int rot1) {
	if(lt->shape=='S') return EnergyAtIsingS(lt,eps,ix,iy,spc1,rot1);
	else if(lt->shape=='H') return EnergyAtIsingH(lt,eps,ix,iy,spc1,rot1);
	else return EnergyAtIsingT(lt,eps,ix,iy,spc1,rot1); }


float EnergyIsingS(Lattice lt,float *eps) {
	float e;
	int ix,iy,nx,ny,bond,is2,spc,rot,*g;
	int k0,k1;

	g=lt->g;
	nx=lt->nx;
	ny=lt->ny;
	spc=lt->spc;
	rot=lt->rot;

	if(rot) {
		spc=(spc+1)*4;
		e=0;
		for(iy=0;iy<ny;iy++)
			for(ix=0;ix<nx;ix++) {
				is2=g[iy*nx+ix];
				k0=g[iy*nx+(ix+1)%nx];
				k1=g[(iy+1)%ny*nx+ix];
				e+=eps[RGTS(is2)*spc+LFTS(k0)];
				e+=eps[DWNS(is2)*spc+UPPS(k1)]; }}
	else if(spc>1) {
		spc++;
		e=0;
		for(iy=0;iy<ny;iy++)
			for(ix=0;ix<nx;ix++)	{
				is2=g[iy*nx+ix]*spc;
				e+=eps[is2+g[iy*nx+(ix+1)%nx]];
				e+=eps[is2+g[(iy+1)%ny*nx+ix]]; }}
	else {
		bond=0;
		for(iy=0;iy<ny;iy++)
			for(ix=0;ix<nx;ix++)
				if(g[iy*nx+ix]) {
					bond+=g[iy*nx+(ix+1)%nx];
					bond+=g[(iy+1)%ny*nx+ix]; }
		e=*eps*bond; }
	return e; }


float EnergyIsingH(Lattice lt,float *eps) {
	float e;
	int ix,iy,nx,ny,bond,is2,spc,rot,*g;
	int k0,k1,k2;

	g=lt->g;
	nx=lt->nx;
	ny=lt->ny;
	spc=lt->spc;
	rot=lt->rot;

	if(rot) {
		spc=(spc+1)*8;
		e=0;
		for(iy=0;iy<ny;iy++)
			for(ix=0;ix<nx;ix++) {
				is2=g[iy*nx+ix];
				if(iy&1) {
					k0=g[iy*nx+(ix+1)%nx];
					k1=g[(iy+1)%ny*nx+(ix+1)%nx];
					k2=g[(iy+1)%ny*nx+ix]; }
				else {
					k0=g[iy*nx+(ix+1)%nx];
					k1=g[(iy+1)%ny*nx+ix];
					k2=g[(iy+1)%ny*nx+(ix-1+nx)%nx]; }
				e+=eps[RGTH(is2)*spc+LFTH(k0)];
				e+=eps[RDWH(is2)*spc+LUPH(k1)];
				e+=eps[LDWH(is2)*spc+RUPH(k2)]; }}
	else if(spc>1) {
		spc++;
		e=0;
		for(iy=0;iy<ny;iy++)
			for(ix=0;ix<nx;ix++) {
				is2=g[iy*nx+ix]*spc;
				e+=eps[is2+g[iy*nx+(ix+1)%nx]];
				e+=eps[is2+g[(iy+1)%ny*nx+ix]];
				e+=eps[is2+g[(iy+1)%ny*nx+(ix+(iy&1?1:-1)+nx)%nx]]; }}
	else {
		bond=0;
		for(iy=0;iy<ny;iy++)
			for(ix=0;ix<nx;ix++)
				if(g[iy*nx+ix]) {
					bond+=g[iy*nx+(ix+1)%nx];
					bond+=g[(iy+1)%ny*nx+ix];
					bond+=g[(iy+1)%ny*nx+(ix+(iy&1?1:-1)+nx)%nx]; }
		e=*eps*bond; }
	return e; }


float EnergyIsingT(Lattice lt,float *eps) {
	float e;
	int ix,iy,nx,ny,bond,is2,spc,rot,*g;
	int k0,k1,k2;

	g=lt->g;
	nx=lt->nx;
	ny=lt->ny;
	spc=lt->spc;
	rot=lt->rot;

	if(rot) {
		spc=(spc+1)*4;
		e=0;
		for(iy=0;iy<ny;iy++)
			for(ix=0;ix<nx;ix++) {
				is2=g[iy*nx+ix];
				if((iy&3)==0) {
					k0=g[(iy-1+ny)%ny*nx+ix];
					k1=g[(iy+1)%ny*nx+ix];
					k2=g[(iy-1+ny)%ny*nx+(ix-1+nx)%nx]; }
				else if((iy&3)==1) {
					k0=g[(iy+1)%ny*nx+ix];
					k1=g[(iy+1)%ny*nx+(ix-1+nx)%nx];
					k2=g[(iy-1+ny)%ny*nx+ix]; }
				else if((iy&3)==2) {
					k0=g[(iy-1+ny)%ny*nx+(ix+1)%nx];
					k1=g[(iy+1)%ny*nx+ix];
					k2=g[(iy-1+ny)%ny*nx+ix]; }
				else {
					k0=g[(iy+1)%ny*nx+(ix+1)%nx];
					k1=g[(iy+1)%ny*nx+ix];
					k2=g[(iy-1+ny)%ny*nx+ix]; }
				if(iy&1) {
					e+=eps[RGTT(is2)*spc+LUPT(k0)];
					e+=eps[LDWT(is2)*spc+RGTT(k1)];
					e+=eps[LUPT(is2)*spc+LDWT(k2)]; }
				else {
					e+=eps[RGTT(is2)*spc+LDWT(k0)];
					e+=eps[LDWT(is2)*spc+LUPT(k1)];
					e+=eps[LUPT(is2)*spc+RGTT(k2)]; }}}
	else if(spc>1) {
		spc++;
		e=0;
		for(iy=0;iy<ny;iy++)
			for(ix=0;ix<nx;ix++) {
				is2=g[iy*nx+ix]*spc;
				e+=eps[is2+g[(iy-1+ny)%ny*nx+ix]];
				e+=eps[is2+g[(iy+1)%ny*nx+ix]];
				e+=eps[is2+g[(iy+(iy&1?1:-1)+ny)%ny*nx+(ix+(iy&2?1:-1)+nx)%nx]]; }}
	else {
		bond=0;
		for(iy=0;iy<ny;iy++)
			for(ix=0;ix<nx;ix++)
				if(g[iy*nx+ix]) {
					bond+=g[(iy-1+ny)%ny*nx+ix];
					bond+=g[(iy+1)%ny*nx+ix];
					bond+=g[(iy+(iy&1?1:-1)+ny)%ny*nx+(ix+(iy&2?1:-1)+nx)%nx]; }
		e=*eps*bond; }
	return e/2; }


float EnergyIsing(Lattice lt,float *eps) {
	if(lt->shape=='S') return EnergyIsingS(lt,eps);
	else if(lt->shape=='H') return EnergyIsingH(lt,eps);
	else return EnergyIsingT(lt,eps); }


int clsizeS(int *g,int nx,int ny,int ix,int iy) {
	int sum;

	if(g[iy*nx+ix]<1) return 0;
	g[iy*nx+ix]*=-1;
	sum=1;
	if(rand()&1) {
		sum+=clsizeS(g,nx,ny,(ix+1)%nx,iy);
		sum+=clsizeS(g,nx,ny,ix,(iy+1)%ny);
		sum+=clsizeS(g,nx,ny,(ix-1+nx)%nx,iy);
		sum+=clsizeS(g,nx,ny,ix,(iy-1+ny)%ny); }
	else {
		sum+=clsizeS(g,nx,ny,ix,(iy-1+ny)%ny);
		sum+=clsizeS(g,nx,ny,(ix-1+nx)%nx,iy);
		sum+=clsizeS(g,nx,ny,ix,(iy+1)%ny);
		sum+=clsizeS(g,nx,ny,(ix+1)%nx,iy); }
	return sum; }


int clsizeH(int *g,int nx,int ny,int ix,int iy) {
	int sum;

	if(g[iy*nx+ix]<1) return 0;
	g[iy*nx+ix]*=-1;
	sum=1;
	if(rand()&1) {
		sum+=clsizeH(g,nx,ny,(ix+1)%nx,iy);
		sum+=clsizeH(g,nx,ny,ix,(iy+1)%ny);
		sum+=clsizeH(g,nx,ny,(ix-1+nx)%nx,iy);
		sum+=clsizeH(g,nx,ny,ix,(iy-1+ny)%ny);
		sum+=clsizeH(g,nx,ny,(ix+(iy&1?1:-1)+nx)%nx,(iy+1)%ny);
		sum+=clsizeH(g,nx,ny,(ix+(iy&1?1:-1)+nx)%nx,(iy-1+ny)%ny); }
	else {
		sum+=clsizeH(g,nx,ny,(ix+(iy&1?1:-1)+nx)%nx,(iy-1+ny)%ny);
		sum+=clsizeH(g,nx,ny,(ix+(iy&1?1:-1)+nx)%nx,(iy+1)%ny);
		sum+=clsizeH(g,nx,ny,ix,(iy-1+ny)%ny);
		sum+=clsizeH(g,nx,ny,(ix-1+nx)%nx,iy);
		sum+=clsizeH(g,nx,ny,ix,(iy+1)%ny);
		sum+=clsizeH(g,nx,ny,(ix+1)%nx,iy); }
	return sum; }


int clsizeT(int *g,int nx,int ny,int ix,int iy) {
	int sum;

	if(g[iy*nx+ix]<1) return 0;
	g[iy*nx+ix]*=-1;
	sum=1;
	if(rand()&1) {
		sum+=clsizeT(g,nx,ny,ix,(iy+1)%ny);
		sum+=clsizeT(g,nx,ny,ix,(iy-1+ny)%ny);
		sum+=clsizeT(g,nx,ny,(ix+(iy&2?1:-1)+nx)%nx,(iy+(iy&1?1:-1)+ny)%ny); }
	else {
		sum+=clsizeT(g,nx,ny,(ix+(iy&2?1:-1)+nx)%nx,(iy+(iy&1?1:-1)+ny)%ny);
		sum+=clsizeT(g,nx,ny,ix,(iy-1+ny)%ny);
		sum+=clsizeT(g,nx,ny,ix,(iy+1)%ny); }
	return sum; }


int CtClustersIsing(Lattice lt,int *count,int ncount) {
	int i,max,ct,nx,ny,*g;
	char shape;

	g=lt->g;
	nx=lt->nx;
	ny=lt->ny;
	shape=lt->shape;

	max=0;
	for(i=0;i<ncount;i++) count[i]=0;
	for(i=0;i<nx*ny;i++)
		if(g[i]>0) {
			if(shape=='S') ct=clsizeS(g,nx,ny,i%nx,i/nx);
			else if(shape=='H') ct=clsizeH(g,nx,ny,i%nx,i/nx);
			else ct=clsizeT(g,nx,ny,i%nx,i/nx);
			if(ct<ncount) count[ct]++;
			if(ct>max) max=ct; }
	for(i=0;i<nx*ny;i++) g[i]*=-1;
	return max; }


void DisplayClustersIsing(Lattice lt) {
	int *count,max,i,nx,ny;

	nx=lt->nx;
	ny=lt->ny;
	count=(int*)calloc(nx*ny,sizeof(int));
	if(!count) return;
	max=CtClustersIsing(lt,count,nx*ny);
	printf("size\tnum\t\tconc.\n");
	for(i=1;i<=max;i++)
		if(count[i]) printf("%i\t\t%i\t\t%f\n",i,count[i],1.0*count[i]/(nx*ny));
	free(count);
	return; }

float MinDistIsing(Lattice lt,int ix,int iy,int jx,int jy) {
	float dist,distx,disty;
	int nx,ny,dx,dy;

	nx=lt->nx;
	ny=lt->ny;
	dx=jx-ix;
	dy=jy-iy;

	if(lt->shape=='S') {
		if(dx<-nx/2) dx+=nx;
		else if(dx>nx/2) dx-=nx;
		if(dy<-ny/2) dy+=ny;
		else if(dy>ny/2) dy-=ny;
		dist=sqrt(dx*dx+dy*dy); }
	else if(lt->shape=='H') {
		if((iy&1)==(jy&1)) distx=dx;
		else if(iy&1) distx=dx-0.5;
		else distx=dx+0.5;
		if(distx<-nx/2) distx+=nx;
		else if(distx>nx/2) distx-=nx;
		if(dy<-ny/2) dy+=ny;
		else if(dy>ny/2) dy-=ny;
		dist=sqrt(distx*distx+3.0/4.0*dy*dy); }
	else {
		if((iy&2)==(jy&2)) distx=sqrt(3)*dx;
		else if(iy&2) distx=sqrt(3)*(dx-0.5);
		else distx=sqrt(3)*(dx+0.5);
		if((iy&3)==(jy&3)) disty=3.0*dy/4.0;
		else if((iy&1)==(jy&1)) disty=1.5+3.0*floor(dy/4.0);
		else if((iy&1)&&(jy&3)==(iy+1&3)) disty=0.5+3.0*floor(dy/4.0);
		else if(iy&1) disty=2.0+3.0*floor(dy/4.0);
		else if((jy&3)==(iy+1&3)) disty=1.0+3.0*floor(dy/4.0);
		else disty=2.5+3.0*floor(dy/4.0);
		if(distx<-nx*sqrt(3)/2.0) distx+=nx*sqrt(3);
		else if(distx>nx*sqrt(3)/2.0) distx-=nx*sqrt(3);
		if(disty<-nx*3.0/8.0) disty+=nx*3.0/4.0;
		else if(disty>nx*3.0/8.0) disty-=nx*3.0/4.0;
		dist=sqrt(distx*distx+disty*disty); }
	return dist; }


void RadCorrFnIsing(Lattice lt,int spc1,int spc2,float *bins,int nbins,float binsize) {
	int i,j,nx,ny,obit,n1,n2;
	int *g;
	float d;
	char shape;

	nx=lt->nx;
	ny=lt->ny;
	g=lt->g;
	shape=lt->shape;
	for(i=0;i<nbins;i++) bins[i]=0;
	if(!lt->rot) obit=0;
	else obit=shape=='H'?3:2;
	for(i=0;i<nx*ny;i++)
		if(g[i]>>obit==spc1)
			for(j=i+1;j<nx*ny;j++)
				if(g[j]>>obit==spc2) {
					d=MinDistIsing(lt,i%nx,i/nx,j%nx,j/nx);
					if(d/binsize<nbins) bins[(int)(d/binsize)]+=2; }
	n1=n2=0;
	for(i=0;i<nx*ny;i++) {
		if(g[i]>>obit==spc1) n1++;
		if(g[i]>>obit==spc2) n2++; }
	if(spc1==spc2) bins[0]+=n1;
	if(!n1||!n2) return;
	if(shape=='S') d=nx*ny;
	else if(shape=='H') d=nx*ny*sqrt(3)/2.0;
	else d=nx*ny*sqrt(3)*3.0/4.0;
	for(i=0;i<nbins;i++)
		bins[i]/=n1*n2*3.1415926535*((i+1)*(i+1)-i*i)*binsize*binsize/d;
	return; }


void DisplayRCFIsing(Lattice lt,int spc1,int spc2,float binsize) {
	float *bins;
	int i,nbins,nx,ny;
	float dx,dy,d;
	char shape;
	
	shape=lt->shape;
	nx=lt->nx;
	ny=lt->ny;
	if(shape=='S') {dx=nx;dy=ny;}
	else if(shape=='H') {dx=nx;dy=sqrt(3)/2.0*ny;}
	else {dx=sqrt(3)*nx;dy=3.0/4.0*ny;}
	d=dx<dy?dx/2.0:dy/2.0;
	nbins=(int)(d/binsize);
	bins=(float*)calloc(nbins,sizeof(float));
	if(!bins) return;
	RadCorrFnIsing(lt,spc1,spc2,bins,nbins,binsize);
	printf("indx\trad\t\tg(r)\n");
	for(i=0;i<nbins;i++)
		printf("%i\t%f\t%f\n",i,binsize*(i+0.5),bins[i]);
	free(bins);
	return; }


float LatticeSiteDistH(float *xptr,float *yptr,int *ixptr,int *iyptr,float sp) {
	float x,y,spy,d;
	int xi,yi;

	x=*xptr;
	y=*yptr;
	spy=sp*0.86602540378;
	sp/=2.0;
	xi=(int)floor(x/sp);
	yi=(int)floor(y/spy);
	if((xi+yi)%2) {
		d=sqrt((x-(xi+1)*sp)*(x-(xi+1)*sp)+(y-yi*spy)*(y-yi*spy));
		if(d<sp) {*xptr=(xi+1)*sp;*yptr=yi*spy;*ixptr=(int)floor(0.5*(xi+1));*iyptr=yi;}
		else {*xptr=xi*sp;*yptr=(yi+1)*spy;d=2.0*sp-d;*ixptr=(int)floor(0.5*xi);*iyptr=yi+1;}}
	else {
		d=sqrt((x-xi*sp)*(x-xi*sp)+(y-yi*spy)*(y-yi*spy));
		if(d<sp) {*xptr=xi*sp;*yptr=yi*spy;*ixptr=(int)floor(0.5*xi);*iyptr=yi;}
		else {*xptr=(xi+1)*sp;*yptr=(yi+1)*spy;d=2.0*sp-d;*ixptr=(int)floor(0.5*(xi+1));*iyptr=yi+1;}}
	return d; }


/* ************************* OpenGL ************************ */
//#ifdef __gl_h_

#include <gl.h>
#include <glut.h>

void DisplayIsingGL(Lattice lt,float *colors) {
	int ix,iy,ic,nx,ny,spc;
	float x,y,color[4];
	int *g;

	nx=lt->nx;
	ny=lt->ny;
	spc=lt->spc;
	g=lt->g;
	if(spc>16) spc=16;
	if(nx<2||ny<2) return;

	if(lt->shape=='S'&&!lt->rot) {
		for(iy=0;iy<ny;iy++) {
			for(ix=0;ix<nx;ix++) {
				if(g[iy*nx+ix]) {
					for(ic=0;ic<4;ic++) color[ic]=colors[g[iy*nx+ix]*4+ic];
					glColor4fv(color);
					glRectf(ix-0.5,iy-0.5,ix+0.5,iy+0.5); }}}}
	else if(lt->shape=='S') {
		for(iy=0;iy<ny;iy++) {
			for(ix=0;ix<nx;ix++) {
				if(g[iy*nx+ix]) {
					for(ic=0;ic<4;ic++) color[ic]=colors[(g[iy*nx+ix]>>2)*4+ic];
					glColor4fv(color);
					glRectf(ix-0.5,iy-0.5,ix+0.5,iy+0.5); }}}}
	else if(lt->shape=='H'&&!lt->rot) {
		for(iy=0;iy<ny;iy++) {
			for(ix=0;ix<nx;ix++) {
				if(g[iy*nx+ix]) {
					x=(iy&1)?ix+0.5:ix;
					y=iy*0.5*sqrt(3);
					for(ic=0;ic<4;ic++) color[ic]=colors[g[iy*nx+ix]*4+ic];
					glColor4fv(color);
					glBegin(GL_POLYGON);
					glVertex2f(x,y-1.0/sqrt(3));
					glVertex2f(x-0.5,y-0.5/sqrt(3));
					glVertex2f(x-0.5,y+0.5/sqrt(3));
					glVertex2f(x,y+1.0/sqrt(3));
					glVertex2f(x+0.5,y+0.5/sqrt(3));
					glVertex2f(x+0.5,y-0.5/sqrt(3));
					glEnd(); }}}}
	else if(lt->shape=='H') {
		for(iy=0;iy<ny;iy++) {
			for(ix=0;ix<nx;ix++) {
				if(g[iy*nx+ix]) {
					x=(iy&1)?ix+0.5:ix;
					y=iy*0.5*sqrt(3);
					for(ic=0;ic<4;ic++) color[ic]=colors[(g[iy*nx+ix]>>3)*4+ic];
					glColor4fv(color);
					glBegin(GL_POLYGON);
					glVertex2f(x,y-1.0/sqrt(3));
					glVertex2f(x-0.5,y-0.5/sqrt(3));
					glVertex2f(x-0.5,y+0.5/sqrt(3));
					glVertex2f(x,y+1.0/sqrt(3));
					glVertex2f(x+0.5,y+0.5/sqrt(3));
					glVertex2f(x+0.5,y-0.5/sqrt(3));
					glEnd(); }}}}

	return; }


//#endif