当前位置：文档之家› 几个fortran程序

几个fortran程序

计算圆周率REAL R,R1,R2,PIISEED=RTC()N0=0N=300000DO I=1,NR1=RAN(ISEED)R2=RAN(ISEED)R=SQRT(R1*R1+R2*R2)IF(R<1.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=1DO I=1,364Y=RAN(ISEED)NUMBER2=INT(365*Y+1)ETR=COUNT(NUMBER1.EQ.NUMBER2)IF (ETR= =1) THENEXITELSEJJJ=JJJ+1M(J)=JJJNUMBER1(I)=NUMBER2END IFEND DOEND DODO I=1,10000IF(M(I).LE.23) SUM=SUM+1END DOPRINT *,SUM/10000END二）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 XXWRITE(*,*) "实验天数JMAX，实验次数IMAX"READ(*,*) JMAX,IMAXISEED=RTC()DO J=1,JMAX !第几天实验X=0 !!!DO I=1,IMAX !第几步跳跃RN=RAN(ISEED)IF(RN<0.5)THENX=X+1ELSEX=X-1END 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 MSFLIBINTEGER XR，YR !在的区域中画一个圆PARAMETER XR=400,YR=400INTEGER 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 MSFLIBINTEGER XR，YR !在的区域中画一个圆PARAMETER XR=400,YR=400INTEGER R, 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 !像素的初始状态（颜色）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 DODO I=1,XR !画晶界DO J=1,YRNDS=0DO K=1,4IF(S(I,J).NE.S(I+XN(K),J+YN(K)))NDS=NDS+1END DOIF(NDS>0)THENIER=SETCOLOR(9)ELSEIER=SETCOLOR(8)END IFIER=SETPIXEL(I,J)END DOEND DOEND五）MC模拟一个晶粒的缩小USE MSFLIBPARAMETER IR=400,JR=400INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IY WRITE(*,*)"PLEASE INPUT THE TIME STEP "READ(*,*)TMAXISEED=RTC()! 定义圆心和半径IRC=IR/2JRC=IR/2R=MIN(IRC,JRC)-10! 定义基体和圆晶粒分别为状态1、状态2IS=1DO I=1,IRDO J=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)END DOEND DOOPEN(1,FILE="E:\LUKE.DAT")! 寻找晶粒边界，计算能量，改变状态。

DO T=1,TMAXDO X=1,IRDO Y=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,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,IY INTEGER 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)+1IS(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,TMAXDO X=1,IRDO Y=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)! 判断与相邻点的能量差，并决定是否改变状态RD=RAN(ISEED)E=COUNT(ISN.NE.NSTATE)IG=IGV(NSTATE)-IGV(IS(IX,JY))DE=E-E0+IG+2.5*RD-1.25IF(DE.LT.0.0)IS(IX,JY)=NSTATEISRE=SETCOLOR(MOD(IS(IX,JY),15))ISRE=SETPIXEL(IX,JY)ENDDOENDDOWRITE(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)+1IS(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))) THENIS1(IX,JY)=1ELSEIS1(IX,JY)=15END IF！画图ISRE=SETCOLOR(IS1(IX,JY))ISRE=SETPIXEL(IX,JY)END DOEND DOIS=IS1END DO九）元胞自动机模拟单晶长大USE MSFLIBPARAMETER IR=400,JR=400INTEGER 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=ISOPEN(1,FILE="E:\LUKE.DAT")DO T=1,TMAX!! 每次循环前重新定义过去状态数组ISIS=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,IRDO JY=1,JRISN=(/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)IG=NSTATE-IS(IX,JY)DE=E-E0+IG+2.5*RD-1.25!! 用现在状态数组IS1记录边界状态的改变IF(DE.LT.0.0)IS1(IX,JY)=NSTATEENDDOEND DO!! 每循环20次在显示屏幕上刷新状态（颜色）DO IX=1,IRDO JY=1,JR! IF(MOD(T,20).EQ.0)THENISN1=(/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 IFENDDOENDDO!!! 记录循环次数和对应的晶粒面积WRITE(1,*)T, SQRT(1.0*COUNT(IS.EQ.1))ENDDOCLOSE(1)END十）元胞自动机模拟多晶长大1.介绍蒙特卡罗和元胞自动机的区别。

e商务文档

几个fortran程序

相关文档推荐：