fortran程序匯總

1、計(jì)算圓周率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(r1.0)n0=n0+1 end do pi=4.0*n0/n write(*,*)pi end 一）蒙特卡洛計(jì)算生日問(wèn)題假設(shè)有 n 個(gè)人在一起，各自的生日為365 天之一，根據(jù)概率理論，與很多人的直覺(jué)相反，只需 23 個(gè)人便有大于 50的幾率人群中至少有2 個(gè)人生日相同。integer m(1:10000), number1(0:364), number2 real x,y ise

2、ed=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

3、 end 二） monte carlo simulation of one dimensional diffusion 蒙特卡羅計(jì)算一維擴(kuò)散問(wèn)題integer x,xx(1:1000,1:1000) real xxm(1:1000) ! x:instantaneous position of atom ! xx(j,i) ：x*x ,j: 第幾天實(shí)驗(yàn) ,i:第幾步跳躍! xxm(i): the mean of xx write(*,*) 實(shí)驗(yàn)天數(shù) jmax ，實(shí)驗(yàn)次數(shù)imax read(*,*) jmax,imax iseed=rtc() do j=1,jmax ! 第幾天實(shí)驗(yàn)x=0 ! do

4、 i=1,imax ! 第幾步跳躍rn=ran(iseed) if(rn0.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 三維的！三）通過(guò)該程序了解fortran 語(yǔ)言如何畫圖（通過(guò)像素畫圖）use msflib integer xr，yr !在的區(qū)域中畫一個(gè)圓 parameter xr=4

5、00,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 !像素的初始狀態(tài)（顏色） 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 四）畫一個(gè)圓（ 1、如何選出晶界區(qū)域；2、進(jìn)一步加深對(duì)畫圖的理解）use msflib integer xr，yr !在的區(qū)域中畫一個(gè)圓parameter xr=400,

6、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 !像素的初始狀態(tài)（顏色）do i=1,xr do j=1,yr if(i-x0)*2+(j-y0)*20)then ier=setcolor(9) else ier=setcolor(8) end if ier=setpixel(i,j) end do end do end 五）mc 模擬一個(gè)晶粒的縮小u

7、se msflib parameter ir=400,jr=400 integer is(0:ir+1,0:jr+1),tmax,isn(1:8),nstate,t,nr,ix,iy write(*,*)please input the time step read(*,*)tmax iseed=rtc() ! 定義圓心和半徑irc=ir/2 jrc=ir/2 r=min(irc,jrc)-10 ! 定義基體和圓晶粒分別為狀態(tài)1、狀態(tài) 2 is=1 do i=1,ir do j=1,jr distance=sqrt(1.0*(i-irc)*2+1.0*(j-jrc)*2) if(distanc

8、e.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) ! 尋找晶粒邊界，計(jì)算能量，改變狀態(tài)。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=

9、count(isn.ne.is(ix,jy) ! 如果不是圓晶粒邊界，則跳出，重新循環(huán)if(e0.eq.0)cycle ! 隨機(jī)尋找一個(gè)相鄰點(diǎn)nr=8*ran(iseed)+1 nstate=isn(nr) ! 判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)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,c

10、ount(is.eq.2) enddo close(1) end 六）蒙特卡羅模擬基體上單晶粒形核長(zhǎng)大use msflib parameter ir=400,jr=400 integer is(0:ir+1,0:jr+1),tmax,isn(1:8),nstate,t,nr,ix,iy write(*,*)please input the time step read(*,*)tmax iseed=rtc() ! 定義圓心和半徑irc=ir/2 jrc=ir/2 ! 定義基體和圓晶粒分別為狀態(tài)10、狀態(tài) 2 is=10 is(irc,jrc)=2 open(1,file=e:luke.dat)

11、 ! 尋找晶粒邊界，計(jì)算能量，改變狀態(tài)。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) ! 如果不是圓晶粒邊界，則跳出，重新循環(huán)if(e0.eq.0)cycle ! 隨機(jī)尋找一個(gè)相鄰點(diǎn)nr=8*ran(iseed)+1 nstat

12、e=isn(nr) ! 判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)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 七）蒙特卡洛模擬基體上多晶粒形核長(zhǎng)大，use msflib ！定義常數(shù)名parameter ir=400

13、,jr=400,nmax=100 ！定義變量：體積能變量（一維數(shù)組） ；基體各像素點(diǎn)變量（二維數(shù)組） ，鄰居的取向號(hào)（一維數(shù)組）等；integer is(0:ir+1,0:jr+1),tmax,isn(1:8),nstate,t,nr,ix0,iy0,ix,iy integer igv(0:nmax) ！讀取晶粒生長(zhǎng)步長(zhǎng)；write(*,*)please input the time step read(*,*)tmax iseed=rtc() ! 定義基體體積能為10,所有晶粒體積能為1 igv=1 igv(0)=10 ! 定義晶粒長(zhǎng)大位置（ 100 個(gè)形核點(diǎn)初始形核位置） ，并附已不同的取

14、向號(hào)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*

15、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) ! 如果不是晶粒邊界，則跳出，重新循環(huán)if(e0.eq.0)cycle ! 隨機(jī)尋找一個(gè)相鄰點(diǎn)nr=8*ran(iseed)+1 nstate=isn(nr) ! 判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)rd=ran(iseed) e=count(isn.ne.nstate) ig=

16、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 八）元胞自動(dòng)機(jī)模擬生命游戲程序（生命永不停止）use msflib parameter ir=400,jr=400,nmax=40000 !nmax-隨機(jī)產(chǎn)生的生命種子integer is(0:1001,0:10

17、01),is1(0:1001,0:1001),isn(1:8), tmax, num ！is-基體的二維數(shù)組print*,please input loop(tmax) read*, tmax iseed=rtc() is=15 ！ “死”的狀態(tài)，基體為白色！賦予生命的種子，“活”的狀態(tài) 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+

18、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 ！判斷鄰居狀態(tài)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

19、=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 do is=is1 end do end 九）元胞自動(dòng)機(jī)模擬單晶長(zhǎng)大use msflib parameter ir=400,jr=400 integer is(0:ir+1,0:jr+1),tmax,isn(1:8),nstate,t,nr,ix0,iy0,ix,jy ! 根據(jù)過(guò)去狀態(tài) is，定義現(xiàn)在狀態(tài)為i

20、s1；為了突出邊界，特別定義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 !=ir*ran(iseed)+1 jrc=jr/2 ! 定義基體體積能為10,晶粒體積能為 1 is=8 is(irc,jrc)=1 ! 在循環(huán)前定義現(xiàn)在狀態(tài)數(shù)組is1 的初始值is1=is open(1,file=e:luke.dat) do t=1,tmax ! 每次循環(huán)前重新定義過(guò)去狀態(tài)數(shù)組is is=is1 ! 邊界條件is(0,0

21、: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) ! 整個(gè)基體上，各個(gè)點(diǎn)上的狀態(tài)都要根據(jù)規(guī)則改變，而非隨機(jī)取點(diǎn)改變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

22、,jy) ! 如果不是晶粒邊界，則跳出，重新循環(huán)if(e0.eq.0)cycle ! 隨機(jī)尋找一個(gè)相鄰點(diǎn)nr=8*ran(iseed)+1 nstate=isn(nr) ! 判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)e=count(isn.ne.nstate) rd=ran(iseed) ig=nstate-is(ix,jy) de=e-e0+ig+2.5*rd-1.25 ! 用現(xiàn)在狀態(tài)數(shù)組 is1 記錄邊界狀態(tài)的改變if(de.lt.0.0)is1(ix,jy)=nstate enddo end do ! 每循環(huán) 20 次在顯示屏幕上刷新?tīng)顟B(tài)（顏色）do ix=1,ir do jy=1,jr

23、! 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,加以區(qū)分if(count(isn1.ne.is1(ix,jy).gt.0) isre=setcolor(15) isre=setpixel(ix,jy) ! end if enddo enddo ! ! 記錄循環(huán)

24、次數(shù)和對(duì)應(yīng)的晶粒面積write(1,*)t, sqrt(1.0*count(is.eq.1) enddo close(1) end 十）元胞自動(dòng)機(jī)模擬多晶長(zhǎng)大1.介紹蒙特卡羅和元胞自動(dòng)機(jī)的區(qū)別。元胞自動(dòng)機(jī)特點(diǎn)：時(shí)間、空間、狀態(tài)都離散的動(dòng)力學(xué)模型，利用確定性或概論性的轉(zhuǎn)換規(guī)則來(lái)描述在離散空間和時(shí)間上復(fù)雜系統(tǒng)演化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) ! 根據(jù)過(guò)去狀態(tài) is，定義現(xiàn)

25、在狀態(tài)為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 ! 定義晶粒長(zhǎng)大位置，并附能量，附帶生長(zhǎng)取向（研究形核數(shù)目與晶粒尺寸的關(guān)系）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 !

26、 在循環(huán)前定義現(xiàn)在狀態(tài)數(shù)組is1 的初始值is1=is ! open(1,file=e:luke.dat) do t=1,tmax ! 每次循環(huán)前重新定義過(guò)去狀態(tài)數(shù)組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) ! 尋找晶粒邊界，計(jì)算能量，改變狀態(tài)。! 整個(gè)基體上，各個(gè)點(diǎn)上的狀態(tài)都要根據(jù)規(guī)則改變，而非隨機(jī)取點(diǎn)改變do ix=1,ir do jy=1,jr ! isn=(/is(ix-

27、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) ! 如果不是晶粒邊界，則跳出，重新循環(huán)if(e0.eq.0)cycle ! 隨機(jī)尋找一個(gè)相鄰點(diǎn)nr=8*ran(iseed)+1 nstate=isn(nr) ! 判斷與相鄰點(diǎn)的能量差，并決定是否改變狀態(tài)e=count(isn.ne.nstate) rd=ran(iseed) ig=igv(nstate)-igv(is(ix,jy)

28、 de=e-e0+ig+2.5*rd-1.25 ! 用現(xiàn)在狀態(tài)數(shù)組 is1 記錄邊界狀態(tài)的改變if(de.lt.0.0)is1(ix,jy)=nstate ! ! isre=setcolor(mod(is1(ix,jy),15) ! isre=setpixel(ix,jy) enddo enddo ! 每循環(huán) 20 次在顯示屏幕上刷新?tīng)顟B(tài)（顏色）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

29、),is1(ix+1,jy-1),is1(ix+1,jy),is1(ix+1,jy+1)/) isre=setcolor(mod(is1(ix,jy),15) ! 如果是邊界，定義特別顏色15,加以區(qū)分if(count(isn1.ne.is1(ix,jy).gt.0) isre=setcolor(15) isre=setpixel(ix,jy) end if enddo enddo ! ! 記錄循環(huán)次數(shù)和對(duì)應(yīng)的晶粒面積write(1,*)t, sqrt(1.0*count(is1.gt.0) enddo close(1) end 十一） 相場(chǎng)方法模擬相變by using boundary tr

30、acking 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 msflib parameter imax=800,jmax=800,nmax=50000,nstep=10 ! 最多晶核數(shù) ,每步形成晶核數(shù)integer is(0:imax+1,0:jmax+1),ist(0:imax+1,0:jmax+1) integer in(1:8),jn(1:8) integer

31、 ni(1:nmax),nj(1:nmax),nt(1:nmax) !核的 i,j 坐標(biāo),形成時(shí)間integer nns,dcn,dcnm ! 運(yùn)行時(shí)間部，最近晶核狀態(tài) ,訪問(wèn) 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.dat) do it=1,tmax is=ist ! 邊界條件is(0,1:jmax)=is(imax,1:jmax

32、) 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)c

33、ycle 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 ! 當(dāng)前深刻形成的核 ,當(dāng)前時(shí)刻不長(zhǎng)大np=np+1 ! 只有 it 時(shí)刻以前形成的核長(zhǎng)大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 !訪問(wèn) 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.