随机行走程序

INTEGERM(1:10000),NUMBER1(0:364),NUMBER2REALX,YISEED=RTC()DOJ=1,10000NUMBER1=0X=RAN(ISEED)NUMBER1(0)=INT(365*X+1)JJJ=1DOI=1,364Y=RAN(ISEED)NUMBER2=INT(365*Y+1)ETR=COUNT(NUMBER1.EQ.NUMBER2)IF(ETR==1)THENEXITELSEJJJ=JJJ+1M(J)=JJJNUMBER1(I)=NUMBER2ENDIFENDDOENDDODOI=1,10000IF(M(I).LE.23)SUM=SUM+1ENDDOPRINT*,SUM/10000ENDMonteCarloSimulationofOneDimensionalDiffusionINTEGERX,XX(1:1000,1:1000)REALXXM(1:1000)!X:INSTANTANEOUSPOSITIONOFATOM!XX(J,I)：X*X,J:第几天实验,I:第几步跳跃!XXM(I):THEMEANOFXXWRITE(*,*)实验天数JMAX，实验次数IMAXREAD(*,*)JMAX,IMAXISEED=RTC()DOJ=1,JMAX!第几天实验X=0!!!DOI=1,IMAX!第几步跳跃RN=RAN(ISEED)IF(RN0.5)THENX=X+1ELSEX=X-1ENDIFXX(J,I)=X*XENDDOENDDOOPEN(1,FILE=C:\DIF1.DAT)DOI=1,IMAXXXM=0.0XXM(I)=1.0*SUM(XX(1:JMAX,I))/JMAX!!WRITE(1,*)I,XXM(I)ENDDOCLOSE(1)END!MonteCarloSimulationofTwoDimensionalDiffusionINTEGERX,Y,XY(1:1000,1:1000)REALXYM(1:1000)!X:INSTANTANEOUSPOSITIONOFATOM!XY(J,I)：X*Y,J:第几天实验,I:第几步跳跃!XYM(I):THEMEANOFXYWRITE(*,*)实验天数JMAX，实验次数IMAXREAD(*,*)JMAX,IMAXISEED=RTC()DOJ=1,JMAX!第几天实验X=0!!!Y=0!!!DOI=1,IMAX!第几步跳跃RN=RAN(ISEED)IF(RN.LT.0.25)THENx=xy=y-1ENDIFIF(RN.LT.0.5.AND.RN.GE.0.25)THENx=xy=y+1ENDIFIF(RN.LT.0.75.AND.RN.GE.0.5)THENx=x-1y=yENDIFIF(RN.GE.0.75)THENx=x+1y=yENDIFXY(J,I)=X*X+Y*YENDDOENDDOOPEN(1,FILE=C:\DIF2.DAT)DOI=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX!!WRITE(1,*)I,XYM(I)ENDDOCLOSE(1)END!MonteCarloSimulationofOneDimensionalDiffusionINTEGERX,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RNREALXYM(1:1000)!X:INSTANTANEOUSPOSITIONOFATOM!XY(J,I)：X*Y,J:第几天实验,I:第几步跳跃!XYM(I):THEMEANOFXYWRITE(*,*)实验天数JMAX，实验次数IMAXREAD(*,*)JMAX,IMAXXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)ISEED=RTC()DOJ=1,JMAX!第几天实验X=0!!!Y=0!!!DOI=1,IMAX!第几步跳跃RN=4*RAN(ISEED)+1X=X+YN(RN)Y=Y+YN(RN)XY(J,I)=X*X+Y*YENDDOENDDOOPEN(1,FILE=C:\DIF2.DAT)DOI=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX!!WRITE(1,*)I,XYM(I)ENDDOCLOSE(1)END做三维空间随机行走？？留作业!MonteCarloSimulationofOneDimensionalDiffusionINTEGERX,XY(1:1000,1:1000),y,XN(1:6),YN(1:6),ZN(1:6),RNREALXYM(1:1000)!X:INSTANTANEOUSPOSITIONOFATOM!XY(J,I)：X*Y,J:第几天实验,I:第几步跳跃!XYM(I):THEMEANOFXYWRITE(*,*)实验天数JMAX，实验次数IMAXREAD(*,*)JMAX,IMAXXN=(/0,0,-1,1,0,0/)YN=(/-1,1,0,0,0,0/)ZN=(/0,0,0,0,1,-1/)ISEED=RTC()DOJ=1,JMAX!第几天实验X=0!!!Y=0!!!Z=0DOI=1,IMAX!第几步跳跃RN=6*RAN(ISEED)+1X=X+XN(RN)Y=Y+YN(RN)Z=Z+ZN(RN)XY(J,I)=X*X+Y*Y+Z*ZENDDOENDDOOPEN(1,FILE=C:\DIF2.DAT)DOI=1,IMAXXYM=0.0XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX!!WRITE(1,*)I,XYM(I)ENDDOCLOSE(1)END1、随机行走程序模拟的意义2、演示二维随机行走一次随机抽样的结果3、如果随机行走是有约束的，计算方法应该如何设计4、讲解三维随机行走通过该程序了解fortran语言如何画图（通过像素画图）USEMSFLIBINTEGERXR，YR!在的区域中画一个圆PARAMETERXR=400,YR=400INTEGERR,S(1:XR,1:YR)X0=XR/2!圆心位置X0，YOY0=YR/2R=MIN(X0-10,Y0-10)!圆半径S=0!像素的初始状态（颜色）DOI=1,XRDOJ=1,YRIF((I-X0)**2+(J-Y0)**2=R**2)S(I,J)=10IER=SETCOLOR(S(I,J))IER=SETPIXEL(I,J)ENDDOENDDOEND!画一个圆（1、如何选出晶界区域；2、进一步加深对画图的理解）USEMSFLIBINTEGERXR，YR!在的区域中画一个圆PARAMETERXR=400,YR=400INTEGERR,S(0:XR+1,0:YR+1),XN(1:4),YN(1:4),SNSXN=(/0,0,-1,1/)YN=(/-1,1,0,0/)X0=XR/2!圆心位置X0，Y0Y0=YR/2R=MIN(X0-10,Y0-10)!圆半径S=0!像素的初始状态（颜色）DOI=1,XRDOJ=1,YRIF((I-X0)**2+(J-Y0)**2=R**2)S(I,J)=10IER=SETCOLOR(S(I,J))IER=SETPIXEL(I,J)ENDDOENDDODOI=1,XR!画晶界DOJ=1,YRNDS=0DOK=1,4IF(S(I,J).NE.S(I+XN(K),J+YN(K)))NDS=NDS+1ENDDOIF(NDS0)THENIER=SETCOLOR(9)ELSEIER=SETCOLOR(8)ENDIFIER=SETPIXEL(I,J)ENDDOENDDOEND如何画有一定宽度的晶界？！MC模拟一个晶粒的缩小（1、如何定义基体和晶粒；（2、如何寻找边界；（3、如何计算能量（构造或描述问题的概率过程）步骤：1、在基体上画一个原始晶粒，并赋予基体和原始晶粒不同的状态（取向号或体积能）2、寻找晶界3、MC的一个时间步（晶粒长大过程中一种随机性）4、计算晶界上网格点的相互作用，每个格点的相互作用导致晶粒长大①随机选取一个初始格点；②若此点属于晶界，那么可以随机转变为它邻居的状态，若不是晶界，则跳出循环，不发生转变；③计算转变前后的能量变化，⊿E=⊿Ev+⊿Es+⊿Eq(自由能=体积能+表面能+能量起伏)④若⊿E小于0，则新晶相被接受，网格状态发生转变。USEMSFLIBPARAMETERIR=400,JR=400INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IYWRITE(*,*)PLEASEINPUTTHETIMESTEPREAD(*,*)TMAXISEED=RTC()!定义圆心和半径IRC=IR/2JRC=JR/2R=MIN(IRC,JRC)-10!定义基体和圆晶粒分别为状态1、状态2IS=1DOI=1,IRDOJ=1,JRDISTANCE=SQRT(1.0*(I-IRC)**2+1.0*(J-JRC)**2)IF(DISTANCE.LT.R)IS(I,J)=2ISE=SETCOLOR(IS(I,J))ISE=SETPIXEL(I,J)ENDDOENDDOOPEN(1,FILE=E:\LUKE.DAT)!寻找晶粒边界，计算能量，改变状态。DOT=1,TMAXDOX=1,IRDOY=1,JRIX=IR*RAN(ISEED)+1JY=JR*RAN(ISEED)+1ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1),IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)E0=COUNT(ISN.NE.IS(IX,JY))!如果不是圆晶粒边界，则跳出，重新循环IF(E0.EQ.0)CYCLE!随机寻找一个相邻点NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态E=COUNT(ISN.NE.NSTATE)RD=RAN(ISEED)DE=E-E0+NSTATE-IS(IX,JY)+2.5*RD-1.25IF(DE.LT.0.0)IS(IX,JY)=NSTATEISRE=SETCOLOR(IS(IX,JY))ISRE=SETPIXEL(IX,JY)ENDDOENDDOWRITE(1,*)T,COUNT(IS.EQ.2)ENDDOCLOSE(1)END（1、演示体积能互换的转变情况；2、演示IX=IR*RAN(ISEED)+1JY=JR*RAN(ISEED)+1的替换；3、能量关系变化；）作业3：基体上单晶粒形核长大（1、选取初始形核点（2、定义体积能的变化）USEMSFLIBPARAMETERIR=400,JR=400INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IYWRITE(*,*)PLEASEINPUTTHETIMESTEPREAD(*,*)TMAXISEED=RTC()!定义圆心和半径IRC=IR/2JRC=IR/2!定义基体和圆晶粒分别为状态10、状态2IS=10IS(IRC,JRC)=2OPEN(1,FILE=E:\LUKE.DAT)!寻找晶粒边界，计算能量，改变状态。DOT=1,TMAXDOX=1,IRDOY=1,JRIX=IR*RAN(ISEED)+1JY=JR*RAN(ISEED)+1ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)!,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)E0=COUNT(ISN.NE.IS(IX,JY))!如果不是圆晶粒边界，则跳出，重新循环IF(E0.EQ.0)CYCLE!随机寻找一个相邻点NR=8*RAN(ISEED)+1NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态E=COUN