mdla.cpp

磁性DLA模拟分形生长

#include <iostream>
#include <math.h>
#include <fstream>
#define rand (float)rand()/RAND_MAX
#define PI 3.1415926
using namespace std;

float J,H,Rmax=0.0;
int N,Nup,Ndown;                                //J是相互作用，H是外场大小 ，N为总粒子个数 ,Rmax用来记录结构的最外围半径 
int** a;                                         //a用来存放粒子的位置信息 ,Up和Down分别表示自旋朝上和朝下的个数 

int** IntialArray()                                                   
{
	int** p;
    p = new int*[200];
	int i,j;
    for(i = 0; i < 200; i ++)
	{
		p[i] = new int[200];
		for(j = 0; j <200; j ++)
		{
			p[i][j] = 0;
		}
	}
    return p;
}

void Aggradation()
{
     int i,j,spin;                             //i,j是随即游走中的粒子坐标 ,spin=1,up;spin=-1,down 
     float temp,dpotential;                        //dpotential是势场的变化 
     for(int n=1;n<N;)
     {
          i=100+int(100*cos(2*PI*rand));           //在边界处随机产生一个粒子 
          j=100+int(100*sin(2*PI*rand));
          for(;;)
          {
               temp=rand;
               if(temp<=0.25)
                    i++;
               if(temp>0.25&&temp<=0.5)
                    j++;
               if(temp>0.5&&temp<=0.75)
                    i--;
               if(temp>0.75)
                    j--;
               if(abs(i-100)<(Rmax+10)&&abs(j-100)<(Rmax+10))                  //判断周围是否有粘附的粒子，并根据势场的变化决定是否粘上 
                    if((a[i+1][j]!=0)||(a[i][j+1]!=0)||(a[i-1][j]!=0)||(a[i][j-1]!=0))
                    {
                         spin=(rand)<0.5?1:-1;
                         dpotential=float(-spin*(a[i+1][j]+a[i][j+1]+a[i-1][j]+a[i][j-1])*J/2-spin*H);
                         if(dpotential<=0)
                              a[i][j]=spin;
                         else if((rand)<exp(-dpotential))
                              a[i][j]=spin;
                         else continue;
                         if(spin>0)
                              Nup++;
                         else Ndown++;
                         if(sqrt((i-100)*(i-100)+(j-100)*(j-100))>Rmax)
                              Rmax=sqrt((i-100)*(i-100)+(j-100)*(j-100));
                         n++;
                         break;
                    }
               if((i-100)*(i-100)+(j-100)*(j-100)>4*100*100)                      //当试验粒子随机游走得太远了，舍弃，重新产生一个 
                    break;
          }
     }
     ofstream fout("1.txt");
     fout<<"x"<<'\t'<<"y"<<'\t'<<"spin"<<endl;
     for(i=0;i<200;i++)
          for(j=0;j<200;j++)
               fout<<i<<'\t'<<j<<'\t'<<a[i][j]<<endl;
     fout.close();
}

void Cal_Dimension()
{
      float scale,Dimension,numerator=0,denominator=0;
      int i,j,k; 
      scale=0.5;
      float b[10][2]={0};                                                      //用于存放计算维数所需数据 
      cout<<"R"<<"   "<<"N"<<endl;
      ofstream fout("d_1.txt");
      for(k=0;k<10;k++)
      {
          for(i=(int)(100-Rmax*scale*(k+1)/10);i<=(int)(100+Rmax*scale*(k+1)/10);i++)
          {
               for(j=int(100-Rmax*scale*(k+1)/10);j<=100+int(Rmax*scale*(k+1)/10);j++)
               {
                    if ((a[i][j]!=0)&&(((i-100)*(i-100)+(j-100)*(j-100))<((Rmax*scale*(k+1)/10)*(Rmax*scale*(k+1)/10))))
                        b[k][1]++;
               }
          }
          b[k][0]=int(Rmax*scale*(k+1)/10)+1;
          fout<<b[k][0]<<"  "<<b[k][1]<<endl;
      }
      for(i=0;i<10;i++)                                                         //利用最小二乘法求分形结构的维数 
          for(j=0;j<10;j++)
          {
               numerator-=log(b[i][0])*log(b[j][1]);
               denominator-=log(b[i][0])*log(b[j][0]);
               if(i==j)
               {
                     numerator+=10*log(b[i][0])*log(b[i][1]);
                     denominator+=10*log(b[i][0])*log(b[i][0]);
               }
          }
      Dimension=numerator/denominator;
      fout<<"The Dimension is:"<<Dimension<<endl; 
      cout<<"The Dimension is:"<<Dimension<<endl; 
      cout<<"自旋朝上的个数为："<<Nup<<"自旋朝下的个数为："<<Ndown<<endl;        
      fout<<"自旋朝上的个数为："<<Nup<<"自旋朝下的个数为："<<Ndown<<endl;
      fout.close();
}

main()
{
     int i,j;
     cout<<"请输入需要的总粒子个数"<<endl;
     cin>>N;
     cout<<"请输入相互作用因子"<<endl;
     cin>>J;
     cout<<"请输入外场的大小"<<endl;
     cin>>H;
     a=IntialArray();
     a[100][100]=1;
     Nup=1;
     Ndown=0; 
     Aggradation();
     Cal_Dimension(); 
     system("pause");
}