fortran程序汇总.doc

计算圆周率REAL R,R1,R2,PIISEED=RTC()N0=0N=300000DO I=1,NR1=RAN(ISEED)R2=RAN(ISEED)R=SQRT(R1*R1+R2*R2)IF(R1.0)N0=N0+1END DOPI=4.0*N0/NWRITE(*,*)PIEND一）蒙特卡洛计算生日问题假设有N个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直觉相反，只需23个人便有大于50的几率人群中至少有2个人生日相同。INTEGER M(1:10000), NUMBER1(0:364), NUMBER2REAL X,YISEED=RTC()DO J=1, 10000NUMBER1=0X=RAN(ISEED)NUMBER1(0)=INT(365*X+1)JJJ=1 DO I=1,364 Y=RAN(ISEED) NUMBER2=INT(365*Y+1) ETR=COUNT(NUMBER1.EQ.NUMBER2) IF (ETR= =1) THEN EXIT ELSE JJJ=JJJ+1 M(J)=JJJ NUMBER1(I)=NUMBER2 END IF END DOEND DO DO I=1,10000 IF(M(I).LE.23) SUM=SUM+1END DO PRINT *,SUM/10000 END 二）MONTE CARLO SIMULATION OF ONE DIMENSIONAL DIFFUSION 蒙特卡罗计算一维扩散问题INTEGER X,XX(1:1000,1:1000)REAL XXM(1:1000)!X:INSTANTANEOUS POSITION OF ATOM!XX(J,I)：X*X ,J:第几天实验,I:第几步跳跃!XXM(I): THE MEAN OF XX WRITE(*,*) 实验天数JMAX，实验次数 IMAXREAD(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第几天实验X=0 !DO I=1,IMAX !第几步跳跃RN=RAN(ISEED)IF(RN0.5)THEN X=X+1ELSE X=X-1 END IFXX(J,I)=X*XEND DOEND DOOPEN(1,FILE=C:DIF1.DAT)DO I=1,IMAXXXM=0.0XXM(I)=1.0*SUM(XX(1:JMAX,I)/JMAX !WRITE(1,*) I, XXM(I)END DOCLOSE(1)END三维的！三）通过该程序了解FORTRAN语言如何画图（通过像素画图）USE MSFLIB INTEGER XR，YR !在的区域中画一个圆PARAMETER XR=400,YR=400 INTEGER R, S(1:XR,1:YR) X0=XR/2 ! 圆心位置X0，YOY0=YR/2R=MIN(X0-10,Y0-10) !圆半径S=0 !像素的初始状态（颜色）DO I=1,XRDO J=1,YRIF(I-X0)*2+(J-Y0)*2=R*2)S(I,J)=10IER=SETCOLOR(S(I,J) IER=SETPIXEL(I,J)END DOEND DOEND四）画一个圆（1、如何选出晶界区域；2、进一步加深对画图的理解）USE MSFLIB INTEGER XR，YR !在的区域中画一个圆PARAMETER XR=400,YR=400 INTEGER R, S(0:XR+1,0:YR+1), XN(1:4), YN(1:4), SNS XN=(/0,0,-1,1/)YN=(/-1,1,0,0/)X0=XR/2 ! 圆心位置X0，Y0Y0=YR/2R=MIN(X0-10,Y0-10) !圆半径S=0 !像素的初始状态（颜色）DO I=1,XRDO J=1,YRIF(I-X0)*2+(J-Y0)*20)THENIER=SETCOLOR(9)ELSEIER=SETCOLOR(8)END IF IER=SETPIXEL(I,J)END DOEND DOEND五）MC模拟一个晶粒的缩小USE MSFLIBPARAMETER IR=400,JR=400 INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IYWRITE(*,*)PLEASE INPUT THE TIME STEP READ(*,*)TMAXISEED=RTC()!定义圆心和半径IRC=IR/2JRC=IR/2R=MIN(IRC,JRC)-10!定义基体和圆晶粒分别为状态1、状态2IS=1 DO I=1,IR DO J=1,JR DISTANCE=SQRT(1.0*(I-IRC)*2+1.0*(J-JRC)*2) IF(DISTANCE.LT.R)IS(I,J)=2 ISE=SETCOLOR(IS(I,J) ISE=SETPIXEL(I,J) END DO END DOOPEN(1,FILE=E:LUKE.DAT)!寻找晶粒边界，计算能量，改变状态。 DO T=1,TMAX DO X=1,IR DO Y=1,JR IX=IR*RAN(ISEED)+1 JY=JR*RAN(ISEED)+1 ISN=(/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)+1 NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态 E=COUNT(ISN.NE.NSTATE) RD=RAN(ISEED) DE=E-E0+NSTATE-IS(IX,JY)+2.5*RD-1.25 IF(DE.LT.0.0)IS(IX,JY)=NSTATE ISRE=SETCOLOR(IS(IX,JY) ISRE=SETPIXEL(IX,JY) ENDDO ENDDO WRITE(1,*)T,COUNT(IS.EQ.2)ENDDO CLOSE(1)END 六）蒙特卡罗模拟基体上单晶粒形核长大 USE MSFLIBPARAMETER IR=400,JR=400 INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IYWRITE(*,*)PLEASE INPUT THE TIME STEP READ(*,*)TMAXISEED=RTC()!定义圆心和半径IRC=IR/2JRC=IR/2 !定义基体和圆晶粒分别为状态10、状态2IS=10IS(IRC,JRC)=2OPEN(1,FILE=E:LUKE.DAT)!寻找晶粒边界，计算能量，改变状态。 DO T=1,TMAX DO X=1,IR DO Y=1,JR IX=IR*RAN(ISEED)+1 JY=JR*RAN(ISEED)+1 ISN=(/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)+1 NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态 E=COUNT(ISN.NE.NSTATE) RD=RAN(ISEED) DE=E-E0+NSTATE-IS(IX,JY)+2.5*RD-1.25 IF(DE.LT.0.0)IS(IX,JY)=NSTATE ISRE=SETCOLOR(IS(IX,JY) ISRE=SETPIXEL(IX,JY) ENDDO ENDDO WRITE(1,*)T,SQRT(1.0*COUNT(IS.EQ.2)ENDDOCLOSE(1)END 七）蒙特卡洛模拟基体上多晶粒形核长大， USE MSFLIB！定义常数名PARAMETER IR=400,JR=400,NMAX=100！定义变量：体积能变量（一维数组）；基体各像素点变量（二维数组），邻居的取向号（一维数组）等； INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,IYINTEGER IGV(0:NMAX)！读取晶粒生长步长；WRITE(*,*)PLEASE INPUT THE TIME STEP READ(*,*)TMAXISEED=RTC()!定义基体体积能为10,所有晶粒体积能为1IGV=1IGV(0)=10!定义晶粒长大位置（100个形核点初始形核位置），并附已不同的取向号DO I=1, NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1 IS(IX0,JY0)=IEND DO!边界条件 IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)OPEN(1,FILE=E:LUKE.DAT)!寻找晶粒边界。 DO T=1,TMAX DO X=1,IR DO Y=1,JR IX=IR*RAN(ISEED)+1 JY=JR*RAN(ISEED)+1 ISN=(/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)+1 NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态 RD=RAN(ISEED) E=COUNT(ISN.NE.NSTATE) IG=IGV(NSTATE)-IGV(IS(IX,JY) DE=E-E0+IG+2.5*RD-1.25 IF(DE.LT.0.0)IS(IX,JY)=NSTATE ISRE=SETCOLOR(MOD(IS(IX,JY),15) ISRE=SETPIXEL(IX,JY) ENDDO ENDDO WRITE(1,*)T,1.0*COUNT(IS.NE.0)ENDDOCLOSE(1)END 八）元胞自动机模拟生命游戏程序（生命永不停止）USE MSFLIBPARAMETER IR=400,JR=400,NMAX=40000 !NMAX-随机产生的生命种子INTEGER IS(0:1001,0:1001),IS1(0:1001,0:1001),ISN(1:8), TMAX, NUM ！IS-基体的二维数组PRINT*,PLEASE INPUT LOOP(TMAX)READ*, TMAXISEED=RTC()IS=15 ！“死”的状态，基体为白色 ！赋予生命的种子，“活”的状态1DO I=1, NMAXIX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1 IS(IX0,JY0)=1END DO!EXECUTE THE RULEDO T=1,TMAXIS1=IS！边界条件 IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)！搜索生命存在的位置DO IX=1,IRDO JY=1,JR！判断邻居状态 ISN=(/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)/)NUM=COUNT(ISN.EQ.1)！赋予生存的条件IF(IS(IX,JY)=15.AND.NUM=3).OR.(IS(IX,JY)=1.AND.(NUM=3.OR.NUM=2) THEN IS1(IX,JY)=1 ELSE IS1(IX,JY)=15 END IF ！画图 ISRE=SETCOLOR(IS1(IX,JY) ISRE=SETPIXEL(IX,JY)END DO END DOIS=IS1END DOEND九）元胞自动机模拟单晶长大 USE MSFLIBPARAMETER IR=400,JR=400 INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY! 根据过去状态IS，定义现在状态为IS1；为了突出边界，特别定义ISN1 INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8) WRITE(*,*)PLEASE INPUT THE TIME STEP READ(*,*)TMAXISEED=RTC() IRC=IR/2 !=IR*ran(iseed)+1JRC=JR/2!定义基体体积能为10,晶粒体积能为1IS=8IS(IRC,JRC)=1! 在循环前定义现在状态数组IS1的初始值 IS1=IS OPEN(1,FILE=E:LUKE.DAT) DO T=1,TMAX! 每次循环前重新定义过去状态数组IS IS=IS1!边界条件 IS(0,0:JR+1)=IS(IR,0:JR+1)IS(IR+1,0:JR+1)=IS(1,0:JR+1)IS(0:IR+1,0)=IS(0:IR+1,JR)IS(0:IR+1,JR+1)=IS(0:IR+1,1)! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变 DO IX=1,IR DO JY=1,JR ISN=(/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)+1 NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态 E=COUNT(ISN.NE.NSTATE) RD=RAN(ISEED) IG=NSTATE-IS(IX,JY) DE=E-E0+IG+2.5*RD-1.25! 用现在状态数组IS1记录边界状态的改变 IF(DE.LT.0.0)IS1(IX,JY)=NSTATE ENDDO END DO! 每循环20次在显示屏幕上刷新状态（颜色） DO IX=1,IR DO JY=1,JR! IF(MOD(T,20).EQ.0)THEN ISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1) & ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/) ISRE=SETCOLOR(IS(IX,JY)! 如果是边界，定义特别颜色15,加以区分 IF(COUNT(ISN1.NE.IS1(IX,JY).GT.0) ISRE=SETCOLOR(15) ISRE=SETPIXEL(IX,JY) ! END IF ENDDO ENDDO! 记录循环次数和对应的晶粒面积 WRITE(1,*)T, SQRT(1.0*COUNT(IS.EQ.1)ENDDOCLOSE(1)END 十）元胞自动机模拟多晶长大1.介绍蒙特卡罗和元胞自动机的区别。n 元胞自动机特点：时间、空间、状态都离散的动力学模型，n 利用确定性或概论性的转换规则来描述在离散空间和时间上复杂系统演化 USE MSFLIB PARAMETER IR=400,JR=400,NMAX=100 INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY INTEGER IGV(0:NMAX)! 根据过去状态IS，定义现在状态为IS1；为了突出边界，特别定义ISN1 INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8) WRITE(*,*)PLEASE INPUT THE TIME STEP READ(*,*)TMAXISEED=RTC() !定义基体体积能为10,晶粒体积能为1IGV=1IGV(0)=10!定义晶粒长大位置，并附能量，附带生长取向（研究形核数目与晶粒尺寸的关系）DO I=1,NMAX IX0=IR*RAN(ISEED)+1JY0=JR*RAN(ISEED)+1IF(IS(IX0,JY0).NE.0)CYCLE IS(IX0,JY0)=IEND DO! 在循环前定义现在状态数组IS1的初始值 IS1=IS! OPEN(1,FILE=E:LUKE.DAT) DO T=1,TMAX! 每次循环前重新定义过去状态数组IS IS=IS1!边界条件 IS(0,0:JR+1)=IS(IR,0:JR+1)IS(IR+1,0:JR+1)=IS(1,0:JR+1)IS(0:IR+1,0)=IS(0:IR+1,JR)IS(0:IR+1,JR+1)=IS(0:IR+1,1)!寻找晶粒边界，计算能量，改变状态。 ! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变 DO IX=1,IR DO JY=1,JR! ISN=(/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)+1 NSTATE=ISN(NR)!判断与相邻点的能量差，并决定是否改变状态 E=COUNT(ISN.NE.NSTATE) RD=RAN(ISEED) IG=IGV(NSTATE)-IGV(IS(IX,JY) DE=E-E0+IG+2.5*RD-1.25! 用现在状态数组IS1记录边界状态的改变 IF(DE.LT.0.0)IS1(IX,JY)=NSTATE! ISRE=SETCOLOR(MOD(IS1(IX,JY),15)! ISRE=SETPIXEL(IX,JY) ENDDO ENDDO! 每循环20次在显示屏幕上刷新状态（颜色） DO IX=1,IR DO JY=1,JR IF(MOD(T,20).EQ.0)THEN ISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1) & ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/) ISRE=SETCOLOR(MOD(IS1(IX,JY),15)! 如果是边界，定义特别颜色15,加以区分 IF(COUNT(ISN1.NE.IS1(IX,JY).GT.0) ISRE=SETCOLOR(15) ISRE=SETPIXEL(IX,JY) END IF ENDDO ENDDO! 记录循环次数和对应的晶粒面积 WRITE(1,*)T, SQRT(1.0*COUNT(IS1.GT.0)ENDDOCLOSE(1)END 十一）相场方法模拟相变 BY USING BOUNDARY TRACKING METHORD,2008,11.23! GRAIN GROWTH VELOCITY IS 1.0*(IT-NT(IST) OR A*(IT-NT(IST)*B!WHICH SHOULD BE NOT GREETER THAN 1 CELL EACH STEP. USE MSFLIBPARAMETER IMAX=800,JMAX=800,NMAX=50000,NSTEP=10 ! 最多晶核数,每步形成晶核数INTEGER IS(0:IMAX+1,0:JMAX+1),IST(0:IMAX+1,0:JMAX+1)INTEGER IN(1:8),JN(1:8)INTEGER NI(1:NMAX),NJ(1:NMAX),NT(1:NMAX) !核的I,J坐标,形成时间INTEGER NNS,DCN,DCNM !运行时间部，最近晶核状态,访问CELL 与晶核的距离平方，最短距离平方 ISEED=RTC()WRITE(*,*)PLEASE INPUT THE TIME STEPREAD(*,*)TMAXIN=(/-1,-1,-1,0,0,1,1,1/)JN=(/-1,0,1,-1,1,-1,0,1/)OPEN(5,FILE=C:IRX.DAT)DO IT=1,TMAXIS=IST!边界条件IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX+1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX+1,0)=IS(0:IMAX+1,JMAX)IS(0:IMAX+1,JMAX+1)=IS(0:IMAX+1,1)!-形核DO I=1,NSTEPINC=RAN(ISEED)*IMAX+1JNC=RAN(ISEED)*JMAX+1IF(IST(INC,JNC).NE.0)CYCLENN=NN+1NI(NN)=INCNJ(NN)=JNCNT(NN)=IT IST(INC,JNC)=NN! END DOC -DO I=1,IMAXDO J=1,JMAXIF(IS(I,J).GT.0)CYCLEDCNM=2*IMAX*JMAXNP=0DO K=1,8II=I+IN(K) JJ=J+JN(K) IF(IS(II,JJ).EQ.0)CYCLE !效率提高 IF(NT(IS(II,JJ).EQ.IT)CYCLE ! 当前深刻形成的核,当前时刻不长大NP=NP+1 ! 只有IT 时刻以前形成的核长大 IDN=I-NI(IS(II,JJ) JDN=J-NJ(IS(II,JJ)IDNABS=ABS(IDN)JDNABS=ABS(JDN)DCN=IDN*IDN+JDN*JDN IF(DCN.LE.DCNM)THEN DCNM=DCN !访问CELL与最近晶核的距离平方NNS=IS(II,JJ) ! END IFEND DOIF(NP.EQ.0)CYCLEGR=1.0*(IT-NT(NNS) !晶粒半径DIS=GR-SQRT(1.0*DCNM)