111.txt

晶体生长蒙特卡洛模拟机制 希望和大家共享

/*	Demo version */
#include	<stdio.h>
#include	<math.h>
#include	<grwin.h>
#include	<time.h>

#define		GR
//#define		DEBUG

#define		M_SIZE	1024
#define		HORIZON 256
#define		BASE	50
#define		PSET_SIZE	2

#define		CRT_WIDE	512
#define		CRT_HIGHT	512

#define		STATE	64
#define		M_SIZE_2	(M_SIZE*M_SIZE)

#define		ulong	unsigned long

unsigned char	Crystal[HORIZON][M_SIZE][M_SIZE];	/*原子种类 */
static unsigned int	Surface[M_SIZE][M_SIZE];		/*表面高*/
static unsigned char	Config[M_SIZE][M_SIZE];			/*表面形状*/
static unsigned int	Surface_pt[M_SIZE][M_SIZE];	/*tbl相反 pt*/
static unsigned int	Surface_tblX[M_SIZE_2*2];		/*集合X*/
static unsigned int	Surface_tblY[M_SIZE_2*2];		/*集合Y*/
static unsigned int	Set_pt[STATE+1];	/* Pt到会集*/
static unsigned int	Nxij[STATE];				/*集合的元素数*/
static double		Prob[STATE][3];			/*蜕变概率*/
static double		Kxij[STATE],PS[STATE];		/*转折频率*/
static double		State_tree[STATE*2];		/*转折果断桌*/
static	int		choice_flg[STATE];		/*Kxij[k]!=0*/

static unsigned int	P_X[M_SIZE+2],P_Y[M_SIZE+2];/*周期条件*/

static	int		Base_h[M_SIZE][M_SIZE];
static	int		Color_Base[M_SIZE];
static  int             Step_N;                         /* Step No. */

static double fai_aa,fai_bb,fai_ab,fai_ab_v,k_a,k_b,alpha,ramda_a,ramda_b;	/*物理変数*/
static	double	F,FM,DMYU;			/*总吸附频率*/
static	double exp3Paa;				//平均吸附频率
static	double		Temp=300;
static 	int		sizex,sizey,SIZE2,Limit;	/*结晶大小*/
static	int		SPI_UP_x,SPI_LOW_x,SPI_UP_y,SPI_LOW_y;
static	long	C_a,C_b,E_a,E_b,kount,NextTime;
static	double	Time;
double	R_RANDOM();
unsigned int I_RANDOM();
unsigned int M_RANDOM();
unsigned char config();
static	int	Disp_mode,timing,Interval,xs1,ys1,xs2,ys2;

static long event;
static  double    INIT,INIT1;
//Gr
static int COLOR_INDEX[16];
static int O_IW,O_IH;

#define		A_atom	1
#define		B_atom	9
#define		CN	7
#define		CNP	CN*16-BASE
int	HEIKOU,TEIJO;			/*整个频率、F是否是锁着的*/
int G[16][2];		//2004/从设置A，用B原子连接了债券的数量得到
FILE	*LG;
int TCa,TCb,TEa,TEb,N_b[5],N_a[5],NB_a,NB_b,kink_C[5],kink_E[5];
static int	MODEL=1,Speed=0;
#define	Bius	4
int Total_AC,Total_BC,Pout;
main()
{	int k,c,i,j,knt,T_knt,count;
	double T_f,kk,Sigma;
	int sp[]={0,1,2},kkk;
	int limit_org;
	double ppp;
	int swt=0,swt1,Inter;
double s1,s2,r1,c1,g1,g2,b1,b2,b3;
	FILE	*fp;	char	bf[64];

	TCa=TCb=TEa=TEb=0;
	data_input();
#ifdef	GR
	graph_init();
#endif
	kkk=Disp_mode;
/*平衡状态界限被做*/
	if(TEIJO==1){limit_org=Limit;Limit=SIZE2/128;}
#ifdef	GR
	repaint();
	BOX_CLS();
#endif
		fp=fopen("Spiral.txt","w");
		fprintf(fp,"%d %d %d %d %f %f %f %f %f\n\n",
			sizex,sizey,Step_N,MODEL,fai_aa,fai_bb,fai_ab,alpha,DMYU);
		fprintf(fp,"Time Layer Ka Kb Sa Sb Ta Tb Ga Gb\n");
		fclose(fp);
	Inter=knt=0;Total_AC=Total_BC=0;Pout=0;Time=0;T_f=0.0;
while(1){
	C_a=C_b=E_a=E_b=0;
	monte();knt++;count++;
	Total_AC += C_a-E_a; Total_BC += C_b-E_b;

	if(Total_AC+Total_BC>Pout){

		TCa=Total_AC+Total_BC;
		for(i=0;i<5;i++)N_a[i]=N_b[i]=0;
		for(i=0;i<5;i++)for(j=0;j<5;j++){
			if(i+j<5){
				k=i*(11-i)/2+j;
				switch(i+j){
					case	2:	N_a[2]+=Nxij[k]+Nxij[k+16];
								N_b[2]+=Nxij[k+32]+Nxij[k+48];
								break;
					case	3:	N_a[3]+=Nxij[k]+Nxij[k+16];
								N_b[3]+=Nxij[k+32]+Nxij[k+48];
								break;
					case	4:	N_a[4]+=Nxij[k]+Nxij[k+16];
								N_b[4]+=Nxij[k+32]+Nxij[k+48];
								break;
				}
			}
		}
		printf("Real Time =%f t=%e G=%d K_A=%d K_B=%d S_A=%d S_B=%d T_A=%d T_B=%d B_A=%d B_B=%d\n",
			(float)clock()/CLOCKS_PER_SEC,Time,Pout/(sizex*sizey),
			N_a[2],N_b[2],N_a[3],N_b[3],N_a[4],N_b[4],Total_AC,Total_BC);
		fp=fopen("Spiral.txt","a");
		fprintf(fp,"%e %f %d %d %d %d %d %d %d %d\n",
			Time,(double)TCa/(sizex*sizey),
			N_a[2],N_b[2],N_a[3],N_b[3],N_a[4],N_b[4],Total_AC,Total_BC);
		fclose(fp);
		Pout += sizex*sizey;
		if(Pout>HEIKOU*sizex*sizey) break;

	}	

//	TCa+=C_a;TCb+=C_b;TEa+=E_a;TEb+=E_b;
	T_f +=F;
	FM=T_f/count;
if(TEIJO==1){		/*做法为了获得平衡状态*/
		if((E_a+E_b)*(C_a+C_b) !=0 ){
			F=F*(double)(E_a+E_b)/(C_a+C_b);
			prob_def();
			tree_set();
		}
		if(Limit<limit_org)Limit *=2;
}
#ifdef	GR

if((knt%10)==0){if(timing==2) repaint();
	if((knt%1000)==0)tree_set();

	fp=fopen("Control.dat","r");fscanf(fp,"%s",bf);	fclose(fp);
	if(bf[0]=='S') break;
	if(bf[0]=='P'){float x,y;
		if(Disp_mode==1)G_CLS();repaint();
		c=GWidle(&k, &x, &y, 0);
		fp=fopen("Control.dat","w");fprintf(fp,"Continue\n");fclose(fp);
	}
	c=INKEY();
	if( c == 'b') break;
		else if(c=='a'){kkk++;Disp_mode=sp[kkk%3];repaint();}
		else if( c=='c') c=wait_key();
}
#endif
//	if(knt>HEIKOU && HEIKOU>0) break;
}
/*
	T_knt=HORIZON;
	for(i=0;i<sizex;i++)for(j=0;j<sizey;j++){if(Surface[i][j]<T_knt)T_knt=Surface[i][j];}
	count=sizex/7;knt=count*7;
	fp=fopen("Spiral.txt","a");
	for(i=0;i<knt;i+=count){
		for(j=0;j<knt;j+=count){int ii,jj;
			C_a=C_b=0;
			for(ii=i;ii<i+count;ii++)for(jj=j;jj<j+count;jj++){
				for(k=T_knt;k>T_knt-50;k--){
					if(Crystal[k][ii][jj]=='A') C_a++; else C_b++;
				}
			}
//			fprintf(fp,"%d %d %e ",C_a,C_b,(double)C_b/(C_a+C_b));
			fprintf(fp,"%e ",(double)C_b/(C_a+C_b));
		}
		fprintf(fp,"\n");
	}
	fclose(fp);
*/
	data_output();
	G_OFF();
}
BOX_CLS(){int i;
return;
	if(Disp_mode==2){
		xs1=481;ys1=100;xs2=481;ys2=310;
	} else {
		xs1=0;ys1=0;xs2=110;ys2=0;
	}
	G_COLOR(0);i=COLOR_INDEX[0];
	GWsrect(xs1,480-ys1,xs1+100,480-(ys1+200),i);
	GWsrect(xs2,480-ys1,xs2+100,480-(ys2+100),i);
	G_COLOR(7);
	G_LINE(xs1,ys1,xs1+100,ys1);
	G_LINE(xs1+100,ys1,xs1+100,ys1+200);
	G_LINE(xs1+100,ys1+200,xs1,ys1+200);
	G_LINE(xs1,ys1+200,xs1,ys1);
	G_LINE(xs2,ys2,xs2+100,ys2);
	G_LINE(xs2+100,ys2,xs2+100,ys2+100);
	G_LINE(xs2+100,ys2+100,xs2,ys2+100);
	G_LINE(xs2,ys2+100,xs2,ys2);
}
EXIT(k)
	int k;
{
	data_output();
	exit(k);
}
double	Fx(double C,double D,double x){
	return(1.0-x-C*x*exp(D*(1.0-2.0*x))-alpha*(1.0-x));
}
double	dFx(double C,double D,double x){
	return(-1.0-C*exp(D*(1.0-2.0*x))*(1-2.0*D*x)+alpha);
}
double	newton(double C,double D){
	double	xo,xi;
	int i;
/*	for(i=0;i<=1000;i++){
		if(Fx(C,D,i*0.001)<0.0) break;
	}
	xo=(i-0.5)*0.001;
printf("i=%d\n",i);
*/
xo=1-alpha;
	for(i=0;i<10;i++){
		xi=xo-Fx(C,D,xo)/dFx(C,D,xo);
		if(fabs((xi-xo)/xo)<1.0E-15) break;
		xo=xi;
printf("xo=%E\n",xo);
	}
	return(xi);
}
data_input()
{//横向φAB纵长和φAB不同得好的那些
	FILE *fp;
	int k,ic;
	double	C,D;

	fp=fopen("ab.in","r");
	fscanf(fp,"%d",&k);
	fscanf(fp,"%d%lf",&TEIJO,&DMYU);
	fscanf(fp,"%d",&HEIKOU);
	fscanf(fp,"%d%d%d",&Disp_mode,&timing,&Interval);
	if(Interval<1) Interval=1;
	fscanf(fp,"%d%d%d",&sizex,&sizey,&Step_N);
	fscanf(fp,"%lf%lf%lf%lf",&fai_aa,&fai_bb,&fai_ab,&fai_ab_v);
	fscanf(fp,"%lf%lf",&ramda_a,&ramda_b);
	fscanf(fp,"%d%lf",&ic,&INIT);
	fscanf(fp,"%lf",&alpha);
	fscanf(fp,"%d",&MODEL);
	fscanf(fp,"%lf",&Temp);
	fscanf(fp,"%d",&Speed);
//	MODEL=1;
//是MODEL:零的情况不是螺形线滑法。是螺形线滑法的情况，Step_N酒馆煤气矢量
	fclose(fp);
	if(k!=1){
		restart();
	} else {
		if(ic==1){INIT=INIT1=200;}else{INIT *=100;INIT1=INIT;}
		if((fai_aa<fai_ab) && (fai_bb<fai_ab)){
			exp3Paa=exp(-3*fai_ab);
		} else {
			if( fabs((fai_aa+fai_bb)*0.5-fai_ab)<0.00001){
				F =1.0/(1.0+alpha/(1.0-alpha)*exp(-2.5*(fai_aa-fai_bb)));
				if(INIT<0) INIT1=(1-F)*100;
				F =1.0/(1.0+alpha/(1.0-alpha)*exp(-3.0*(fai_aa-fai_bb)));
				if(INIT<0) INIT=(1-F)*100;
//				INIT=INIT1;
				exp3Paa=(F*exp(-3.0*fai_aa)+(1.0-F)*exp(-3.0*fai_bb));
			} else {
				if(0/*fai_aa>fai_ab &&fai_bb>fai_ab*/){
					if(fai_aa>fai_bb)exp3Paa=exp(-3*fai_aa); else 
							exp3Paa=exp(-3*fai_bb);
				} else {
				if(alpha==0){exp3Paa=exp(-3*fai_aa);
					if(INIT<0.0) INIT=INIT1=0.0;
				}else if(alpha>=1.0){exp3Paa=exp(-3*fai_bb);
					if(INIT<0.0) INIT=INIT1=100.0;
				}else{
					C=alpha*exp(-3.0*(fai_aa-fai_bb));
					D=6.0*(fai_aa+fai_bb-2.0*fai_ab)*0.5;
					F=newton(C,D);if(F<1.E-16) F=0.0;
					if(INIT <0.0)	INIT=INIT1=(1.0-F)*100.0;
					exp3Paa=F*exp(-3.0*fai_aa)*exp(D*(1-F)*(1-F))
					+(1.0-F)*exp(-3.0*fai_bb)*exp(D*F*F);
				}}
			}
		}
		F=exp(DMYU)*exp3Paa;
//printf("表面%e大块%e涨潮%e\n",INIT1*0.01,INIT*0.01,F);getch();
		initialize();
	}
	LG=fopen("ab.log","w");
	fprintf(LG,"Size_x Size_y  %d %d\n",sizex,sizey);
	fprintf(LG,"# of Step      %d\n",Step_N);
	fprintf(LG,"φaa,φbb,φab %lf %lf %lf %lf\n",fai_aa,fai_bb,fai_ab,fai_ab_v);
	fprintf(LG,"λa,λb        %lf %lf\n",ramda_a,ramda_b);
	fprintf(LG,"δ             %lf\n",alpha);
	fprintf(LG,"F   DMYU       %E  %lf\n",exp3Paa,DMYU);
	fprintf(LG,"# of caluc     %d\n",HEIKOU);
	fclose(LG);
/*
#ifdef	GR
	printf("\033[2J\033[1;1H");
#endif
	printf("1 ： 新2 ： 继续 ");scanf("%d",&k);
	printf("1:平衡状态被做2:它计算以它给的F ");scanf("%d",&TEIJO);
	printf("重复频率 ?");scanf("%d",&HEIKOU);
	printf("水晶大小X - Y? ");scanf("%d%d",&sizex,&sizey);
	printf("No. Step ? "); scanf("%d",&Step_N);
	SIZE2=sizex*sizey;
	Limit=SIZE2*Interval;
if(timing==1) Limit= 100;
	printf("φaa,φbb,φab ? ");scanf("%lf%lf%lf",&fai_aa,&fai_bb,&fai_ab);
	printf("温度 ? ");scanf("%d",&Temp);
	fai_aa *=(800+273.15)/(Temp+273.15);
	fai_bb *=(800+273.15)/(Temp+273.15);
	fai_ab *=(800+273.15)/(Temp+273.15);
	printf("Diffusion length for A atom and B atom? ");
		scanf("%lf%lf",&ramda_a,&ramda_b);
	printf("早密度(%)更多眼睛？");scanf("%lf",&INIT1);
	INIT=INIT1;
*/
}
data_output()
{
	int	z,y,x,f;
	FILE	*fp;

	RANDOM_CLOSE();

	fp=fopen("ab.dat","wb");
	fprintf(fp,"%d\n",Step_N);
	fprintf(fp,"%d %d\n",sizex,sizey);
	fprintf(fp,"%E %E %E\n",fai_aa,fai_bb,fai_ab);
	fprintf(fp,"%E %E\n",ramda_a,ramda_b);
	fprintf(fp,"%E\n",alpha);
if(TEIJO!=1) FM=FM/exp(DMYU);
	fprintf(fp,"%E ",FM);
	fprintf(fp,"%d\n",MODEL);
	for(z=0;z<HORIZON;z++){ f=0;
		for(x=0;x<sizex;x++){
			fwrite(Crystal[z][x],1,sizey,fp);
			if( f==0 ){
				for(y=0;y<sizey;y++)
				  if(Crystal[z][x][y]!=0){f=1;break;}
			}
		}
		if( f==0 ) break;
	}
	fclose(fp);
}
restart(){
	double br,Tfai_aa,Tfai_bb,Tfai_ab,Tramda_a,Tramda_b,Talpha;
	int	z,y,x,f,z_min;
	FILE	*fp;
	fp=fopen("ab.dat","rb");
	fscanf(fp,"%d",&Step_N);
	fscanf(fp,"%d%d",&sizex,&sizey);
	fscanf(fp,"%lf%lf%lf",&Tfai_aa,&Tfai_bb,&Tfai_ab);
	fscanf(fp,"%lf%lf",&Tramda_a,&Tramda_b);
	fscanf(fp,"%lf",&Talpha);
	fscanf(fp,"%lf",&F);F=F*exp(DMYU);
	fscanf(fp,"%d\n",&MODEL);
	f=1;z=0;
	z_min=0;
	while(f==1){ f=0;
		for(x=0;x<sizex;x++){
			fread(Crystal[z][x],1,sizey,fp);
				for(y=0;y<sizey;y++)
				  if(Crystal[z][x][y]!=0){f=1;
					Surface[x][y]=z;
				  } else if(z_min==0) z_min=z;
		}
		z++;
	}
	if(z_min>BASE){
		for(f=0;f<z-z_min+BASE;f++)
			for(x=0;x<sizex;x++)for(y=0;y<sizey;y++){
				Crystal[f][x][y]=Crystal[f+z_min-BASE][x][y];
			}
		for(x=0;x<sizex;x++)for(y=0;y<sizey;y++)Surface[x][y] -=z_min-BASE;
		z=z-z_min+BASE;
	}
	
	for(;z<HORIZON;z++){
		for(x=0;x<sizex;x++)for(y=0;y<sizey;y++){
			Crystal[z][x][y]=0;
		}
	}
	fclose(fp);
	G_SET();
	SIZE2=sizex*sizey;
	Limit=SIZE2*Interval;
if(timing==1) Limit= 100;
	RANDOM_OPEN();
	perio_set();
	prob_def();
	config_set();
	nxij_set();
	tree_set();
	table_set();
	kount=0;
	Time=0.0;
}
G_SET(){
	int i,j,k;// i ：连接了到A原子和横向债券的数字&j ：B原子
	for(i=0;i<5;i++)for(j=0;j<5;j++){
		if(i+j<5){
			k=i*(11-i)/2+j;
			G[k][0]=i;
			G[k][1]=j;
		}
	}
}
initialize()
{
	G_SET();
	SIZE2=sizex*sizey;
	Limit=SIZE2*Interval;
if(timing==1) Limit=Limit/=sizex;
	RANDOM_OPEN();
	perio_set();
	prob_def();
	surface_init();
	config_set();
	nxij_set();
	tree_set();
	table_set();
	kount=0;
	Time=0.0;
}
table_set()
{
	int	i,j,k,n,sz_share,sz;
	double	ritu;
	
	n=0;
	sz_share=M_SIZE_2/60/2;
	sz=M_SIZE_2*2-sz_share*60-24;	/*	剩下的領域	*/
	ritu = (double)sz / (double)(sizex*sizey);
	for(k=0;k<STATE;k++){
		Set_pt[k]=n;
		if( Kxij[k] != 0.0 ) n+= sz_share+Nxij[k]*ritu-1.0;
		n++;
	}
	for(k=0;k<STATE;k++) Nxij[k]=0;
	for(i=0;i<sizex;i++) for(j=0;j<sizey;j++){
		k=Config[i][j];
		n=Surface_pt[i][j]=Set_pt[k]+Nxij[k];
		Surface_tblX[n]=i;
		Surface_tblY[n]=j;
		Nxij[k]++;
	}
}
tree_set()
{
	int	i,k;
	for(i=0;i<STATE;i++) {
		State_tree[i+STATE]=Nxij[i]*Kxij[i];
	}
	for(k=32;k>0;k/=2){
		for(i=k;i<2*k;i++){
			State_tree[i]=State_tree[2*i]+State_tree[2*i+1];
		}
	}
}
nxij_set()
{
	int	i,j,k;

	for(k=0;k<STATE;k++) Nxij[k]=0;
	for(i=0;i<sizex;i++) for(j=0;j<sizey;j++){
		Nxij[Config[i][j]]++;
	}
}
config_set()
{//check check
	int	i,j,k;