计算材料学Fortran程序汇总
INTEGER M(1:10000), NUMBER1(0:364), NUMBER2
REAL X,Y
ISEED=RTC()
DO J=1, 10000
NUMBER1=0
X=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 DO
END DO
DO I=1,10000
IF(M(I).LE.23) SUM=SUM+1
END 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，实验次数IMAX"
READ(*,*) JMAX,IMAX
ISEED=RTC()
DO J=1,JMAX !第几天实验
X=0 !!!
DO I=1,IMAX !第几步跳跃
RN=RAN(ISEED)
IF(RN<0.5)THEN
X=X+1
ELSE
X=X-1
END IF
XX(J,I)=X*X
END DO
END DO
OPEN(1,FILE="C:\DIF1.DAT")
DO I=1,IMAX
XXM=0.0
XXM(I)=1.0*SUM(XX(1:JMAX,I))/JMAX !!
WRITE(1,*) I, XXM(I)
END DO
CLOSE(1)
END
============================================================ ! Monte Carlo Simulation of Two Dimensional Diffusion
INTEGER X,Y,XY(1:1000,1:1000)
REAL XYM(1:1000)
! X:INSTANTANEOUS POSITION OF ATOM
! XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃
! XYM(I): THE MEAN OF XY
WRITE(*,*) "实验天数JMAX，实验次数IMAX"
READ(*,*) JMAX,IMAX
ISEED=RTC()
DO J=1,JMAX !第几天实验
X=0 !!!
Y=0 !!!
DO I=1,IMAX !第几步跳跃
RN=RAN(ISEED)
IF(RN.LT.0.25)THEN
x=x
y=y-1
END IF
IF(RN.LT.0.5.AND.RN.GE.0.25)THEN
x=x
y=y+1
END IF
IF(RN.LT.0.75.AND.RN.GE.0.5)THEN
x=x-1
y=y
END IF
IF(RN.GE.0.75)THEN
x=x+1
y=y
END IF
XY(J,I)=X*X+Y*Y
END DO
END DO
OPEN(1,FILE="C:\DIF2.DAT")
DO I=1,IMAX
XYM=0.0
XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX !!
WRITE(1,*) I, XYM(I)
END DO
CLOSE(1)
END
===============================================
! Monte Carlo Simulation of One Dimensional Diffusion INTEGER X,XY(1:1000,1:1000),y,XN(1:4),YN(1:4),RN
REAL XYM(1:1000)
! X:INSTANTANEOUS POSITION OF ATOM
! XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃
! XYM(I): THE MEAN OF XY
WRITE(*,*) "实验天数JMAX，实验次数IMAX"
READ(*,*) JMAX,IMAX
XN=(/0,0,-1,1/)
YN=(/-1,1,0,0/)
ISEED=RTC()
DO J=1,JMAX !第几天实验
X=0 !!!
Y=0 !!!
DO I=1,IMAX !第几步跳跃
RN=4*RAN(ISEED)+1
X=X+YN(RN)
Y=Y+YN(RN)
XY(J,I)=X*X+Y*Y
END DO
END DO
OPEN(1,FILE="C:\DIF2.DAT")
DO I=1,IMAX
XYM=0.0
XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX !!
WRITE(1,*) I, XYM(I)
END DO
CLOSE(1)
END
做三维空间随机行走？？留作业
! Monte Carlo Simulation of One Dimensional Diffusion INTEGER X,XY(1:1000,1:1000),y,XN(1:6),YN(1:6),ZN(1:6),RN REAL XYM(1:1000)
! X:INSTANTANEOUS POSITION OF ATOM
! XY(J,I)：X*Y ,J:第几天实验,I:第几步跳跃
! XYM(I): THE MEAN OF XY
WRITE(*,*) "实验天数JMAX，实验次数IMAX"
READ(*,*) JMAX,IMAX
XN=(/0,0,-1,1,0,0/)
YN=(/-1,1,0,0,0,0/)
ZN=(/0,0,0,0,1,-1/)
ISEED=RTC()
DO J=1,JMAX !第几天实验
X=0 !!!
Y=0 !!!
Z=0
DO I=1,IMAX !第几步跳跃
RN=6*RAN(ISEED)+1
X=X+XN(RN)
Y=Y+YN(RN)
Z=Z+ZN(RN)
XY(J,I)=X*X+Y*Y+Z*Z
END DO
END DO
OPEN(1,FILE="C:\DIF2.DAT")
DO I=1,IMAX
XYM=0.0
XYM(I)=1.0*SUM(XY(1:JMAX,I))/JMAX !!
WRITE(1,*) I, XYM(I)
END DO
CLOSE(1)
==============================================
通过该程序了解fortran语言如何画图（通过像素画图）
USE MSFLIB
INTEGER XR，YR !在的区域中画一个圆
PARAMETER XR=400,YR=400
INTEGER R, S(1:XR,1:YR)
X0=XR/2 ! 圆心位置X0，YO
Y0=YR/2
R=MIN(X0-10,Y0-10) !圆半径
S=0 !像素的初始状态（颜色）
DO I=1,XR
DO J=1,YR
IF((I-X0)**2+(J-Y0)**2<=R**2)S(I,J)=10
IER=SETCOLOR(S(I,J))
IER=SETPIXEL(I,J)
END DO
END DO
END
=================================================
! 画一个圆（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，Y0
Y0=YR/2
R=MIN(X0-10,Y0-10) !圆半径
S=0 !像素的初始状态（颜色）
DO I=1,XR
DO J=1,YR
IF((I-X0)**2+(J-Y0)**2<=R**2)S(I,J)=10
IER=SETCOLOR(S(I,J))
IER=SETPIXEL(I,J)
END DO
END DO
DO I=1,XR !画晶界
DO J=1,YR
NDS=0
DO K=1,4
IF(S(I,J).NE.S(I+XN(K),J+YN(K)))NDS=NDS+1
END DO
IF(NDS>0)THEN
IER=SETCOLOR(9)
ELSE
IER=SETCOLOR(8)
END IF
IER=SETPIXEL(I,J)
END DO
END DO
END
如何画有一定宽度的晶界？
====================================================
！MC模拟一个晶粒的缩小
（1、如何定义基体和晶粒；
（2、如何寻找边界；
（3、如何计算能量（构造或描述问题的概率过程）
步骤：1、在基体上画一个原始晶粒，并赋予基体和原始晶粒不同的状态（取向号或体积能）
2、寻找晶界
3、MC的一个时间步（晶粒长大过程中一种随机性）
4、计算晶界上网格点的相互作用，每个格点的相互作用导致晶粒长大
①随机选取一个初始格点；
②若此点属于晶界，那么可以随机转变为它邻居的状态，
若不是晶界，则跳出循环，不发生转变；
③计算转变前后的能量变化，⊿E=⊿E v +⊿E s +⊿E q
(自由能=体积能+表面能+能量起伏)
④若⊿E小于0，则新晶相被接受，网格状态发生转变。
USE MSFLIB
PARAMETER IR=400,JR=400
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTA TE,T,NR,IX,IY
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
! 定义圆心和半径
JRC=JR/2
R=MIN(IRC,JRC)-10
! 定义基体和圆晶粒分别为状态1、状态2
IS=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 DO
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)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTA TE)
RD=RAN(ISEED)
DE=E-E0+NSTA TE-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
（1、演示体积能互换的转变情况；
2、演示IX=IR*RAN(ISEED)+1
JY=JR*RAN(ISEED)+1的替换；
3、能量关系变化；
）
============================================
作业3：基体上单晶粒形核长大
（1、选取初始形核点
（2、定义体积能的变化）
USE MSFLIB
PARAMETER IR=400,JR=400
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTA TE,T,NR,IX,IY WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
! 定义圆心和半径
IRC=IR/2
JRC=IR/2
! 定义基体和圆晶粒分别为状态10、状态2
IS=10
IS(IRC,JRC)=2
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)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTA TE)
RD=RAN(ISEED)
DE=E-E0+NSTA TE-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
基体上多晶粒形核长大，
（1、如何生成多晶粒形核核心
（2、如何定义、并区别多晶粒
3、为什么要有边界条件）
步骤：1、定义变量：体积能变量（一维数组）；基体各像素点变量（二维数组），邻居的取向号（一维数组）等；
2、读取晶粒生长步长；
3、寻找晶界
4、计算晶界上网格点的相互作用，每个格点的相互作用
导致晶粒长大或缩小：
①随机选取一个初始格点；
②若此点属于晶界，那么可以随机转变为它邻居的状态，
若不是晶界，则跳出循环，不发生转变；
③计算转变前后的能量变化，⊿E=⊿E v +⊿E s +⊿E q
(自由能=体积能+表面能+能量起伏)
④若⊿E小于0，则新晶相被接受，网格状态发生转变。
=============================================================
USE MSFLIB
! 定义常数名
PARAMETER IR=400,JR=400,NMAX=100
! 定义变量：体积能变量（一维数组）；基体各像素点变量（二维数组），邻居的取向号（一维数组）等；
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTA TE,T,NR,IX0,IY0,IX,IY
INTEGER IGV(0:NMAX)
! 读取晶粒生长步长；
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
! 定义基体体积能为10,所有晶粒体积能为1
IGV=1
IGV(0)=10
! 定义晶粒长大位置（100个形核点初始形核位置），并附已不同的取向号
DO I=1, NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IS(IX0,JY0)=I
END 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.NSTA TE)
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
! iii= MOD(IS(IX,JY),15)
! If(iii==0) iii=iii+1物理意义
1晶粒长大与时间的关系
2 形核数目与的关系
3、与EBSD 比较
4、为什么圆弧会长成近似的直线
====================================================================
Monte-Carlo 方法模拟晶粒长大过程中的新研究进展
第一个问题：采用Monte-Carlo 方法存在以下几方面的问题
（1）模拟过程只计及最近邻格点的相互作用，
难于考虑畸变场等物理场以及晶界形态对晶粒长大过程的影响，不利于在唯象层次上深入开展对晶粒
长大过程的研究。
应当指出，大部分的计算机模拟主要集中于二维系统，虽然能够在一定程度上提供晶粒长大的有关动力学和拓扑学的信息，但是，任何多晶体材料中的晶粒都是三维的，所有二维晶粒长大的模拟都是基于某种特定的假设或简化，因此模拟三维晶粒的长大更具有实际意义。同时，在目前的研究中，对晶粒生长过程的形貌演化模拟的逼真度比较好，
但在生长模型的边界处理等方面还有待完善，特别是晶粒生长与晶化温度、时间、气氛等参数密切相关，这种复杂工艺条件下的晶粒生长过程的模拟是当前该领域正待解决的难题。另外，国内外用各种方法模拟晶粒长大，主要集中在对方法本身的讨论、模拟晶粒分布函数特点及晶粒形貌等方面，而将模拟方法与实际工艺联系起来的研究者比较少。
第二个问题：算法种类：
1、标准算法
2、改进型算法
3、其它算法
第三个问题：模拟晶粒长大的一些应用
1、金属凝固过程的形核；
2、再结晶及存在织构的再结晶（形核位置时考虑能量关系）；
3、第二相粒子的影响（第二相粒子对晶粒长大有阻碍作用. 这种情况下模拟较复杂. 由于无法用解析解定量描述, 用MC 方法模拟显得
尤为重要；
5、异常晶粒的长大；
6、三维晶粒的长大（模拟不同的点阵结构生长过程）；
山东大学中国有色金属学报，2008，4
7、蒙特卡洛方法不是数学的方法，而是一种多次模拟实验的方法。
===========================================================
元胞自动机程序（生命永不停止）
use msflib
PARAMETER 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*, Tmax
ISEED=RTC()
is=15 ！“死”的状态，基体为白色
！赋予生命的种子，“活”的状态1
DO I=1, NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IS(IX0,JY0)=1
END DO
!execute the rule
do t=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)
！搜索生命存在的位置
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)/)
num=COUNT(ISN.eq.1)
！赋予生存的条件
if((is(ix,jy)==15.and.num==3).or.(is(ix,jy)==1.and.(num==3.or.num== 4))) then is1(ix,jy)=1
else
is1(ix,jy)=15
end if
！画图
ISRE=SETCOLOR(IS1(IX,JY))
ISRE=SETPIXEL(IX,JY)
end do
end do
is=is1
end do
end
======================================================
元胞自动机—产生和发展. 四个阶段:
1940s 诞生：Von Neumann 自我复制机.
1960-70s起步：JH.Conway 生命游戏.
1980s 理论研究：S.Wolfram CA分类.
1980-90s 应用：HPP-FHP格子气自动机.
https://www.360docs.net/doc/6413463407.html,ngton N.Packard 人工生命
元胞自动机（Cellular Automata，简称CA，也有人称其为细胞自动机、
点格自动机、分子自动机或单元自动机）是一种建立在离散的时间和
空间上的动力学系统。散布在规则格网（Lattice Grid）中的每一元胞
（Cell）取有限的离散状态，遵循同样的作用规则，依据确定的局部
规则作同步更新。大量元胞通过简单的相互作用而构成动态系统的演化。
与一般的动力学模型不同，元胞自动机不是由严格定义的物理方程
或函数确定，而是由一系列模型构造的规则构成，凡是满足这些规则的模
型都可以算做是元胞自动机模型。因此，元胞自动机是一类模型的总称，或
者是一个方法框架。其特点是时间、空间、状态都是离散的，每个变量只
取有限多个状态，且其状态改变的规则在时间和空间上都是局部的。
元胞自动机是一种对具有局域连通性的格点，应用局部（有时为中等
范围）确定性或概论性的转换规则来描述在离散空间和时间上复杂系统演化
规律的同步算法。
元胞自动机可以采取任意的维数，自动机的格点描述被认为与所选模型
密切相关的基本的系统实体。独立的格点可以表示连续的体积单元，原子颗粒
大原子团、汽车、人口、晶格缺陷等。
通过作用于每个格点的确定性或概率性转换规则，自动机的演化得以实施。这些规则决定了每个元胞的状态，并将其描述为其前一状态和邻近元胞的函数。在
计算状态转变过程中所考虑的邻近位置的数目确定了相互作用的范围和自动机
的局部演化。
！转变概率
元胞自动机特点：时间、空间、状态都离散的动力学模型，
利用确定性或概论性的转换规则来描述在离散空间和时间上复杂系统演化
============================================================= 元胞自动机模拟单晶长大
USE MSFLIB
PARAMETER IR=400,JR=400
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTA TE,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(*,*)TMAX
ISEED=RTC()
IRC=IR/2
JRC=JR/2
!!
! 定义基体体积能为10,晶粒体积能为1
IS=10
IS(IRC,JRC)=1
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
! 边界条件
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,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.NSTA TE)
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)=NSTA TE
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))
ENDDO
CLOSE(1)
END
===================================================
元胞自动机模拟单晶长大，再附加一个规则
USE MSFLIB
PARAMETER IR=400,JR=400,NMAX=100
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTA TE,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(*,*)TMAX
ISEED=RTC()
IRC=IR/2
JRC=JR/2
! 定义基体体积能为10,晶粒体积能为1
IS=10
IS(IRC,JRC)=1
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
!!
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
!!
! 边界条件
IS(0,1:JR)=IS(IR,1:JR)
IS(IR,1:JR)=IS(iR+1,1:JR)
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
! 随机寻找一个相邻点
sss=0
do iii=1,8
NSTATE=ISN(iii)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTA TE)
RD=RAN(ISEED)
IG=NSTATE-IS(IX,JY)
DE=E-E0+IG+2.5*RD-1.25
!! 用现在状态数组IS1记录边界状态的改变
IF(DE.LT.0.0) then
sss=sss+1
end if
end do
if(sss.gt.5) IS1(IX,JY)=11-is(ix,jy)
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(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))
ENDDO
CLOSE(1)
END
USE MSFLIB
PARAMETER IR=400,JR=400,NMAX=100
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTA TE,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(*,*)TMAX
ISEED=RTC()
! 定义基体体积能为10,晶粒体积能为1
IGV=1
! 定义晶粒长大位置，并附能量
DO I=1,NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IF(IS(IX0,JY0).NE.0)CYCLE
IS(IX0,JY0)=I
END DO
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
!!
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
!!
! 边界条件
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,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.NSTA TE)
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)=NSTA TE
!!
!! ISRE=SETCOLOR(MOD(IS(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(IS(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))
ENDDO
CLOSE(1)
END
错误的多晶长大
USE MSFLIB
PARAMETER IR=400,JR=400,NMAX=100
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTA TE,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(*,*)TMAX
ISEED=RTC()
!!
! 定义基体体积能为10,晶粒体积能为1
IGV=1
IGV(0)=10
! 定义晶粒长大位置，并附能量
DO I=1,NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IS(IX0,JY0)=I
END DO
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
!!
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
!!
! 边界条件
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,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
! 随机寻找一个相邻点
sss=0
do iii=1,8
! NR=8*RAN(ISEED)+1
NSTATE=ISN(iii)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTA TE)
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
IF(DE.LT.0.0) then
sss=sss+1
nstate1=nstate
if(sss.gt.5) IS1(IX,JY)=nstate1
end if
end do
ENDDO
ENDDO
!! 每循环20次在显示屏幕上刷新状态（颜色）
DO IX=1,IR
DO JY=1,JR
IF(MOD(T,40).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(IS(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, COUNT(IS1.GT.0)
ENDDO
CLOSE(1)
END
USE MSFLIB
PARAMETER IR=400,JR=400,NMAX=100
INTEGER IS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTA TE,T,NR,IX0,IY0,IX,JY INTEGER IGV(0:NMAX)
!! 根据过去状态IS，定义现在状态为IS1；为了突出边界，特别定义ISN1 INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8),ISN2(1:8)
WRITE(*,*)"PLEASE INPUT THE TIME STEP "
READ(*,*)TMAX
ISEED=RTC()
!!
! 定义基体体积能为10,晶粒体积能为1
IGV=1
IGV(0)=10
! 定义晶粒长大位置，并附能量
DO I=1,NMAX
IX0=IR*RAN(ISEED)+1
JY0=JR*RAN(ISEED)+1
IS(IX0,JY0)=I
END DO
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
!!
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
!!
! 边界条件
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,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
! 随机寻找一个相邻点
iiii=0
do iii=1,8
! NR=8*RAN(ISEED)+1
NSTATE=ISN(iii)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTA TE)
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
IF(DE.LT.0.0) then
iiii=iiii+1
isn2(iiii)=nstate
id1=iiii*ran(iseed)+1
if(iiii.gt.5) IS1(IX,JY)=isn2(id1)
else
end if
end do
ENDDO
ENDDO
!! 每循环20次在显示屏幕上刷新状态（颜色）
DO IX=1,IR
DO JY=1,JR
IF(MOD(T,40).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, COUNT(IS1.GT.0)
ENDDO
CLOSE(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.
步骤：
1、程序的主题语句；
2、邻居关系的不同定义方法；
3、形核，每个时间步形核数10；
4、获得转变点与晶粒核心的距离；
5、比较与相场公式所得距离之差；
USE MSFLIB
PARAMETER 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 STEP"
READ(*,*)TMAX
IN=(/-1,-1,-1,0,0,1,1,1/)
JN=(/-1,0,1,-1,1,-1,0,1/)
OPEN(5,FILE="C:\IRX.DA T")
DO IT=1,TMAX
IS=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,NSTEP
INC=RAN(ISEED)*IMAX+1
JNC=RAN(ISEED)*JMAX+1
IF(IST(INC,JNC).NE.0)CYCLE
NN=NN+1
NI(NN)=INC
NJ(NN)=JNC
NT(NN)=IT
IST(INC,JNC)=NN !
END DO
C -------------------
DO I=1,IMAX
DO J=1,JMAX
IF(IS(I,J).GT.0)CYCLE
DCNM=2*IMAX*JMAX
NP=0
DO K=1,8
II=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 IF
END DO
IF(NP.EQ.0)CYCLE
GR=1.0*(IT-NT(NNS)) !晶粒半径
DIS=GR-SQRT(1.0*DCNM)
IF(DIS.GE.0.0)IST(i,j)=NNS
III= MOD(IST(I,J),15)
If(III==0) III=III+1
IER=SETCOLOR(III)
IER=SETPIXEL(I,J)
END DO
END DO
! 计算相变分数，
IRX=COUNT(IST.NE.0)
WRITE(5,*)IT,1.0*IRX/(IMAX*JMAX)
END DO
CLOSE(5)
END
单晶
! 相变by Using Boundary Tracking Methord,2008,11.13
! 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 MSFLIB
PARAMETER IMAX=400,JMAX=400
INTEGER IST(0:IMAX+1,0:JMAX+1)
INTEGER DCNM !访问CELL 与晶核的距离平方
ISEED=RTC()
WRITE(*,*)"PLEASE INPUT THE TIME STEP"
READ(*,*)TMAX
OPEN(5,FILE="C:\IRX.DA T")
IRC=IMAX/2
JRC=JMAX/2
IST=1
IST(IRC,JRC)=2
DO IT=1,TMAX
! 边界条件
IST(0,1:JMAX)=IST(IMAX,1:JMAX)
IST(IMAX+1,1:JMAX)=IST(1,1:JMAX)
IST(0:IMAX+1,0)=IST(0:IMAX+1,JMAX)
IST(0:IMAX+1,JMAX+1)=IST(0:IMAX+1,1)
FORTRAN 90 程序编程规范
FORTRAN 90 程序编程规范 Fortran 90 编程规范，使程序代码高度组织化，更加易读、易懂、易于维护，程序更加高效。使编出的程序更易懂、易于维护。 1 语言选择 数值预报创新系统软件开发应避免使用Fortran77 的某些过时特征以Fortran 90不一致的特征。选择Fortran 90 作为开发语言，并采用Fortran 90 的新功能，如动态内存的分配(dynamic memory allocation)、递归（recursion ）, 模块(modules)、POINTER 、长变量名、自由格式等。 Fortran 77其中某些只是一些冗余的功能，这些功能已经过时，另外，还有一些在Fortran90 中被证明是不好的用法，建议不要使用。 2 Fortran 90 的新特性 2.1.1 建议使用的Fortran 90 新特性 建议使用Fortran 90 提供的模块（module ），并用Use ONLY 指定module 中哪些变量或派生类型定义可用于调用程序。 尽量使用数组下标三元组，这样可优化并减少所需的代码行数。为提高可读性，要在括号内表明数组的维数，例如： 1dArrayA(:) = 1dArrayB(:) + 1dArrayC(:) 2dArray(: , :) = scalar * Another2dArray(: , :) 当访问数组的子集时，例如在有限差分等式中，可以通过使用下标三元组实现。例如：2dArray(: , 2:len2) = scalar *（ & Another2dArray(:, 1:len2 -1) & - Another2dArray(:, 2:len2) & ) 对程序单元（program units ）命名，并使用End program ，End subroutine ，End interface ，End module 等结构再次指定“program unit ”的名称。 在逻辑表达式中使用>、 >=、 ==、 <、 <=、 /=，它们分别代 替.gt.、.ge.、.eq.、.lt.、.le.、.ne. 。新的表示方法更接近标准的数学符号 在变量定义中始终使用“：：”；始终用“DIMENSION ”定义数组形状；始终用（len=）的语法格式声明字符变量的长度。
Fortran复习题讲解
FORTRAN程序设计复习题 选择题 B (1)下列各FORTRAN表达式中合法的是 A) S+T*2P >= 0.0 B) .NOT. (A*B+C) C) A2+B2/(C+D) <= 0.0 D) (A+B).NOT．A*B．GT．0.0(.NOT.只跟一个表达式) C (2)数学式(3/5)ex+y的FORTRAN表达式是 A) 3*EXP(X+Y)/5 B) 3*E* *(X+Y)/5.0 C) (3/5)*EXP(X+Y) D) (3.0/5.0)EXP(X+Y) D (3)下列FORTRAN77表达式中不合法的是 A) A．GT．B．EQV．C．GT．D B) A．AND．B．AND．C．AND．D C) .NOT．(X．LE．0.0) D) A．LT．B．LT．C．LT．D D (4)下列叙述中不正确的是 A) FORTRAN子程序可以单独编译 B) 对一个FORTRAN源程序进行编译和连接无误后可生成可执行文件 C) 即使编译和连接都正确无误，FORTRAN程序运行时仍可能出错 D) FORTRAN连接的主要任务是把函数库中的函数翻译成机器指令(主要任务为连接目标文件) B (5)在下列FORTRAN77运算符中，优先级最高的是 A) .AND. B) .NOT. C) .OR. D) .EQ. B (6)FORTRAN表达式"6/5+9/2**3/2"的值为 A) 33 B) 1 C) 5 D) 3 A (7)下列FORTRAN77 表达式中，合法的是：A）相当于0.0.LE.X .AND.（.NOT.（X.GE.10.0 ））A) 0.0.LE.X .AND..NOT.X.GE.10.0 B) 0.0.GE.X.LT.10.0 C) 0.0.LE.X.GE.10. D) 10.0.LE.X.OR..AND.X.GT.0.0 D (8)关于编译一个FORTRAN源程序文件，下列说法中错误的是 A) 允许编译只有一个主程序而没有子程序的源文件 B) 允许编译有多个子程序的源文件 C) 允许编译只有一个子程序而没有主程序的源文件 D) 允许编译有多个主程序的源文件 C (9)在FORTRAN77源程序中，续行标志符必须放在 A) 第1列 B) 第1-6列 C) 第6列 D) 第5列 D (10)下列关于"SUBROUTIN E MAP(X，Y)"语句行的叙述中，不正确的是 A) 这是子程序的第一个语句 B) 字符串"MAP"是子程序名 C) 变量X是子程序的形参 D) 子程序执行后，MAP将返回整型数据（子程序无返回值，自定义函数才有） A (11)FORTRAN表达式"2/4+0.5"的值是 A) 0.5 B) 1 C) 1.0 D) 0 D (12)FORTRAN表达式"MOD(4.9，1.3)"的值是 A)1.3 B)0.0 C)4.9 D)1.0 A (13下列FORTRAN运算符中，优先级最低的是 A）逻辑运算符.AND. B）算术运算符* C)关系运算符 >= D)算术运算符＋ A (14下列语句函数的定义中正确的是 A）F(X，Y)＝(X＋Y)/(X*Y)＋7.0 B）FUNCTION FUN(I，J，K)＝3*I＋2*J＋0.5*K
fortran程序案例题汇编(14道)
1．Fibonacci数列定义如下： F1=1 F2=1 F n=F n-1+F n-2 (n>2) 求Fibonacci数列的前30项。 integer f(30),i f(1)=1 f(2)=2 do i=3,30 f(i)=f(i-1)+f(i-2) enddo print*,f end 2．输入10个学生的总分，求每个学生的名次integer s(10),a(10),i,j do i=1,10 read*,s(i) enddo do i=1,10 a(i)=1 do j=1,10 if(s(i) 3．给定一组数，按照从小到大的顺序输出。integer a(10) integer p do i=1,10 read *,a(i) enddo do j=1,9 p=j do i=j+1,10 if (a(i) 4．输入若干名学生的学号和三门课程（语数英）的成绩，要求从键盘输入一个学生的学号，能打印出该学生的三门功课成绩和总分。 character*6,dimension(:),allocatable::xue integer,dimension(:,:),allocatable::g integer,dimension(:),allocatable::zong integer i,j,n character*6,xh print *,"请输入学生的个数" read *,n allocate(xue(n)) allocate(g(n,3)) allocate(zong(n)) do i=1,n read *,xue(i),(g(i,j),j=1,3) enddo do i=1,n zong(i)=0 do j=1,3 zong(i)=zong(i)+g(i,j) enddo enddo print *,"请输入你要打印的学生的学号" read *,xh do i=1,n if(xue(i)==xh)then print *,(g(i,j),j=1,3),zong(i) exit endif enddo end FORTRAN笔记 2014.10.07 目录 第七讲_FORTRAN的基本知识.ppt (2) FORTRAN语言程序设计初步 (2) FORTRAN源程序的书写格式（以77为例） (2) 变量 (2) 变量类型 (2) 算术运算符和运算优先级 (3) 赋值语句 (3) 参数语句（PARAMETER语句） (3) END语句 (3) PAUSE语句 (3) 逻辑运算和选择结构 (4) 关系表达式 (4) FORTRAN中数组的定义及使用 (4) 其他 (5) 1. fortran语言定义CHARACTER*6 TTL(14,3),CNAM(400)是什么意思? (5) 2. fortran里character*10 是什么意思 (5) 3. Fortran中kind是什么函数？ (5) 第七讲_FORTRAN的基本知识.ppt FORTRAN语言程序设计初步 FORTRAN是Formula Translation的缩写，意为“公式翻译”，它是为科学、工程问题或企事业管理中的那些能够用 数学公式表达的问题而设计的，其数值计算的功能较强。 常用的是FORTRAN77和FORTRAN90两种标准。 1、一个程序由若干个程序单位组成。主程序和每一个子程序分别是一个独立的程序单位。 2、每一个程序单位都是以“END”结束的。 3、一个程序单位包括若干行。 1）语句行。由一个FORTRAN语句组成。 2）非语句行，即注释行。 4、FORTRAN程序中的语句可以没有标号，也可以有标号，根据需要而定。标号的作用是标志一个语句以便被其 他语句引用。 5、一个程序单位中各类语句的位置是有一定规定的。 6、FORTRAN源程序必须按一定的格式书写。 FORTRAN源程序的书写格式（以77为例） 每一行有80列，分别如下： 1、第1-5列为标号区。一行中第一列为“C”或“*”，该行即被认为是注释行。 2、第6列为“续行标志区”，如果在一行的第6列上写一个非空格和非零的字符，则该行作为其上一行的续行。 3、第7-72列为语句区。 4、第73-80列，注释区。 变量 变量名：一个变量需要用一个名字（变量名）来识别。在同一个程序单位中不能用同一个变量名代表两个不同的变 量。 FORTRAN的变量名按以下规则选定： 1）第一个字符必须是字母，即变量名必须以字母开头； 2）在一个字母后面可以跟1-5为数字或字母。 如果选定的变量名超过6个字符，则只有前面6个字符有效。 注：在变量名中大写与小写字母是等价的。 变量类型 整型变量Integer、实型变量Real、双精度变量Double Precision、复型变量Complex、逻辑型变量Logical和字符型变量Character。 1、隐含约定（I-N规则） FORTRAN规定：在程序中的变量名，凡以字母I，J，K，L，M，N六个字母开头的，即认为该变量为整型变量。 在程序中，凡是变量名以字母I,J,K,L,M,N,i,j,k,l,m,n开头的变量被默认为整型变量，以其他字母开头的变量被 默认为实型变量。 2、用类型说明语句确定变量类型 1）INTEGER语句（整型说明语句） 2）REAL语句（实型说明语句） 3）DOUBLE PRECISION语句（双精度说明语句） 4）COMPLEX语句（复型说明语句） 5）LOGICAL语句（逻辑型说明语句） 1、找出100-150之间和400-450之间能被9整除的数 2、找出100-999之间的回文数 3、找出水仙花数 4、输出“*”图形 5、输入n个数，找出大于平均值的数和最小数 6、把一个数值型数组的相同数删除到只剩一个 7、形成一个5*5矩阵，对角线元素为“i”，其余为“j”program juzhen implicit none integer I,J integer,parameter::size=5 integer::a(size,size) forall(I=1:size,J=1:size,I>J) a(I,J)=j forall(I=1:size,J=1:size,I==J) a(I,J)=i forall(I=1:size,J=1:size,I integer,allocatable::c(:) write(*,*) "输入数组A的数据个数" read(*,*) n write(*,*) "输入数组B的数据个数" read(*,*) m l=m+n allocate (a(n)) allocate (b(m)) allocate (c(l)) write(*,*) "从小到大输入A的元素" do i=1,n read(*,*) a(i) end do write(*,*) "从小到大输入B的元素" do i=1,m read(*,*) b(i) end do do i=1,n c(i)=a(i) end do i=1 do while(i 1)实例3—求多个半径下的圆周长 ! z3.f90 --Fortran95 ! FUNCTIONS: ! z3 - Entry point of console application. !************************************************************************* ! PROGRAM: z3 ! PURPOSE: Entry point for the console application. !************************************************************************ program z3 ! 求多个半径下的圆周长 ! 主程序 ! PROGRAM Z3 PRINT *, 'R=',1.2,'C=',C(1.2) PRINT *, 'R=',3.4,'C=',C(3.4) PRINT *, 'R=',15.6,'C=',C(15.6) PRINT *, 'R=',567.3,'C=',C(567.3) END program z3 !子程序 FUNCTION C(R) PI=3.1415926 C=2*PI*R RETURN ! Body of z3 end 2)实例4—键盘与显示器输入/输出 a)Fortran 基本操作 b)程序指令 ! ZXZ_I_O.f90 ! FUNCTIONS: ! ZXZ_I_O - Entry point of console application. ! PROGRAM: ZXZ_I_O ! PURPOSE: Entry point for the console application. !***************输入、输出样式种种************************** program ZXZ_I_O implicit none !变量声明的位置 INTEGER(2) i; INTEGER(4) j; INTEGER(4) m; REAL n INTEGER A,B ! Variables PRINT*,'输入整数A'; READ*, A PRINT*,'输入整数B'; READ*, B B=A+B PRINT*,'B=A+B=',B WRITE(*,*) 'A*B=',A*B PRINT* ,'以上为计算机的计算结果，注意B的值' FORTRAN程序设计复习题 一、选择题 B (1)下列各FORTRAN表达式中合法的是 A) S+T*2P >= B) .NOT. (A*B+C) C) A2+B2/(C+D) <= D) (A+B).NOT．A*B．GT．(.NOT.只跟一个表达式) C (2)数学式(3/5)ex+y的FORTRAN表达式是 A) 3*EXP(X+Y)/5 B) 3*E* *(X+Y)/ C) (3/5)*EXP(X+Y)D) EXP(X+Y) D (3)下列FORTRAN77表达式中不合法的是 A) A．GT．B．EQV．C．GT．D B) A．AND．B．AND．C．AND．D C) .NOT．(X．LE．D) A．LT．B．LT．C．LT．D D(4)下列叙述中不正确的是 A) FORTRAN子程序可以单独编译 B) 对一个FORTRAN源程序进行编译和连接无误后可生成可执行文件 C) 即使编译和连接都正确无误，FORTRAN程序运行时仍可能出错 D) FORTRAN连接的主要任务是把函数库中的函数翻译成机器指令(正确描述：主要任务为连接目标文件) B (5)在下列FORTRAN77运算符中，优先级最高的是 A) .AND. B) .NOT. C) .OR. D) .EQ. B (6)FORTRAN表达式"6/5+9/2**3/2"的值为 A) 33 B) 1 C) 5 D) 3 A (7)下列FORTRAN77表达式中，合法的是： A) .AND.. B) 10.0 C) D) 提示：A）相当于 .AND.（.NOT.（）） D (8)关于编译一个FORTRAN源程序文件，下列说法中错误的是 A) 允许编译只有一个主程序而没有子程序的源文件 B) 允许编译有多个子程序的源文件 C) 允许编译只有一个子程序而没有主程序的源文件 D) 允许编译有多个主程序的源文件 C (9)在FORTRAN77源程序中，续行标志符必须放在 A) 第1列 B) 第1-6列C) 第6列D) 第5列 D (10)下列关于"SUBROUTIN E MAP(X，Y)"语句行的叙述中，不正确的是 A) 这是子程序的第一个语句 B) 字符串"MAP"是子程序名 C) 变量X是子程序的形参D) 子程序执行后，MAP将返回整型数据 提示：子程序无返回值，自定义函数才有） A (11)FORTRAN表达式"2/4+"的值是 A) B) 1 C) D) 0 提示：2/4默认等于整型，=》 D (12)FORTRAN表达式"MOD，"的值是 A) B)0.0 C) D) A (13下列FORTRAN运算符中，优先级最低的是 A）逻辑运算符.AND. B）算术运算符* F O R T R A N 第2章FORTRAN90基础知识： 程序单元的概念: fortran90程序是一种分块形式的程序，整个程序由若干程序单元块组成。每个程序只有一个主程序单元。各单元体的程序体形式上相同。 程序单元可以是主程序、子程序（外部过程或内部过程）、模块MODULE （供其他程序单元引用即把该程序单元内的全部语句复制到程序单元中）或块数据程序单元BLOCK 。 语言元素：指在程序设计中用到的基本成分，如字符集、常量、变量、记号（标号、关键字、名字、常数、运算符和定界符）以及其他的基本结构等。 字符集：英文字母、阿拉伯数字、下划线、21个特殊字符（$和？没有规定用法）。 数据结构： 整型INTEGER （34-2下划线后面是种别参数），n 位字长的计算机的数据表示范围一般为12~211-+---n n ，种别参数由内在询问函数KIND 送回，种别值提供的最小范围由内在函数SELECTED-INT-KIND 送回，表示十进制幂的范围由内在函数RANGE 送回； 实型REAL ，小数形式和指数形式；复型COMPLEX （种别类 型参数值取实部与虚部中较大者）； 字符型CHARACTER，由一对单撇号或双撇号之间的字符序列组成； 逻辑型LOGICAL。 派生数据类型TYPE； 数组INTEGER,DIMENSION(1,50)::A，可直接对数组元素进行运算如3*A,SQRT(A)； 字符子串，在字符串CHARACTER(LEN=80)::ROW中，ROW(1:3)就表示字符串ROW中第1到第3个元素组成的子串。 变量名的命名规则：不能超过31个字符；组成成分是字母、数字、下划线；第一个字符必须是字母。 第3章基本语句： 类型说明语句：类型说明（种别说明），属性说明：：变量名表 尽量避免把两个实数作相等或不相等的比较。淘汰隐式说明IMPLICIT NONE 种别说明：种别参数即对可移植数据精度和范围进行选择的机制 KIND(X) 返回变元X的种别参数值 SELECTED-REAL-KIND(n,m) 产生一个种别值，它表示某一精度和范围。N指十进制有效位数，m指明值范围内以10为底的幂次。 计算圆周率 REAL R,R1,R2,PI ISEED=RTC() N0=0 N=300000 DO I=1,N R1=RAN(ISEED) R2=RAN(ISEED) R=SQRT(R1*R1+R2*R2) IF(R<1.0)N0=N0+1 END DO PI=4.0*N0/N WRITE(*,*)PI END 一）蒙特卡洛计算生日问题 假设有N个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直觉相反，只需23个人便有大于50％的几率人群中至少有2个人生日相同。 INTEGER M(1:10000), NUMBER1(0:364), NUMBER2 REAL X,Y ISEED=RTC() DO J=1, 10000 NUMBER1=0 X=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 DO END DO DO I=1,10000 IF(M(I).LE.23) SUM=SUM+1 END 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，实验次数IMAX" READ(*,*) JMAX,IMAX ISEED=RTC() DO J=1,JMAX !第几天实验 X=0 !!! DO I=1,IMAX !第几步跳跃 RN=RAN(ISEED) IF(RN<0.5)THEN X=X+1 ELSE X=X-1 END IF XX(J,I)=X*X END DO END DO OPEN(1,FILE="C:\DIF1.DAT") DO I=1,IMAX XXM=0.0 XXM(I)=1.0*SUM(XX(1:JMAX,I))/JMAX !! WRITE(1,*) I, XXM(I) END DO CLOSE(1) END 三维的！ 三）通过该程序了解FORTRAN语言如何画图（通过像素画图）USE MSFLIB INTEGER XR，YR !在的区域中画一个圆 PARAMETER XR=400,YR=400 INTEGER R, S(1:XR,1:YR) X0=XR/2 ! 圆心位置X0，YO Y0=YR/2 R=MIN(X0-10,Y0-10) !圆半径 S=0 !像素的初始状态（颜色） DO I=1,XR DO J=1,YR IF((I-X0)**2+(J-Y0)**2<=R**2)S(I,J)=10 IER=SETCOLOR(S(I,J)) IER=SETPIXEL(I,J) END DO END DO END 程序1 PROGRAM chaper1272 INTEGER :: A(10)=(/5,7,4,8,12,2,10,3,9,11/) INTEGER :: sum=0,ave !打开一数据文件,设置一个有格式直接存取文件,将10个数分2个记录写入文件。 !数据文件生成2个记录,每个记录行长度相同,记录长度为25 OPEN(1,FILE='input21.dat',FORM='FORMATTED',ACCESS='D IRECT',RECL=25) WRITE(1,"(5I5)",REC=1)(A(I)+10,I=1,5) !按格式说明将头5个数写入第1个记录 WRITE(1,"(5I5)",REC=2)(A(I)+10,I=6,10)!按格式说明将后5个数写入第2个记录 READ(1,"(5I5)",REC=2)(A(I),I=6,10) !按格式说明从第2个记录中读取后5个数 READ(1,"(5I5)",REC=1)(A(I),I=1,5) !按格式说明从第1个记录中读取头5个数 DO I=1,10 sum=sum+A(I) ENDDO ave=sum/10 !打开一个最大记录长度为22的有格式顺序存取文件 OPEN(2,FILE='input22.dat',FORM='FORMATTED',ACCESS='D IRECT',RECL=22) WRITE(2,"('10个数之和为：',I5)",REC=1) sum !输出1记录行,记录长度为22 WRITE(2,"('10个数平均值为：',I5)",REC=2)ave !输出1记录行,记录长度为22 WRITE(2,"(A)",REC=3) '程序运行正常结束。' !输出一个记录行,记录长度为22 END 程序2 PROGRAM average REAL sum,ave INTEGER n OPEN(1,file='score.dat') PRINT*,'正在统计平均成绩,请等待。' sum=0.0;n=0 DO READ(1,*,END=100) S sum=sum+s n=n+1 ENDDO 100 ave=sum/n 第一次课绘制简单分子 一、画一个苯酰胺 目的：介绍Materials Visualizer中画结构的工具。 用到的模块：Materials Visualizer 化学家每天都要处理很多种类的小分子和中间物。所以容易的创建模型对建模环境都是很重要的。苯酰胺是典型的小分子结构。以下通过建立他的结构来学习Materials Studio。下面是要建立的苯酰胺的结构： Benzamide 1.创建3D文档： 从菜单中选择File | New...打开New Document对话框。选择3D Atomistic Document（三维原子文档），按OK。建立了一个三维窗口，工程管理器中显示建立了名为3D Atomistic Document.xsd的文件。在工程管理器这个文件名上右击鼠标，选择Rename改名。键入my_benzamide的新名字，按回车。选择File | Save 命令，或按标准工具条中的按钮。在my quickstart文件夹（每个工程都对应一个同名的文件夹）中建立了名为my_benzamide.xsd的文件。 2. 改变到Ball and Stick球棍模型显示方式。 三维窗口中右击鼠标，选择Display Style，打开Display Style对话框，在Atom 选项卡上设置。Materials Studio能在任何显示方式下添加原子。 3. 画环和原子链。 在草画工具条上单击Sketch Ring 按钮，鼠标移到三维窗口。鼠标变为 铅笔行状提示你处于草画模式。鼠标榜的数字表示将要画的环包括的原子数目。可以通过按3-8的数字键改变。确保这个数字为6，三维窗口中单击。画出了一个6个C原子的环。如果安装ALT键单击，产生共振键。 现在单击草画工具条Sketch Atom 按钮，这是通用添加原子工具，可加 入任何元素，默认加入C原子。如下在环上加入两个C原子。在环上移动鼠标，当一个原子变为绿色时单击，键的一端就在这个原子上，移动鼠标再单击就加入了一个C原子，再移动，并双击。这样在环上加入了两个原子。另一种结束添加原子的方法是在最后一个原子位置单击，然后按ESC键。注意，新加入的原子的化学键已经自动加上。 注意：你可以按Undo 按钮取消错误操作。 4、加入氧原子。 按Sketch Atom按钮旁的向下按钮，显示可选元素，选择氧Oxygen，在支链上移动鼠标，当变为蓝色显示时单击，这个原子就有了一个化学键，移动鼠标并双击。加入了O原子。在3D窗口工具条上按按钮，进入了选择模式。 5. 编辑元素类型。 单击链末端的C原子，选定它。选定的对象用黄色显示。按Modify Element按钮旁的箭头，显示元素列表，选择Nitrogen氮，选定的原子就变为了氮原子。单击三维窗口中空白地方，取消选择，就可以看到这种变化了。 6．编辑键类型。 在三维窗口中在C和O原子中间单击选定C－O键。选定的键以黄色显示。按下SHIFT键，单击其它三个相间的键。现在选定了三个C－C键和一个C－O键。 单击Modify Bond 按钮旁的向下按钮，显示键类型的下拉列表，选择Double Bond双键。取消选定。 7. 调整氢原子和结构 现在可以给结构自动加氢。单击Adjust Hydrogen 按钮，自动给模型加入 计算材料学（Computational Materials Science），是材料科学与计算机科学的交叉学科，是一门正在快速发展的新兴学科，是关于材料组成、结构、性能、服役性能的计算机模拟与设计的学科，是材料科学研究里的“计算机实验”。它涉及材料、物理、计算机、数学、化学等多门学科。 计算材料学- 学科介绍 计算材料学（Computational Materials Science），是材料科学与计算机科学的交叉学科，是一门正在快速发展的新兴学科，是关于材料组成、结构、性能、服役性能的计算机模拟与设计的学科，是材料科学研究里的“计算机实验”。它涉及材料、物理、计算机、数学、化学等多门学科。计算材料学主要包括两个方面的内容：一方面是计算模拟，即从实验数据出发，通过建立数学模型及数值计算，模拟实际过程；另一方面是材料的计算机设计，即直接通过理论模型和计算，预测或设计材料结构与性能。前者使材料研究不是停留在实验结果和定性的讨论上，而是使特定材料体系的实验结果上升为一般的、定量的理论，后者则使材料的研究与开发更具方向性、前瞻性，有助于原始性创新，可以大大提高研究效率。因此，计算材料学是连接材料学理论与实验的桥梁。 计算材料学- 研究领域 材料的组成、结构、性能、服役性能是材料研究的四大要素，传统的材料研究以实验室研究为主，是一门实验科学。但是，随着对材料性能的要求不断的提高，材料学研究对象的空间尺度在不断变小，只对微米级的显微结构进行研究不能揭示材料性能的本质，纳米结构、原子像已成为材料研究的内容，对功能材料甚至要研究到电子层次。因此，材料研究越来越依赖于高端的测试技术，研究难度和成本也越来越高。另外，服役性能在材料研究中越来越受到重视，服役性能的研究就是要研究材料与服役环境的相互作用及其对材料性能的影响。随着材料应用环境的日益复杂化，材料服役性能的实验室研究也变得越来越困难。总之，仅仅依靠实验室的实验来进行材料研究已难以满足现代 新材料研究和发展的要求。然而计算机模拟技术可以根据有关的基本理论，在计算机虚拟环境下从纳观、微观、介观、宏观尺度对材料进行多层次研究，也可以模拟超高温、超高压等极端环境下的材料服役性能，模拟材料在服役条件下的性能演变规律、失效机理，进而实现材料服役性能的改善和材料设计。因此，在现代材料学领域中，计算机“实验”已成为与实验室的实验具有同样重要地位的研究 手段，而且随着计算材料学的不断发展，它的作用会越来越大。 计算材料学的发展是与计算机科学与技术的迅猛发 展密切相关的。从前，即便使用大型计算机也极为困难的一些材料计算，如材料的量子力学计算等，现在使用微机就能够完成，由此可以预见，将来计算材料学必将有更加迅速的发展。另外，随着计算材料学的不断进步与成熟，材料的计算机模拟与设计已不仅仅是材料物理以及材料 计算理论学家的热门研究课题，更将成为一般材料研究人员的一个重要研究工具。由于模型与算法的成熟，通用软件的出现，使得材料计算的广泛应用成为现实。因此，计算材料学基础知识的掌握已成为现代材料工作者必备的 技能之一。 计算材料学涉及材料的各个方面，如不同层次的结构、各种性能等等，因此，有很多相应的计算方法。在进行材料计算时，首先要根据所要计算的对象、条件、要求等因素选择适当的方法。要想做好选择，必须了解材料计算方法的分类。目前，主要有两种分类方法：一是按理论模型和方法分类，二是按材料计算的特征空间尺寸(Characterist ic space scale)分类。材料的性能在很大程度上取决于材料的微结构，材料的用途不同，决定其性能的微结构尺度会有很大的差别。例如，对结构材料来说，影响其力学性能的结构尺度在微米以上，而对于电、光、磁等功能材料来说可能要小到纳米，甚至是电子结构。因此，计算材料学的研究对象的特征空间尺度从埃到米。时间是计算材料学的另一个重要的参量。对于不同的研究对象或计算方法，材料计算的时间尺度可从10-15秒（如分子动力学方法等）到年（如对于腐蚀、蠕变、疲劳等的模拟）。对于具有不同特征空间、时间尺度的研究对象，均有相应的材料计算方法。 目前常用的计算方法包括第一原理从头计算法，分子动力学方法，蒙特卡洛方法，有限元分析等。 计算材料学是目前材料科学中发展最快的科目。随着大量的论文发表和科研材料的科学家数量的快速增加，现在已有条件检验计算材料学对材料科学的影响。可以回顾以下计算机模拟在新材料的发展上起过什么作用，或许更重要的是，总结在材料性能的基础研究方面有哪些突破应该归公与计算机模拟。 有意思的是，人们对模拟方法的期望竟产会超过实际的结果，然而这些过高的期望往往更多来自非专业认识而不是那些正在做这些工作的饿人。如果问一个从事计算材料学的人关于现在使用的方法，他回非茶馆内谨慎地回答：尽管我们的方法很可靠，但仍需要大的发展。这些方法有欠缺并不奇怪，计算材料学只有几十年的历史。因此，我们在关注它现今的地位时，必须同样关注仿镇与建模的可靠性。这门年轻的学科，已经有诸多长足的发展，涉及到许多包含多种距离尺度的现象。 也许最原始的计算材料学是计算固体的电子结构。这些计算显然已经非常成功地表述了材料的结构和性质。现在，对于许多晶体材料，预计的点阵常数和实验值仅相差百分之几。最近的弹性常数计算方法得到了与实验值非常吻合的结果，而且实行起来也比实验容易得多。多体理论的发展，使得目前已能对简单半导体禁带宽度进行预测。 但是，基于这样一些成果，电子结构计算往往表述得似乎比实验值更精确。实际上，着些计算含有很多近似，而且很容易发生误导而得到错误的结论。此外，近似法限 Monte Carlo方法模拟及其应用 沈合平 （上海大学材料学院学院，上海200072） 摘要：扫描电子显微学中使用二次电子和背散射电子作为成像信号时遇到了一些问题，通过计算机模拟可以很好地解决这些问题.本文首先介绍了Monte Carlo方法，再介绍了Monte Carlo方法在扫描电镜模拟中的应用，并且总结了计算机模拟的发展方向。 关键词：扫描电子显微镜；衬度；Monte Carlo；计算机；模拟 Abstract:Scanning electron microscopy using the secondary electron and backscattered electron imaging signal as many problems encountered,which can be solved by computer simulation.This paper describes the Monte Carlo method, and then introduced the Monte Carlo method in the simulation of a scanning electron microscope, and summed up the direction of the computer simulation. Keywords:SEM;Contrast;Monte Carlo;computer;simulation 1. 引言 扫描电子显微学中使用二次电子和背散射电子作为成像信号。对于研究材料的表面形貌非常重要。低能二次电子主要反映试样的表面形貌特征,而较高能量的背散射电子既可在一定程度上反映试样的表面特征,也可表征试样的内部成分和结构差异。 多数二次电子的能量很小,从表面发射时的峰值能量仅为数eV,故其在材料内部的运动范围有限,只有那些在表面附近产生的二次电子才能从试样表面发射出来。二次电子主要用于表征试样的表面形貌特征。而具有较高能量的背散射电子则是入射电子在深入试样的内部后由于多次散射效应再从表面发射出来的那些电子,它们既包含试样的表面信息,也含有试样结构差异和内部成分的信息。 当用场发射扫描电镜观察数十纳米尺度以下的小颗粒时,衬度与大尺度颗粒的情形相差很大,二次电子图像仅仅呈现出一些亮点和较弱的光点,而背散射电子图像则显示大量的亮点,因此难以判定颗粒在基底表面的分布情况。由于该颗粒/基底体系的扫描电镜图像衬度的形成机制较为复杂[1,2],尚未得到理解,因而限制了从SEM图像中提取出有用信息。 如果能模拟计算二次电子和背散射电子信号产生的过程,将有助于理解扫描电子显微镜的成像和图像衬度机理。因此,研究者们利用电子散射轨迹模拟的Monte Carlo方法已做了一些研究,如Gauvin[3]模拟了一个嵌入到均匀基底内的球的背散射电子像和X射线像;Radzimski和Russ[4]基于利用Rutherford散射截面和Bethe阻止本领的单散射模型,模拟了多层多元素试样在二维方向上表面形貌的背散射电子像;而Yan和Gomati[5]则开发了一个三维的Monte Carlo程序用来模拟一些比较复杂试样的背散射电子和Auger电子像,但这个程序也要求试样的几何结构必须能被解析地表达出来,所以能模拟的情形仍然有限。特别是这些研究中还不能得到二次电子像,主要原因是因为计算二次电子发射的产额相当困难。 因此,无论从计算方法还是从应用前景,模拟成分非均匀和形貌特殊的试样成像(特别是二次电子像)均成为有重要意义的研究工作。 2.Monte Carlo方法 Monte Carlo 方法是在二战期间产生和发展起来的。他的奠基者是美籍匈牙利人数学家冯诺伊曼（J.Von Neumann 1903-1957）。作为Monte Carlo 方法的最初应用，是解决蒲丰氏 1. Fibonacci数列定义如下： F1=1 F2=1 F n=F n-1+F n-2 (n>2) 求Fibonacci数列的前30项。integer f(30),i f(1)=1 f(2)=2 do i=3,30 f(i)=f(i-1)+f(i-2) enddo print*,f end 2.输入10个学生的总分，求每个学生的名次integer s(10),a(10),i,j do i=1,10 read*,s(i) enddo do i=1,10 a(i)=1 do j=1,10 if(s(i) 3．给定一组数，按照从小到大的顺序输出。integer a(10) integer p do i=1,10 read *,a(i) enddo do j=1,9 p=j do i=j+1,10 if (a(i) 4．输入若干名学生的学号和三门课程 (语数英) 的成绩，要求从键盘输入一个学生的学号，能打印出该学生的三门功课成绩和总分。 character*6,dimension(:),allocatable::xue integer,dimension(:,:),allocatable::g integer,dimension(:),allocatable::zong integer i,j,n character*6,xh print *," 请输入学生的个数" read *,n allocate(xue(n)) allocate(g(n,3)) allocate(zong(n)) do i=1,n read *,xue(i),(g(i,j),j=1,3) enddo do i=1,n zong(i)=0 do j=1,3 zong(i)=zong(i)+g(i,j) enddo enddo print *," 请输入你要打印的学生的学号" read *,xh do i=1,n if(xue(i)==xh)then print *,(g(i,j),j=1,3),zong(i) exit endif enddo end 对于FORTRAN的初学者。这些例子可作为小练习。 1.例题：计算工资问题。每小时工资为RATE,如果工作超过40小时，加班呢部分工资是正常时间工资的1.5倍。 C Payroll with overtime program payroll real rate, hours, pay read (*,*) rate, hours if (hours>40) then pay=40*rate+(hours-40)*1.5*rate else pay=hours*rate END IF print *,"rate=", rate print *, "hours=", hours print *,"pay=",pay end 2.学生成绩问题。大于80为A级。大于60小于80为B级。小于60为C级。IF的嵌套。注意空格可以看清楚else if ,end if,print的内容. PROGRAM GRADESTUDENT READ *,GRADE IF (GRADE .GE. 80) THEN PRINT *,GRADE,"GRADE=>A" ELSE IF (GRADE.LT.60) THEN PRINT*,GRADE,"GRADE=>C" ELSE PRINT*,GRADE,"GRADE=>B" END IF END IF END 3.三个数按从小到大排序。 PROGRAM MAXMIN REAL A,B,C,T READ *,A,B,C IF (A.GT.B) THEN T=A A=B B=T ELSE END IF IF (B.GT.C) THEN T=B B=C PROGRAM EQUATION READ (*,*) X IF (X .GE. 0.0) Y=3*X+6 IF (X .LT. 0.0) Y=-X**2+2*X-8 PRINT *,"X=",X, "Y=",Y END 6.CONTINUE语句。为了使循环的起止范围清晰，使用终端语句 例：求5！ program ex1 FACT=1 do 10, N=1,5,1 ##DO 语句标号，循环不变量=初值，终值，步长FACT=FACT*N 10continue print*,FACT end 求1！+2！+3！+···+N！ READ *,N FACT=1(完整)Fortran经典编程语言笔记(你值得拥有)
中国海洋大学Fortran程序考试典型题
fortran程序实例
FORTRAN程序设计复习题及答案
FORTRAN知识点总结
几个fortran程序
fortran程序30个
计算材料学讲稿-6
计算材料学
计算材料学
fortran程序案例题汇编(14道)
FORTRAN经典入门程序20例