2 D四角格子Isingモデルコード

#include 
#include 
#include 
#include 
#include

using namespace std;

ofstream out("data.dat");
const int L = 50;                //                 
const int MCS = 20000;           //        
const int MS = 10000;             //          
const double J = -1;             //J=1AFM,J=-1FM
            
double H = 0;

double InitialSpin(int s[][L])
{
	for (int i = 0; i < L; i++)
	{
		for (int j = 0; j < L; j++)
		{
			//s[i][j] = 1; 
			s[i][j]= 2 * (rand() % 2) - 1;
		}
	}
	return 0;
}
//          
//----------------------------------------

double ExchangeEnergy(int Spin[][L])
{
	int down,up,left,right;
	double e = 0;
	for (int i = 0; i < L; i++)
	{
		for (int j = 0; j < L; j++)
		{
			
			e += J*Spin[i][j]*(Spin[(i+1)%L][j]+Spin[(i-1+L)%L][j]+Spin[i][(j+1)%L]+Spin[i][(j-1+L)%L]);
		}
	}
	return e;
}
//          

//          
double Magnetization(int s[][L])
{
	double m = 0;
	for (int i = 0; i < L; i++)
	{
		for (int j = 0; j < L; j++)
		{
			m += s[i][j];
		}
	}
	return m;
}
//        
void MonteCarlo(int Spin[][L],double T)
{
	int iaccept = 0;//       
	double se = 0;
	double se2 = 0;
	double sm = 0;
	double sm2 = 0;
	double de = 0;
	double m;
	int up,down,right,left;
	int E1,E2;
	double p,r;
	for (int mcs = MCS; mcs > 0; mcs--)
	{
		for(int i=0;ip)
                 {
                 	Spin[i][j]=-Spin[i][j];
                 	iaccept++;
				 }
                
                
		       }
           
			}
		}
				if (mcs < MS)
		{
			m=0;
			 m = Magnetization(Spin) / L / L;         //      
			sm += m;         //            
			sm2 += m*m;      //            
			double e = ExchangeEnergy(Spin) / L / L - H*m;  //      
			se += e;        //            
			se2 += e*e;    //              
		}
		}
	se/=MS;
	se2 /= MS;
	sm /= MS;
	sm2 /= MS;
	double Capacity = (se2 - se*se) / T / T;//      
	double Manet = (sm2 - sm*sm) / T;   //      
	out << T << '\t'  << '\t' << se << '\t' << Capacity << '\t' << sm << '\t' <