111.txt

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

do{
	rnd=R_RANDOM()*State_tree[1];
	k=1;
	while( k < 64 ){
		k*=2;
		if( rnd > State_tree[k] ) {rnd-=State_tree[k];k++;}
	}
	k-=64;
	if( choice_flg[k] == 0 || Nxij[k]<=0) tree_set();
} while( choice_flg[k] == 0 || Nxij[k]<=0 );
	return(k);
}
unsigned char config(a,b,c,i,j)
	int	*a,*b,*c,i,j;
{
	int	I,J,hight,hbius,z;

	I=J=0;
	hight=Surface[i][j];
	if( Crystal[hight-1][i][j]=='A' ) 		*c=0; else *c=16;// 0:在A之下; 16 ：在B之下
	z=Base_h[i][j]-Base_h[i][P_Y[j+2]];
	if(z/3==0){hbius=hight;}
	else if(z>0){hbius=hight+(z+1)/4;} else {hbius=hight+(z-1)/4;}
	switch( Crystal[hbius][i][P_Y[j+2]] ){
		case 'A'	: I++;break;
		case 'B'	: J++;break;
	}
	z=Base_h[i][j]-Base_h[i][P_Y[j]];
	if(z/3==0){hbius=hight;}
	else if(z>0){hbius=hight+(z+1)/4;} else {hbius=hight+(z-1)/4;}
	switch( Crystal[hbius][i][P_Y[j]] ){
		case 'A'	: I++;break;
		case 'B'	: J++;break;
	}
	z=Base_h[i][j]-Base_h[P_X[i+2]][j];
	if(z/3==0){hbius=hight;}
	else if(z>0){hbius=hight+(z+1)/4;} else {hbius=hight+(z-1)/4;}
	switch( Crystal[hbius][P_X[i+2]][j] ){
		case 'A'	: I++;break;
		case 'B'	: J++;break;
	}
	z=Base_h[i][j]-Base_h[P_X[i]][j];
	if(z/3==0){hbius=hight;}
	else if(z>0){hbius=hight+(z+1)/4;} else {hbius=hight+(z-1)/4;}
	switch( Crystal[hbius][P_X[i]][j] ){
		case 'A'	: I++;break;
		case 'B'	: J++;break;
	}
	*a=I;
	*b=J;
	return( Crystal[hight][i][j] );
}
#ifdef	GR
wait_key()
{
	int c;
	while( (c=INKEY()) == 255 );
	return(c);
}
#endif
/*在惯例之下为图表*/
#ifdef	GR
static int POINT_SIZE;
#define LX 6
#define LY 4
#define LZ 4
#define Hd 2
static int Lx=6;
static int Ly=4;
static int Lz=4;
static int Lxx,Lyy,Lzz,Lzs,Lshift;
static int Show_sizex;
static int Screen[1280];
static char OLD_COLOR[M_SIZE][M_SIZE];

PSETWI(int n)
{int IU,IV;
	if( n<2 ) n=1;
    POINT_SIZE=n;
	G_COLOR_INDEX(COLOR_INDEX);
// GWldcsiz((float)(n), (float)(n), &IU, &IV);
//printf("Size= %d %d %d\n",n,IU,IV);
//GWsetpen(-1,-1,IU, -1);
}
PSETW(int x,int y,int c)
{ int i;
if(Disp_mode !=2 )return;
    x *=POINT_SIZE;
    y *=POINT_SIZE;
	GWsrect((float)x,(float)(CRT_HIGHT-y),(float)(x+POINT_SIZE),
			(float)(CRT_HIGHT-(y+POINT_SIZE)),COLOR_INDEX[c%16]);
}

PSET_BI(int n)
{
int i;
if(sizex>=sizey)i=sizex;else i=sizey;
PSETWI(CRT_HIGHT/i);
for(i=0;i<1280;Screen[i]=4000,i++);
	if( sizex > 31 ) Show_sizex=31; else Show_sizex=sizex-1;
	Lx=LX*n; Ly=LY*n; Lz=LZ*n;
	Lxx=Lx-1; Lyy=Ly-1; Lzz=Lz-1;
	Lzs=Lz*Show_sizex;Lshift=CRT_HIGHT-1-Show_sizex*Lz-3*Lz;
}
int kabeSW=1;
repaint()
{
	int j,k,h,KLim;
	int c;
	KLim=sizey;
	if(sizey>=CRT_HIGHT){KLim=CRT_HIGHT;}
//GWmode(1,1);
//GWsetogn(0,-1);
if( Disp_mode == 2){	
	for(j=0;j<sizex;j++) for(k=0;k<KLim;k++) {
		h=Surface[j][k];
		if( Crystal[h][j][k] == 'A' ){
				c=(h+CNP+Color_Base[j])%CN + A_atom;
        } else {
				c=(h+CNP+Color_Base[j])%CN + B_atom;
        }
		if(c!=OLD_COLOR[j][k]){PSETW(j,k,c);OLD_COLOR[j][k]=c;}
	}
}
if( Disp_mode == 1 ){
	for(j=0;j<=Show_sizex;j++) for(k=0;k<=Show_sizex;k++) {
		h=Surface[j][k];
		if( Crystal[h][j][k] == 'A' ){
			PSET_BN(j,k,h,A_atom);
		} else {
			PSET_BN(j,k,h,B_atom);
		}
	}
	kabe_nuri();
}
//GWflush(-1);
}
kabe_nuri()
{
	int	i,j;
	for(i=0;i<=Show_sizex;i++){
		for(j=BASE-3;j<Surface[Show_sizex][i]+Base_h[Show_sizex][i];j++)
			{kabex(Show_sizex,i,j);}
		for(j=BASE-3;j<Surface[i][Show_sizex]+Base_h[i][Show_sizex];j++)
			{kabey(i,Show_sizex,j);}
	}
}
float PL[4][2];
kabex(x,y,z)
{
	int h,xx,yy,i,xs,ys,hh;

	h=z-BASE;
	if( Crystal[z-Base_h[x][y]][x][y] == 'A' ) hh= (z+CNP)%CN + A_atom;
				else  hh= (z+CNP)%CN + B_atom;
	xx=x*Lx+Lzs-y*Lz;
	yy=y*Lz-h*Ly+Lshift;
	G_COLOR(hh);hh=COLOR_INDEX[hh];
	G_COLOR(0);
	GWsetbrs(hh,1,0);
	PL[0][0]=xx+Lx;PL[0][1]=CRT_HIGHT-yy;
	PL[1][0]=xx+Lx;PL[1][1]=CRT_HIGHT-(yy-Ly);
	PL[2][0]=xx+Lx+Lz;PL[2][1]=CRT_HIGHT-(yy-Ly-Lz);
	PL[3][0]=xx+Lx+Lz;PL[3][1]=CRT_HIGHT-(yy-Lz);
	GWpolygon(PL,4,0);	
}
kabey(x,y,z)
{
	int h,xx,yy,i,xs,ys,hh;

	h=z-BASE;
	if( Crystal[z-Base_h[x][y]][x][y] == 'A' ) hh= (z+CNP)%CN + A_atom;
				else  hh= (z+CNP)%CN + B_atom;
	xx=x*Lx+Lzs-y*Lz;
	yy=y*Lz-h*Ly+Lshift;
	G_COLOR(hh);hh=COLOR_INDEX[hh];
	G_COLOR(0);
	GWsetbrs(hh,1,1);
	PL[0][0]=xx;PL[0][1]=CRT_HIGHT-yy;
	PL[1][0]=xx+Lx;PL[1][1]=CRT_HIGHT-yy;
	PL[2][0]=xx+Lx;PL[2][1]=CRT_HIGHT-(yy-Ly);
	PL[3][0]=xx;PL[3][1]=CRT_HIGHT-(yy-Ly);
	GWpolygon(PL,4,0);	
}
PSET_BN(int x,int y,int z,int c)
{
	int h,hh,xx,yy,i,xs,ys;
	if(Disp_mode !=1 ) return;
	if(x<0 || x>Show_sizex || y<0 || y>Show_sizex) return;
	h=z-BASE; hh=(z+CNP)%CN+c;
	xx=x*Lx+Lzs-y*Lz;
	yy=y*Lz-h*Ly+Lshift;
	G_COLOR(hh);hh=COLOR_INDEX[hh];
	G_COLOR(0);
	GWsetbrs(hh,1,1);
	PL[0][0]=xx;PL[0][1]=CRT_HIGHT-yy;
	PL[1][0]=xx+Lx;PL[1][1]=CRT_HIGHT-yy;
	PL[2][0]=xx+Lx;PL[2][1]=CRT_HIGHT-(yy-Ly);
	PL[3][0]=xx;PL[3][1]=CRT_HIGHT-(yy-Ly);
	GWpolygon(PL,4,0);
	PL[0][0]=xx+Lx;PL[0][1]=CRT_HIGHT-yy;
	PL[1][0]=xx+Lx;PL[1][1]=CRT_HIGHT-(yy-Ly);
	PL[2][0]=xx+Lx+Lz;PL[2][1]=CRT_HIGHT-(yy-Ly-Lz);
	PL[3][0]=xx+Lx+Lz;PL[3][1]=CRT_HIGHT-(yy-Lz);
	GWpolygon(PL,4,0);
	PL[0][0]=xx;PL[0][1]=CRT_HIGHT-(yy-Ly);
	PL[1][0]=xx+Lx;PL[1][1]=CRT_HIGHT-(yy-Ly);
	PL[2][0]=xx+Lx+Lz;PL[2][1]=CRT_HIGHT-(yy-Ly-Lz);
	PL[3][0]=xx+Lz;PL[3][1]=CRT_HIGHT-(yy-Ly-Lz);
	GWpolygon(PL,4,0);
}
PSET_H_check(int x,int y,int L){
	int xs,ys,SW,i;
	SW=0;
	xs=x;ys=y;
	for(i=0;i<L;i++){
		if(Screen[xs]>=ys){SW=1;break;}
		xs++;ys--;
	}
	if(SW==0) return;
	x=xs;y=ys;
	for(;i<L;i++){
		xs++;ys--;
		if(Screen[xs]<ys){xs--;ys++;break;}
	}
	G_LINE(x,y,xs,ys);
}
PSET_BNN(int x,int y,int z,int c)
{
	int h,hh,xx,yy,i,xs,ys,k;
	h=z-BASE; hh=(z+CNP)%CN+c;
	xx=x*Lx+Lzs-y*Lz;
	yy=y*Lz-h*Ly+Lshift;xs=hh;
	G_COLOR(hh);hh=COLOR_INDEX[hh];
	GWsetbrs(hh,1,1);
	k=0;
		if(Screen[xx+2]>=yy-Ly-1)
			{PL[k][0]=xx;PL[k][1]=CRT_HIGHT-(yy-Ly);k++;}
		if(Screen[xx+Lx]>=yy-Ly)
			{PL[k][0]=xx+Lx;PL[k][1]=CRT_HIGHT-(yy-Ly);k++;}
		else
			{for(i=xx+Lx-1;i>xx;i--){
					if(Screen[i]>=yy-Ly){
						PL[k][0]=i;PL[k][1]=CRT_HIGHT-(yy-Ly);k++;
						break;
					}
				}
			}
		if(Screen[xx+Lx+Lz-2]>=yy-Ly-Lz+1)
			{PL[k][0]=xx+Lx+Lz;PL[k][1]=CRT_HIGHT-(yy-Ly-Lz);k++;}
		else
			{for(i=xx+Lx+Lz-2;i>xx+Lz;i--){
					if(Screen[i]>=yy-Lz-Ly+1){
						PL[k][0]=i;PL[k][1]=CRT_HIGHT-(yy-Ly-Lz);k++;
						break;
					}
				}
			}
		if(Screen[xx+Lz]>=yy-Ly-Lz+1)
			{PL[k][0]=xx+Lz;PL[k][1]=CRT_HIGHT-(yy-Ly-Lz);k++;}
		else
			{for(i=xx+Lz+1;i<xx+Lx+Lz;i++){
					if(Screen[i]>=yy-Lz-Ly+1){
						PL[k][0]=i;PL[k][1]=CRT_HIGHT-(yy-Ly-Lz);k++;
						break;
					}
				}
			}
	if(k>2)GWpolygon(PL,k,0);//上面
	k=0;
		if(Screen[xx+Lx]>=yy)
			{PL[k][0]=xx+Lx;PL[k][1]=CRT_HIGHT-yy;k++;}
		if(Screen[xx+Lx+1]>=yy-Ly)
			{PL[k][0]=xx+Lx;PL[k][1]=CRT_HIGHT-(yy-Ly);k++;}
		if(Screen[xx+Lx+Lz-1]>=yy-Ly-Lz+2)
			{PL[k][0]=xx+Lx+Lz;PL[k][1]=CRT_HIGHT-(yy-Ly-Lz);k++;}
		if(Screen[xx+Lx+Lz-1]>=yy-Lz)
			{PL[k][0]=xx+Lx+Lz;PL[k][1]=CRT_HIGHT-(yy-Lz);k++;}
		if(k>2) GWpolygon(PL,k,0);//右側面
	k=0;
		if(Screen[xx]>=yy)//{
			{PL[k][0]=xx;PL[k][1]=CRT_HIGHT-yy;k++;}
		if(Screen[xx+Lx-1]>=yy)
			{PL[k][0]=xx+Lx;PL[k][1]=CRT_HIGHT-yy;k++;}
		if(Screen[xx+Lx-1]>=yy-Ly+1)
			{PL[k][0]=xx+Lx;PL[k][1]=CRT_HIGHT-(yy-Ly);k++;}
		if(Screen[xx+1]>=yy-Ly+1)
			{PL[k][0]=xx;PL[k][1]=CRT_HIGHT-(yy-Ly);k++;}
	if(k>2) GWpolygon(PL,k,0);//前面
	G_H_H_line(xx+Lz,yy-Ly-Lz,0,Lx);
	G_H_H_line(xx,yy,0,Lx);
	G_H_H_line(xx,yy-Ly,0,Lx);
	G_H_V_line(xx,yy-Ly,0,Ly+1);
	G_H_V_line(xx+Lx,yy-Ly,0,Ly+1);
	G_H_V_line(xx+Lx+Lz,yy-Ly-Lz,0,Ly+1);
	G_COLOR(0);
	PSET_H_check(xx+Lx,yy-Ly,Ly);
	PSET_H_check(xx+Lx,yy,Ly);
	PSET_H_check(xx,yy-Ly,Ly);
}
PSET_BH(int x,int y,int z,int c)
{
	int h,hh,xx,yy,i,xs,ys;
	if(Disp_mode != 1) return;
	if(x<0 || x>Show_sizex || y<0 || y>Show_sizex) return;
	Make_Screen(x,y);
	PSET_BNN(x,y,z,c);
}
Make_Screen(int x,int y)
{
int i,j,y1,x2,xs,xe;
	y1=y+Hd;if(y1>Show_sizex) y1=Show_sizex;
	x2=x+Hd;if(x2>Show_sizex) x2=Show_sizex;
	xs=x*Lx+Lzs-y1*Lz;
	xe=(x2+1)*Lx+Lzs-(y-1)*Lz;
	for(i=xs;i<=xe;i++) Screen[i]=4000;
	for(i=x+1;i<=x2;i++){Mask_Screen(i,y);}
	for(j=y+1;j<=y1;j++) for(i=x;i<=x2;i++){Mask_Screen(i,j);}
}
Mask_Screen(int x,int y)
{
int	xx,yy,h,i;
	if( x<0 || x>Show_sizex || y<0 || y>Show_sizex ) return;
	h=Surface[x][y]-BASE+Base_h[x][y];
	xx=x*Lx+Lzs-y*Lz;
	yy=y*Lz-(h+1)*Ly+Lshift+1;
	for(i=xx;i<=xx+Lz;i++){
		yy--;
		if( yy<Screen[i] ) Screen[i] =yy;
	}	
	for(;i<=xx+Lzz+Lx;i++){
		if( yy<Screen[i] ) Screen[i] =yy;
	}
}
PSET_BE(int x,int y,int z,int c)
{
	int i,j,a,b,zz;
	if(Disp_mode !=1 ) return;
	if(x<0 || x>Show_sizex || y<0 || y>Show_sizex) return;
	if( y<1 || x<1 ) {PSET_cl(x,y,z+1);}
	for(i=x-Hd;i<=x;i++) if( i>=0 && i<=Show_sizex ) {
 		for(j=y-Hd;j<=y;j++){
		if( j>=0 && j<=Show_sizex ) {
		if( config(&a,&b,&zz,i,j) == 'A' ) c=A_atom; else c=B_atom;
//getch();
			PSET_BH(i,j,Surface[i][j]+Base_h[i][j],c);
		}
	}}
}
PSET_cl(int x,int y,int z)
{
int h,hh,xx,yy,i,xs,ys;
Make_Screen(x,y);
h=z-BASE;
	xx=x*Lx+Lzs-y*Lz;
	yy=y*Lz-h*Ly+Lshift;
	G_COLOR(0);hh=COLOR_INDEX[0];
	GWsetbrs(hh,1,1);
/*
	PL[0][0]=xx;PL[0][1]=CRT_HIGHT-yy;
	PL[1][0]=xx+Lx+Lz;PL[1][1]=CRT_HIGHT-yy;
	PL[2][0]=xx+Lx+Lz;PL[2][1]=CRT_HIGHT-(yy-Ly-Lz);
	PL[3][0]=xx;PL[3][1]=CRT_HIGHT-(yy-Ly-Lz);
	GWpolygon(PL,4,0);
	
*/
	PL[0][0]=xx;PL[0][1]=CRT_HIGHT-(yy-Ly);
	PL[1][0]=xx+Lx;PL[1][1]=CRT_HIGHT-(yy-Ly);
	PL[2][0]=xx+Lx+Lz;PL[2][1]=CRT_HIGHT-(yy-Ly-Lz);
	PL[3][0]=xx+Lz;PL[3][1]=CRT_HIGHT-(yy-Ly-Lz);
	GWpolygon(PL,4,0);
	PL[0][0]=xx+Lx;PL[0][1]=CRT_HIGHT-yy;
	PL[1][0]=xx+Lx;PL[1][1]=CRT_HIGHT-(yy-Ly);
	PL[2][0]=xx+Lx+Lz;PL[2][1]=CRT_HIGHT-(yy-Ly-Lz);
	PL[3][0]=xx+Lx+Lz;PL[3][1]=CRT_HIGHT-(yy-Lz);
	GWpolygon(PL,4,0);
	PL[0][0]=xx;PL[0][1]=CRT_HIGHT-yy;
	PL[1][0]=xx+Lx;PL[1][1]=CRT_HIGHT-yy;
	PL[2][0]=xx+Lx;PL[2][1]=CRT_HIGHT-(yy-Ly);
	PL[3][0]=xx;PL[3][1]=CRT_HIGHT-(yy-Ly);
	GWpolygon(PL,4,0);
	return;
}
G_H_line(xs,ys,color,length)
{
	if( ys<0 || ys>=CRT_HIGHT || xs>=CRT_WIDE ) return;
	if(xs<0) {
		xs=0;length+=xs;
		if( length<=0 ) return;
	}
	if( xs+length > CRT_WIDE ) length=CRT_WIDE-xs;
	G_horizontal_line(xs,ys,color,length);
}
G_H_H_line(xs,ys,color,length)
{
int	i,j;

	if( Screen[xs]>=ys ){
		if( Screen[xs+length-1] >= ys ) G_H_line(xs,ys,color,length);
			else {
				for(i=xs+length-2;i>=xs;i--){
					length--;
					if( Screen[i] >= ys ) break;
				}
				G_H_line(xs,ys,color,length);
			}
	} else {
		if( Screen[xs+length-1] >= ys ){
			j=xs+length-1;
			for(i=xs+1;i<j;i++){
				length--;
				if( Screen[i] >= ys ) break;
			}
			xs=i;
			G_H_line(xs,ys,color,length);
		}
	}
}
G_V_line(xs,ys,color,length)
{
	if( xs<0 || xs>=CRT_WIDE || ys>=CRT_HIGHT ) return;
	if(ys<0) {
		ys=0;length+=ys;
		if(length<=0) return;
	}
	if( ys+length > CRT_HIGHT ) length=CRT_HIGHT-ys;
	G_vertical_line(xs,ys,color,length);
}
G_H_V_line(xs,ys,color,length)
{
int	i;

	if( (i=Screen[xs-1]) >= ys ){
		if( i < ys+length ) length=i-ys+1;
		G_V_line(xs,ys,color,length);
	}
}
G_H_PSET(xs,ys,color)
{
	if( Screen[xs] >= ys ) G_PSET(xs,ys,color);
}
static	int		crt_shift;
static	double	norm;
graph_init()
{
	int i;
G_WINDOW(CRT_WIDE,CRT_HIGHT);
	G_INIT();
//getch();
//GWmode(1,1);
	GWarrange(6);
//	GWsize(3,&O_IW,&O_IH);
        G_PALETTE(3,0,0,0);
        G_PALETTE(1,0,255,0);
        G_PALETTE(2,255,0,0);
        G_PALETTE(0,255,255,255);
        G_PALETTE(4,0,0,255);
        G_PALETTE(5,255,255,0);
        G_PALETTE(6,0,255,255);
        G_PALETTE(7,255,0,255);

        G_PALETTE(8,200,200,200);
        G_PALETTE(9,0,127,0);
        G_PALETTE(10,127,0,0);
        G_PALETTE(11,100,100,100);
        G_PALETTE(12,0,0,127);
        G_PALETTE(13,127,127,0);
        G_PALETTE(15,127,0,127);
        G_PALETTE(14,0,127,127);
	GWmode(1,0);
	G_CLS();
	PSET_BI(PSET_SIZE);
	norm=(double)(CRT_WIDE-crt_shift)/sizex/sizex;
	crt_shift=PSET_SIZE*sizex*4;
	G_CLS();
}
#endif