有一个程序利用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;
}

连蒙带猜，问题很多啦，只说几个主要的
1. possibility中，你求的是 e^(-con/T)，可是要求中写的是 e^(-deltaE*con/T)
2. tempAngle=rand()%1000; 中tempAngle的单位是什么？角度还是弧度？cos函数的参数应当是弧度。%1000反正很奇怪，既不是360的整数倍，又不是PI的整数倍。
顺便问一下，从#include "stdafx.h"看出你用的是VC，但不知道你用的是哪个版本的VC，是VC6，还是VC2017？因为我怕我写的标准C++代码，你用VC等不标准的编译器编译失败。

能不能在40半径内出现你老师说的“突然下降到0附近”，与温度及con相关吧，结果还与一开始粒子所处位置相关吧
我是一窍不通，只根据你的要求写

#define _USE_MATH_DEFINES
#include <cmath>
#include <random>
#include <cstdlib>
#include <ctime>
#include <iostream>

double dist( double x, double y )
{
    return hypot( x, y ); // 如果你的VC不认识 hypot 就用 _hypot
}

double random( double range )
{
    static bool bfirst = true;
    if( bfirst )
    {
        bfirst = false;
        srand( unsigned(time(NULL)) );
    }
    return rand()/(1.0+RAND_MAX) * range;
}

double energy( double x, double y )
{
    double d = dist( x, y );
    return -0.5*cos(0.5*M_PI*d)-cos(1.5*M_PI*d); // 假如你的公式中用的是角度，那么这里得先乘以M_PI/180
}

double possibility( double x_from, double y_from, double x_to, double y_to, double con, double temper )
{
    double e_from = energy( x_from, y_from );
    double e_to = energy( x_to, y_to );
    return e_from>e_to ? 1 : exp((e_from-e_to)*con/temper);
}

double calculate( double temper )
{
    const int con = 200;
    const double radius = 40; // 势阱半径40
    const unsigned step = 10000; // 模拟一万次
    unsigned count = 0; // 基态数量

    double x, y; // 当前位置
    {
        double a = random( 2*M_PI ); // 随机角度
        double r = random( radius ); // 随机半径
        x = r * cos( a );
        y = r * sin( a );
    }
    for( unsigned i=0; i!=step; ++i )
    {
        double a = random( 2*M_PI ); // 随机角度
        double x_new = x + cos( a ); // 新位置
        double y_new = y + sin( a );

        if( dist(x_new,y_new) < radius ) // 不越界
        {
            double p = possibility( x,y, x_new,y_new, con, temper ); // 从(x,y)到(x_new,y_new)的跳跃几率
            if( random(1) < p ) // 跳
            {
                x = x_new;
                y = y_new;
            }
        }

        if( dist(x,y) < 14 )
            ++count;
    }

    return count*1.0/step;
}

using namespace std;

int main( void )
{
    for( int temper=100; temper<1000; temper+=10 )
    {
        cout << "Temperature: " << temper
             << " Possibility: " << calculate(temper) << '\n';
    }
}

原来这里有错！ possobility里没写了deltaE。
非常感谢！
然后tempAngle的单位是radian。啊！原来这里也不对。受到的是随机整数。
我用的是vs2010里的c++控制台程序。保存一次后用vs2017打开的。 因为vs2017里没有c++。
我晚上看代码了。
很头疼的问题找到了问题点了。
感谢大佬！

如果要我说的话，
E=-0.5cos(0.5pi*d)-cos(1.5*pi*d)
就不对，这是一个周期震荡函数，震荡周期大约在 4 左右。而理论上它应该是d越大则E越大的函数。
因为这个E公式不对，所以接下来的几率也就不对，也就没有粒子趋向于圆心的行为，导致最终结果完全是随机的。

“vs2010里的c++控制台程序”
------ VS 包含 VC++、、VC# 等等，所以不要说你用的是VS，M$说过“VS不是编译器、VC才是”

“vs2017里没有c++”
------ 是没有C++，但有VC++2017

E函数长得这样。
温度上升时，跃迁到更高的能量的地方的概率接近于1，粒子的行为越来越接近随机行走。
所以处于基态的统计概率随温度上升，逐渐下降（过某一个点突然下降），下降到某个值。
（基态的面积除以势阱的面积，比这个值稍微高一丢丢。）
大概是这么个思路。

学到了您的random函数用法。
代码主要我自己修改了。结果虽然不能说完美，但是已经很不错了。
谢谢！我要结贴。 （以后的日子。。。）