求大佬！ MonteCarlo算法
有一个程序利用MonteCarlo行走方式，模拟势阱里的粒子二维势阱的半径为40
粒子的行走方式：
1.每一次只走一步
2.每一次移动距离为固定值s
3.每一次行走的方向随机
4.行走概率：如果由行走所带来的能量差为 deltaE
则这种改变的发生的概率p为：
e^(-deltaE*con/T)) ，deltaE>=0
1 , deltaE<0
这里E是对距离（从原点到粒子）的函数
E=-0.5cos(0.5pi*d)-cos(1.5*pi*d) ,d:距离，pi:圆周率
然后设置步骤数（1000），统计这个粒子处于基态的概率（基态是从原点的距离0-14的区域）
最后的目的是得出不同温度下的概率。（对每一个温度做统计计算）
这个问题我大概就结了一个星期。（我太菜了）
写出了代码，但是结果不太理想。 有大佬了解这样的问题吗？ 我的代码哪里有问题？
老师说，正确的结果是 概率随温度慢慢下降然后过某一个点，突然下降到0附近（发生相变）。
我的代码是这样的。不到100行。
我拿出我的全部积分奖励大佬。
#include "stdafx.h"
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
using namespace std;
double dist(double x,double y)
{
double result;
result=sqrt(x*x+y*y);
return result;
}
double energy(double x,double y)
{
double result;
result=-0.5*cos(0.5*3.141592*dist(x,y)-cos(1.5*3.141592*dist(x,y)));
return result;
}
double possibility(double con,double temper)
{
return pow(2.71828,-con/temper);
}
double calculate (double temper)
{
int i,step,con=200; //概率比例系数设为40
double x[10000],y[10000],tempAngle,tempDist,tempX,tempY,p,count=0;
step=1000; //统计一万步的结果
srand((unsigned)time(NULL));
tempAngle=rand()%1000;
tempDist=rand()%40; //设势阱的半径为步伐的40倍
x[0]=tempDist*cos((double)tempAngle);
y[0]=tempDist*sin((double)tempAngle);
for (i=1;i<step;i++)
{
tempAngle=rand()%1000;
tempX=x[i-1]+cos(tempAngle);
tempY=y[i-1]+sin(tempAngle);
if (dist(tempX,tempY)>=40) //如果超过了边界
{
x[i]=x[i-1];
y[i]=y[i-1];
}
else
{
if (energy(tempX,tempY)<energy(x[i-1],y[i-1])) //如果预测位置的能量比原来位置的能量小
{
x[i]=tempX;
y[i]=tempY;
}
else
{
p=possibility(con,temper);
p*=1000;
if ((rand()%1000)<p) //以概率p
{
x[i]=tempX;
y[i]=tempY;
}
else
{
x[i]=x[i-1];
y[i]=y[i-1];
}
}
}
}
for (i=0;i<step;i++)
{
if (dist(x[i],y[i])<14)
count++;
}
return count/(double)step;
}
int main()
{
double temper;
int i;
for (temper=100;temper<1000;temper+=10)
{
cout << "Temperature:" << temper << ' ' << "Possibility"
<< calculate(temper) << endl;
}
scanf("%d",&i);
return 0;
}
