fortran程序汇总(1)

计算圆周率  REALR,R1,R2,PIISEED=RTC()N0=0N=300000DOI=1,NR1=RAN(ISEED)R2=RAN(ISEED)R=SQRT(R1*R1R2*R2)IF(R<1.0)N0=N01  ENDDOPI=4.0*N0/NWRITE(*,*)PIEND一）蒙特卡洛计算生日问题假设有N个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直觉相反，只需23个人便有大于50％的几率人群中至少有2个人生日相同。INTEGER M(1:10000),NUMBER1(0:364),NUMBER2REALX,YISEED=RTC()DOJ=1,10000NUMBER1=0X=RAN(ISEED)NUMBER1(0)=INT(365*X1)JJJ=1DOI=1,364Y=RAN(ISEED)NUMBER2=INT(365*Y1)ETR=COUNT(NUMBER1.EQ.NUMBER2)IF(ETR==1)THENEXITELSEJJJ=JJJ1M(J)=JJJNUMBER1(I)=NUMBER2ENDIFENDDO  ENDDODOI=1,10000IF(M(I).LE.23)SUM=SUM1ENDDOPRINT*,SUM/10000END    二）MONTECARLOSIMULATIONOFONEDIMENSIONALDIFFUSION蒙特卡罗计算一维扩散问题INTEGERX,XX(1:1000,1:1000)REALXXM(1:1000)!  X:INSTANTANEOUSPOSITIONOFATOM!  XX(J,I)：X*X,J:第几天实验,I:第几步跳跃!  XXM(I):THEMEANOFXXWRITE(*,*)"实验天数JMAX，实验次数IMAX"READ(*,*)JMAX,IMAXISEED=RTC()DO J=1,JMAX!第几天实验X=0!!!DO I=1,IMAX!第几步跳跃RN=RAN(ISEED)IF(RN<0.5)THENX=X1ELSEX=X-1  ENDIF      XX(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三维的！三）通过该程序了解FORTRAN语言如何画图（通过像素画图）USEMSFLIB  INTEGER XR，YR !在的区域中画一个圆PARAMETERXR=400,YR=400INTEGER R,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、进一步加深对画图的理解）USEMSFLIBINTEGER XR，YR !在的区域中画一个圆PARAMETERXR=400,YR=400INTEGERR,S(0:XR1,0:YR1),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(IXN(K),JYN(K)))NDS=NDS1ENDDO    IF(NDS>0)THENIER=SETCOLOR(9)ELSEIER=SETCOLOR(8)ENDIFIER=SETPIXEL(I,J)ENDDOENDDOEND五）MC模拟一个晶粒的缩小USEMSFLIBPARAMETERIR=400,JR=400INTEGER IS(0:IR1,0:JR1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IYWRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()!  定义圆心和半径IRC=IR/2JRC=IR/2R=MIN(IRC,JRC)-10  !  定义基体和圆晶粒分别为状态1、状态2IS=1DOI=1,IRDOJ=1,JRDISTANCE=SQRT(1.0*(I-IRC)**21.0*(J-JRC)**2)IF(DISTANCE.LT.R)IS(I,J)=2ISE=SETCOLOR(IS(I,J))ISE=SETPIXEL(I,J)ENDDO     ENDDOOPEN(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,JY1),IS(IX,JY-1),IS(IX,JY1),IS(IX1,JY-1),IS(IX1,JY),IS(IX1,JY1)/)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-E0NSTATE-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   六）蒙特卡罗模拟基体上单晶粒形核长大USEMSFLIBPARAMETERIR=400,JR=400INTEGER IS(0:IR1,0:JR1),TMAX,ISN(1:8),NSTATE,T,NR,IX,IYWRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)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,JY1),IS(IX,JY-1),IS(IX,JY1),IS(IX1,JY-1),IS(IX1,JY),IS(IX1,JY1)/)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-E0NSTATE-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  七）蒙特卡洛模拟基体上多晶粒形核长大，USEMSFLIB！定义常数名PARAMETERIR=400,JR=400,NMAX=100！定义变量：体积能变量（一维数组）；基体各像素点变量（二维数组），邻居的取向号（一维数组）等；INTEGER IS(0:IR1,0:JR1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,IYINTEGERIGV(0:NMAX)！读取晶粒生长步长；WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()!  定义基体体积能为10,所有晶粒体积能为1IGV=1IGV(0)=10!  定义晶粒长大位置（100个形核点初始形核位置），并附已不同的取向号DOI=1,NMAXIX0=IR*RAN(ISEED)1JY0=JR*RAN(ISEED)1  IS(IX0,JY0)=IENDDO  !  边界条件IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX1,0)=IS(0:IMAX1,JMAX)IS(0:IMAX1,JMAX1)=IS(0:IMAX1,1)OPEN(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,JY1),IS(IX,JY-1)&      ,IS(IX,JY1),IS(IX1,JY-1),IS(IX1,JY),IS(IX1,JY1)/)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-E0IG2.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   八）元胞自动机模拟生命游戏程序（生命永不停止）USEMSFLIBPARAMETERIR=400,JR=400,NMAX=40000!NMAX-随机产生的生命种子INTEGERIS(0:1001,0:1001),IS1(0:1001,0:1001),ISN(1:8),TMAX,NUM！IS-基体的二维数组PRINT*,'PLEASEINPUTLOOP(TMAX)'READ*,TMAXISEED=RTC()IS=15！“死”的状态，基体为白色！赋予生命的种子，“活”的状态1DOI=1,NMAXIX0=IR*RAN(ISEED)1JY0=JR*RAN(ISEED)1  IS(IX0,JY0)=1ENDDO  !EXECUTETHERULEDOT=1,TMAXIS1=IS！边界条件IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX1,0)=IS(0:IMAX1,JMAX)IS(0:IMAX1,JMAX1)=IS(0:IMAX1,1)！搜索生命存在的位置  DOIX=1,IRDOJY=1,JR！判断邻居状态ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY1),IS(IX,JY-1)&    ,IS(IX,JY1),IS(IX1,JY-1),IS(IX1,JY),IS(IX1,JY1)/)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)=1  ELSEIS1(IX,JY)=15ENDIF！画图ISRE=SETCOLOR(IS1(IX,JY))ISRE=SETPIXEL(IX,JY)ENDDOENDDOIS=IS1ENDDOEND九）元胞自动机模拟单晶长大USEMSFLIBPARAMETERIR=400,JR=400INTEGERIS(0:IR1,0:JR1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY!! 根据过去状态IS，定义现在状态为IS1；为了突出边界，特别定义ISN1INTEGERIS1(0:IR1,0:JR1),ISN1(1:8)WRITE(*,*)"PLEASEINPUTTHETIMESTEP"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")DOT=1,TMAX!! 每次循环前重新定义过去状态数组ISIS=IS1!  边界条件IS(0,0:JR1)=IS(IR,0:JR1)IS(IR1,0:JR1)=IS(1,0:JR1)IS(0:IR1,0)=IS(0:IR1,JR)IS(0:IR1,JR1)=IS(0:IR1,1)!!   整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变DOIX=1,IRDOJY=1,JRISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY1),IS(IX,JY-1)&      ,IS(IX,JY1),IS(IX1,JY-1),IS(IX1,JY),IS(IX1,JY1)/)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-E0IG2.5*RD-1.25!!    用现在状态数组IS1记录边界状态的改变IF(DE.LT.0.0)IS1(IX,JY)=NSTATEENDDOENDDO!!  每循环20次在显示屏幕上刷新状态（颜色）DOIX=1,IRDOJY=1,JR!      IF(MOD(T,20).EQ.0)THENISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY1),IS1(IX,JY-1)&   ,IS1(IX,JY1),IS1(IX1,JY-1),IS1(IX1,JY),IS1(IX1,JY1)/)ISRE=SETCOLOR(IS(IX,JY))!    如果是边界，定义特别颜色15,加以区分        IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)!      ENDIFENDDOENDDO!!!  记录循环次数和对应的晶粒面积WRITE(1,*)T,SQRT(1.0*COUNT(IS.EQ.1))ENDDOCLOSE(1)END   十）元胞自动机模拟多晶长大1.介绍蒙特卡罗和元胞自动机的区别。⏹元胞自动机特点：时间、空间、状态都离散的动力学模型，⏹利用确定性或概论性的转换规则来描述在离散空间和时间上复杂系统演化USEMSFLIBPARAMETERIR=400,JR=400,NMAX=100INTEGERIS(0:IR1,0:JR1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JYINTEGERIGV(0:NMAX)!! 根据过去状态IS，定义现在状态为IS1；为了突出边界，特别定义ISN1INTEGERIS1(0:IR1,0:JR1),ISN1(1:8)WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXISEED=RTC()!  定义基体体积能为10,晶粒体积能为1IGV=1IGV(0)=10!  定义晶粒长大位置，并附能量，附带生长取向（研究形核数目与晶粒尺寸的关系）DOI=1,NMAX   IX0=IR*RAN(ISEED)1JY0=JR*RAN(ISEED)1  IF(IS(IX0,JY0).NE.0)CYCLEIS(IX0,JY0)=IENDDO  !! 在循环前定义现在状态数组IS1的初始值IS1=IS!!OPEN(1,FILE="E:\LUKE.DAT")DOT=1,TMAX!! 每次循环前重新定义过去状态数组ISIS=IS1!!!  边界条件IS(0,0:JR1)=IS(IR,0:JR1)IS(IR1,0:JR1)=IS(1,0:JR1)IS(0:IR1,0)=IS(0:IR1,JR)IS(0:IR1,JR1)=IS(0:IR1,1)!  寻找晶粒边界，计算能量，改变状态。  !!   整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变DOIX=1,IRDOJY=1,JR!!ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY1),IS(IX,JY-1)&      ,IS(IX,JY1),IS(IX1,JY-1),IS(IX1,JY),IS(IX1,JY1)/)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=IGV(NSTATE)-IGV(IS(IX,JY))    DE=E-E0IG2.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)ENDDOENDDO!!  每循环20次在显示屏幕上刷新状态（颜色）DOIX=1,IRDOJY=1,JRIF(MOD(T,20).EQ.0)THENISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY1),IS1(IX,JY-1)&   ,IS1(IX,JY1),IS1(IX1,JY-1),IS1(IX1,JY),IS1(IX1,JY1)/)ISRE=SETCOLOR(MOD(IS1(IX,JY),15))!    如果是边界，定义特别颜色15,加以区分        IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0)ISRE=SETCOLOR(15)ISRE=SETPIXEL(IX,JY)ENDIFENDDOENDDO!!!  记录循环次数和对应的晶粒面积WRITE(1,*)T,SQRT(1.0*COUNT(IS1.GT.0))ENDDOCLOSE(1)END   十一）  相场方法模拟相变BYUSINGBOUNDARYTRACKINGMETHORD,2008,11.23!  GRAINGROWTHVELOCITYIS "1.0*(IT-NT(IST))"OR"A*(IT-NT(IST))**B"!  WHICHSHOULDBENOTGREETERTHAN1CELLEACHSTEP.USEMSFLIBPARAMETERIMAX=800,JMAX=800,NMAX=50000,NSTEP=10!最多晶核数,每步形成晶核数INTEGERIS(0:IMAX1,0:JMAX1),IST(0:IMAX1,0:JMAX1)INTEGER IN(1:8),JN(1:8)INTEGERNI(1:NMAX),NJ(1:NMAX),NT(1:NMAX)!核的I,J坐标,形成时间INTEGERNNS,DCN,DCNM!运行时间部，最近晶核状态,访问CELL与晶核的距离平方，最短距离平方ISEED=RTC()WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXIN=(/-1,-1,-1,0,0,1,1,1/)JN=(/-1,0,1,-1,1,-1,0,1/)OPEN(5,FILE="C:\IRX.DAT")DOIT=1,TMAXIS=IST!  边界条件IS(0,1:JMAX)=IS(IMAX,1:JMAX)IS(IMAX1,1:JMAX)=IS(1,1:JMAX)IS(0:IMAX1,0)=IS(0:IMAX1,JMAX)IS(0:IMAX1,JMAX1)=IS(0:IMAX1,1)!--------------------------形核DOI=1,NSTEPINC=RAN(ISEED)*IMAX1JNC=RAN(ISEED)*JMAX1IF(IST(INC,JNC).NE.0)CYCLENN=NN1NI(NN)=INCNJ(NN)=JNCNT(NN)=ITIST(INC,JNC)=NN  !ENDDOC   -------------------DOI=1,IMAXDOJ=1,JMAXIF(IS(I,J).GT.0)CYCLEDCNM=2*IMAX*JMAXNP=0DOK=1,8II=IIN(K)JJ=JJN(K)IF(IS(II,JJ).EQ.0)CYCLE   !效率提高IF(NT(IS(II,JJ)).EQ.IT)CYCLE !当前深刻形成的核,当前时刻不长大NP=NP1           !只有IT时刻以前形成的核长大 IDN=I-NI(IS(II,JJ))JDN=J-NJ(IS(II,JJ))IDNABS=ABS(IDN)JDNABS=ABS(JDN)  DCN=IDN*IDNJDN*JDNIF(DCN.LE.DCNM)THENDCNM=DCN !访问CELL与最近晶核的距离平方NNS=IS(II,JJ)!ENDIFENDDOIF(NP.EQ.0)CYCLEGR=1.0*(IT-NT(NNS))!晶粒半径DIS=GR-SQRT(1.0*DCNM)IF(DIS.GE.0.0)IST(I,J)=NNSIII=MOD(IST(I,J),15)IF(III==0)III=III1IER=SETCOLOR(III)IER=SETPIXEL(I,J)ENDDOENDDO!  计算相变分数，IRX=COUNT(IST.NE.0)WRITE(5,*)IT,1.0*IRX/(IMAX*JMAX)ENDDOCLOSE(5)  END!  相变byUsingBoundaryTrackingMethord,2008,11.13!  graingrowthvelocityis "1.0*(IT-NT(IST))"or"a*(IT-NT(IST))**b"!  whichshouldbenotgreeterthan1celleachstep.USEMSFLIBPARAMETERIMAX=400,JMAX=400INTEGERIST(0:IMAX1,0:JMAX1)INTEGERDCNM !访问CELL与晶核的距离平方ISEED=RTC()WRITE(*,*)"PLEASEINPUTTHETIMESTEP"READ(*,*)TMAXOPEN(5,FILE="C:\IRX.DAT")IRC=IMAX/2JRC=JMAX/2IST=1IST(IRC,JRC)=2DOIT=1,TMAX!  边界条件IST(0,1:JMAX)=IST(IMAX,1:JMAX)IST(IMAX1,1:JMAX)=IST(1,1:JMAX)IST(0:IMAX1,0)=IST(0:IMAX1,JMAX)IST(0:IMAX1,JMAX1)=IST(0:IMAX1,1)C   -------------------DOI=1,IMAXDOJ=1,JMAXIDN=I-IRCJDN=J-JRCIDNABS=ABS(IDN)