fortran程序汇总.docx
- 文档编号:6184565
- 上传时间:2023-01-04
- 格式:DOCX
- 页数:21
- 大小:1.99MB
fortran程序汇总.docx
《fortran程序汇总.docx》由会员分享,可在线阅读,更多相关《fortran程序汇总.docx(21页珍藏版)》请在冰豆网上搜索。
fortran程序汇总
计算圆周率
REALR,R1,R2,PI
ISEED=RTC()
N0=0
N=300000
DOI=1,N
R1=RAN(ISEED)
R2=RAN(ISEED)
R=SQRT(R1*R1+R2*R2)
IF(R<1.0)N0=N0+1
ENDDO
PI=4.0*N0/N
WRITE(*,*)PI
END
一)蒙特卡洛计算生日问题
假设有N个人在一起,各自的生日为365天之一,根据概率理论,与很多人的直觉相反,只需23个人便有大于50%的几率人群中至少有2个人生日相同。
INTEGERM(1:
10000),NUMBER1(0:
364),NUMBER2
REALX,Y
ISEED=RTC()
DOJ=1,10000
NUMBER1=0
X=RAN(ISEED)
NUMBER1(0)=INT(365*X+1)
JJJ=1
DOI=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
ENDIF
ENDDO
ENDDO
DOI=1,10000
IF(M(I).LE.23)SUM=SUM+1
ENDDO
PRINT*,SUM/10000
END
二)MONTECARLOSIMULATIONOFONEDIMENSIONALDIFFUSION蒙特卡罗计算一维扩散问题
INTEGERX,XX(1:
1000,1:
1000)
REALXXM(1:
1000)
!
X:
INSTANTANEOUSPOSITIONOFATOM
!
XX(J,I):
X*X,J:
第几天实验,I:
第几步跳跃
!
XXM(I):
THEMEANOFXX
WRITE(*,*)"实验天数JMAX,实验次数IMAX"
READ(*,*)JMAX,IMAX
ISEED=RTC()
DOJ=1,JMAX!
第几天实验
X=0!
!
!
DOI=1,IMAX!
第几步跳跃
RN=RAN(ISEED)
IF(RN<0.5)THEN
X=X+1
ELSE
X=X-1
ENDIF
XX(J,I)=X*X
ENDDO
ENDDO
OPEN(1,FILE="C:
\DIF1.DAT")
DOI=1,IMAX
XXM=0.0
XXM(I)=1.0*SUM(XX(1:
JMAX,I))/JMAX!
!
WRITE(1,*)I,XXM(I)
ENDDO
CLOSE
(1)
END
三维的!
三)通过该程序了解FORTRAN语言如何画图(通过像素画图)
USEMSFLIB
INTEGERXR,YR!
在的区域中画一个圆
PARAMETERXR=400,YR=400
INTEGERR,S(1:
XR,1:
YR)
X0=XR/2!
圆心位置X0,YO
Y0=YR/2
R=MIN(X0-10,Y0-10)!
圆半径
S=0!
像素的初始状态(颜色)
DOI=1,XR
DOJ=1,YR
IF((I-X0)**2+(J-Y0)**2<=R**2)S(I,J)=10
IER=SETCOLOR(S(I,J))
IER=SETPIXEL(I,J)
ENDDO
ENDDO
END
四)画一个圆(1、如何选出晶界区域;2、进一步加深对画图的理解)
USEMSFLIB
INTEGERXR,YR!
在的区域中画一个圆
PARAMETERXR=400,YR=400
INTEGERR,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,Y0
Y0=YR/2
R=MIN(X0-10,Y0-10)!
圆半径
S=0!
像素的初始状态(颜色)
DOI=1,XR
DOJ=1,YR
IF((I-X0)**2+(J-Y0)**2<=R**2)S(I,J)=10
IER=SETCOLOR(S(I,J))
IER=SETPIXEL(I,J)
ENDDO
ENDDO
DOI=1,XR!
画晶界
DOJ=1,YR
NDS=0
DOK=1,4
IF(S(I,J).NE.S(I+XN(K),J+YN(K)))NDS=NDS+1
ENDDO
IF(NDS>0)THEN
IER=SETCOLOR(9)
ELSE
IER=SETCOLOR(8)
ENDIF
IER=SETPIXEL(I,J)
ENDDO
ENDDO
END
五)MC模拟一个晶粒的缩小
USEMSFLIB
PARAMETERIR=400,JR=400
INTEGERIS(0:
IR+1,0:
JR+1),TMAX,ISN(1:
8),NSTATE,T,NR,IX,IY
WRITE(*,*)"PLEASEINPUTTHETIMESTEP"
READ(*,*)TMAX
ISEED=RTC()
!
定义圆心和半径
IRC=IR/2
JRC=IR/2
R=MIN(IRC,JRC)-10
!
定义基体和圆晶粒分别为状态1、状态2
IS=1
DOI=1,IR
DOJ=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)
ENDDO
ENDDO
OPEN(1,FILE="E:
\LUKE.DAT")
!
寻找晶粒边界,计算能量,改变状态。
DOT=1,TMAX
DOX=1,IR
DOY=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
六)蒙特卡罗模拟基体上单晶粒形核长大
USEMSFLIB
PARAMETERIR=400,JR=400
INTEGERIS(0:
IR+1,0:
JR+1),TMAX,ISN(1:
8),NSTATE,T,NR,IX,IY
WRITE(*,*)"PLEASEINPUTTHETIMESTEP"
READ(*,*)TMAX
ISEED=RTC()
!
定义圆心和半径
IRC=IR/2
JRC=IR/2
!
定义基体和圆晶粒分别为状态10、状态2
IS=10
IS(IRC,JRC)=2
OPEN(1,FILE="E:
\LUKE.DAT")
!
寻找晶粒边界,计算能量,改变状态。
DOT=1,TMAX
DOX=1,IR
DOY=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))
ENDDO
CLOSE
(1)
END
七)蒙特卡洛模拟基体上多晶粒形核长大,
USEMSFLIB
!
定义常数名
PARAMETERIR=400,JR=400,NMAX=100
!
定义变量:
体积能变量(一维数组);基体各像素点变量(二维数组),邻居的取向号(一维数组)等;
INTEGERIS(0:
IR+1,0:
JR+1),TMAX,ISN(1:
8),NSTATE,T,NR,IX0,IY0,IX,IY
INTEGERIGV(0:
NMAX)
!
读取晶粒生长步长;
WRITE(*,*)"PLEASEINPUTTHETIMESTEP"
READ(*,*)TMAX
ISEED=RTC()
!
定义基体体积能为10,所有晶粒体积能为1
IGV=1
IGV(0)=10
!
定义晶粒长大位置(100个形核点初始形核位置),并附已不同的取向号
DOI=1,NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IS(IX0,JY0)=I
ENDDO
!
边界条件
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")
!
寻找晶粒边界。
DOT=1,TMAX
DOX=1,IR
DOY=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)
ENDDO
CLOSE
(1)
END
八)元胞自动机模拟生命游戏程序(生命永不停止)
USEMSFLIB
PARAMETERIR=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*,TMAX
ISEED=RTC()
IS=15!
“死”的状态,基体为白色
!
赋予生命的种子,“活”的状态1
DOI=1,NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IS(IX0,JY0)=1
ENDDO
!
EXECUTETHERULE
DOT=1,TMAX
IS1=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)
!
搜索生命存在的位置
DOIX=1,IR
DOJY=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
ENDIF
!
画图
ISRE=SETCOLOR(IS1(IX,JY))
ISRE=SETPIXEL(IX,JY)
ENDDO
ENDDO
IS=IS1
ENDDO
END
九)元胞自动机模拟单晶长大
USEMSFLIB
PARAMETERIR=400,JR=400
INTEGERIS(0:
IR+1,0:
JR+1),TMAX,ISN(1:
8),NSTATE,T,NR,IX0,IY0,IX,JY
!
!
根据过去状态IS,定义现在状态为IS1;为了突出边界,特别定义ISN1
INTEGERIS1(0:
IR+1,0:
JR+1),ISN1(1:
8)
WRITE(*,*)"PLEASEINPUTTHETIMESTEP"
READ(*,*)TMAX
ISEED=RTC()
IRC=IR/2!
=IR*ran(iseed)+1
JRC=JR/2
!
定义基体体积能为10,晶粒体积能为1
IS=8
IS(IRC,JRC)=1
!
!
在循环前定义现在状态数组IS1的初始值
IS1=IS
OPEN(1,FILE="E:
\LUKE.DAT")
DOT=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)
!
!
整个基体上,各个点上的状态都要根据规则改变,而非随机取点改变
DOIX=1,IR
DOJY=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
ENDDO
!
!
每循环20次在显示屏幕上刷新状态(颜色)
DOIX=1,IR
DOJY=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)
!
ENDIF
ENDDO
ENDDO
!
!
!
记录循环次数和对应的晶粒面积
WRITE(1,*)T,SQRT(1.0*COUNT(IS.EQ.1))
ENDDO
CLOSE
(1)
END
十)元胞自动机模拟多晶长大
1.介绍蒙特卡罗和元胞自动机的区别。
⏹元胞自动机特点:
时间、空间、状态都离散的动力学模型,
⏹利用确定性或概论性的转换规则来描述在离散空间和时间上复杂系统演化
USEMSFLIB
PARAMETERIR=400,JR=400,NMAX=100
INTEGERIS(0:
IR+1,0:
JR+1),TMAX,ISN(1:
8),NSTATE,T,NR,IX0,IY0,IX,JY
INTEGERIGV(0:
NMAX)
!
!
根据过去状态IS,定义现在状态为IS1;为了突出边界,特别定义ISN1
INTEGERIS1(0:
IR+1,0:
JR+1),ISN1(1:
8)
WRITE(*,*)"PLEASEINPUTTHETIMESTEP"
READ(*,*)TMAX
ISEED=RTC()
!
定义基体体积能为10,晶粒体积能为1
IGV=1
IGV(0)=10
!
定义晶粒长大位置,并附能量,附带生长取向(研究形核数目与晶粒尺寸的关系)
DOI=1,NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IF(IS(IX0,JY0).NE.0)CYCLE
IS(IX0,JY0)=I
ENDDO
!
!
在循环前定义现在状态数组IS1的初始值
IS1=IS
!
!
OPEN(1,FILE="E:
\LUKE.DAT")
DOT=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)
!
寻找晶粒边界,计算能量,改变状态。
!
!
整个基体上,各个点上的状态都要根据规则改变,而非随机取点改变
DOIX=1,IR
DOJY=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次在显示屏幕上刷新状态(颜色)
DOIX=1,IR
DOJY=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)
ENDIF
ENDDO
ENDDO
!
!
!
记录循环次数和对应的晶粒面积
WRITE(1,*)T,SQRT(1.0*COUNT(IS1.GT.0))
ENDDO
CLOSE
(1)
END
十一)相场方法模拟相变BYUSINGBOUNDARYTRACKINGMETHO
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- fortran 程序 汇总