[Fortran] 纯文本查看 复制代码
program sandpile
integer Z,h(3000,3000),L,t,ix,iy,an(3000,3000),a,n,na,aa
integer ax1,ay1,ah1(3000,3000),ax2,ay2,ah2(3000,3000)
integer ax3,ay3,ah3(3000,3000)
integer tn(100000),tx(100000),ty(100000),tx1(100000)
integer ty1(1000000),m,i,j,noapi,k,m1,n1,n2,n3,n4
integer a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14
integer a15,a16,a17,a18
integer na1,na2,na3,na4,na5,na6,na7,na8,na9,na10,na11
integer na12,na13,na14,na15,na16,na17,na18,na19,na20
integer na21,na22,na23
call CPU_time(time_begin)
Z=2
L=32
km=8
noap=1000000
Do ix=1,L !設ix平面從一到L
Do iy=1,L !設iy平面從一到L
h(ix,iy)=rand()*Z+1 !晶格內的高度1~4隨機
enddo
enddo
p1=0
p2=0
p3=0
p4=0
Do noapi=1,noap
if(km.gt.1)then
Do ix=1,L
Do iy=1,L
ax1=ix
ay1=iy
ah1(ax1,ay1)=h(ix,iy)
enddo
enddo
mmx=0
mmx1=0
na=0
aa=0
endif
do mm=1,km
n=0
a=0
Do x=1,L
Do y=1,L
an(x,y)=0
enddo
enddo
ix=rand()*L+1
iy=rand()*L+1
t=0
h(ix,iy)=h(ix,iy)+1
if(km.gt.1)then
if(h(ix,iy).le.Z.and.mmx.eq.1) then
kfc=0
kfc=rand()*2+1
if(kfc.eq.2) mmx=0
endif
if(h(ix,iy).le.Z.and.mmx.eq.0) then
do ix=1,L
do iy=1,L
ax2=ix
ay2=iy
ah2(ax2,ay2)=h(ix,iy)
enddo
enddo
mmx=1
endif
endif
if(h(ix,iy).gt.Z) then
n=n+1
an(ix,iy)=an(ix,iy)+1
if(an(ix,iy).eq.1) then
a=a+1
else
a=a+0
end if
h(ix,iy)=h(ix,iy)-Z
t=1
tn(1)=1
tx(1)=ix
ty(1)=iy
101 m=0
DO k=1,tn(t)
ixtmp=tx(k)
iytmp=ty(k)
DO 90 kn=1,Z
call neighbor(noapi,ixtmp,iytmp,kn,L,i,j,nd)
if(nd.eq.1) goto 90
h(i,j)=h(i,j)+1
if(h(i,j).gt.Z) then
n=n+1
an(i,j)=an(i,j)+1
if(an(i,j).eq.1) then
a=a+1
else
a=a+0
end if
m=m+1
tx1(m)=i
ty1(m)=j
h(i,j)=h(i,j)-Z
endif
90 continue
enddo
if(m.ge.1) then
t=t+1
tn(t)=m
DO m1=1,m
tx(m1)=tx1(m1)
ty(m1)=ty1(m1)
enddo
goto 101
endif
endif
if(km.gt.1)then
if(mmx1.eq.1.and.mmx.eq.0) then
if(aa.gt.a.and.mmx1.eq.1)mmx1=0
if(aa.eq.a.and.mmx1.eq.1)then
kfc1=0
kfc1=rand()*2+1
if(kfc1.eq.2) mmx1=0
endif
endif
if(mmx1.eq.0.and.mmx.eq.0) then
do ix=1,L
do iy=1,L
ax3=ix
ay3=iy
ah3(ax3,ay3)=h(ix,iy)
enddo
enddo
na=n
aa=a
mmx1=1
endif
Do ax1=1,L
Do ay1=1,L
ix=ax1
iy=ay1
h(ix,iy)=ah1(ax1,ay1)
enddo
enddo
endif
enddo
if(km.gt.1)then
if(mmx.eq.1) then
do ax2=1,L
do ay2=1,L
ix=ax2
iy=ay2
h(ix,iy)=ah2(ax2,ay2)
enddo
enddo
a=0
n=0
endif
if(mmx1.eq.1.and.mmx.eq.0) then
do ax3=1,L
do ay3=1,L
ix=ax3
iy=ay3
h(ix,iy)=ah3(ax3,ay3)
enddo
enddo
a=aa
n=na
endif
endif
n1=0
n2=0
n3=0
n4=0
DO i1=1,L
DO j1=1,L
if(h(i1,j1).eq.1) n1=n1+1
if(h(i1,j1).eq.2) n2=n2+1
if(h(i1,j1).eq.3) n3=n3+1
if(h(i1,j1).eq.4) n4=n4+1
enddo
enddo
if(2**(0).le.a.and.a.lt.2**(1)) a1=a1+1
if(2**(1).le.a.and.a.lt.2**(2)) a2=a2+1
if(2**(2).le.a.and.a.lt.2**(3)) a3=a3+1
if(2**(3).le.a.and.a.lt.2**(4)) a4=a4+1
if(2**(4).le.a.and.a.lt.2**(5)) a5=a5+1
if(2**(5).le.a.and.a.lt.2**(6)) a6=a6+1
if(2**(6).le.a.and.a.lt.2**(7)) a7=a7+1
if(2**(7).le.a.and.a.lt.2**(8)) a8=a8+1
if(2**(8).le.a.and.a.lt.2**(9)) a9=a9+1
if(2**(9).le.a.and.a.lt.2**(10)) a10=a10+1
if(2**(10).le.a.and.a.lt.2**(11)) a11=a11+1
if(2**(11).le.a.and.a.lt.2**(12)) a12=a12+1
if(2**(12).le.a.and.a.lt.2**(13)) a13=a13+1
if(2**(13).le.a.and.a.lt.2**(14)) a14=a14+1
if(2**(14).le.a.and.a.lt.2**(15)) a15=a15+1
if(2**(15).le.a.and.a.lt.2**(16)) a16=a16+1
if(2**(16).le.a.and.a.lt.2**(17)) a17=a17+1
if(2**(17).le.a.and.a.lt.2**(18)) a18=a18+1
if(2**(0).le.n.and.n.lt.2**(1)) na1=na1+1
if(2**(1).le.n.and.n.lt.2**(2)) na2=na2+1
if(2**(2).le.n.and.n.lt.2**(3)) na3=na3+1
if(2**(3).le.n.and.n.lt.2**(4)) na4=na4+1
if(2**(4).le.n.and.n.lt.2**(5)) na5=na5+1
if(2**(5).le.n.and.n.lt.2**(6)) na6=na6+1
if(2**(6).le.n.and.n.lt.2**(7)) na7=na7+1
if(2**(7).le.n.and.n.lt.2**(8)) na8=na8+1
if(2**(8).le.n.and.n.lt.2**(9)) na9=na9+1
if(2**(9).le.n.and.n.lt.2**(10)) na10=na10+1
if(2**(10).le.n.and.n.lt.2**(11)) na11=na11+1
if(2**(11).le.n.and.n.lt.2**(12)) na12=na12+1
if(2**(12).le.n.and.n.lt.2**(13)) na13=na13+1
if(2**(13).le.n.and.n.lt.2**(14)) na14=na14+1
if(2**(14).le.n.and.n.lt.2**(15)) na15=na15+1
if(2**(15).le.n.and.n.lt.2**(16)) na16=na16+1
if(2**(16).le.n.and.n.lt.2**(17)) na17=na17+1
if(2**(17).le.n.and.n.lt.2**(18)) na18=na18+1
if(2**(18).le.n.and.n.lt.2**(19)) na19=na19+1
if(2**(19).le.n.and.n.lt.2**(20)) na20=na20+1
if(2**(20).le.n.and.n.lt.2**(21)) na21=na21+1
if(2**(21).le.n.and.n.lt.2**(22)) na22=na22+1
if(2**(22).le.n.and.n.lt.2**(23)) na23=na23+1
p1=p1+real(n1)/real(L*L)
p2=p2+real(n2)/real(L*L)
p3=p3+real(n3)/real(L*L)
p4=p4+real(n4)/real(L*L)
enddo
write(*,*) 'n1=',na1,',','a1=',a1
write(*,*) 'n2=',na2,',','a2=',a2
write(*,*) 'n3=',na3,',','a3=',a3
write(*,*) 'n4=',na4,',','a4=',a4
write(*,*) 'n5=',na5,',','a5=',a5
write(*,*) 'n6=',na6,',','a6=',a6
write(*,*) 'n7=',na7,',','a7=',a7
write(*,*) 'n8=',na8,',','a8=',a8
write(*,*) 'n9=',na9,',','a9=',a9
write(*,*) 'n10=',na10,',','a10=',a10
write(*,*) 'n11=',na11,',','a11=',a11
write(*,*) 'n12=',na12,',','a12=',a12
write(*,*) 'n13=',na13,',','a13=',a13
write(*,*) 'n14=',na14,',','a14=',a14
write(*,*) 'n15=',na15,',','a15=',a15
write(*,*) 'n16=',na16,',','a16=',a16
write(*,*) 'n17=',na17,',','a17=',a17
write(*,*) 'n18=',na18,',','a18=',a18
write(*,*) 'n19=',na19
write(*,*) 'n20=',na20
write(*,*) 'n21=',na21
write(*,*) 'n22=',na22
write(*,*) 'n23=',na23
write(*,*) '1的平均',p1/real(noap)
write(*,*) '2的平均',p2/real(noap)
write(*,*) '3的平均',p3/real(noap)
write(*,*) '4的平均',p4/real(noap)
call CPU_time(time_end)
write(*,*) 'time=',time_end-time_begin,'s'
read(*,*)
stop
end
subroutine neighbor(noapi,i1,j1,k,L,i,j,nd)
k=rand()*4+1
nd=0
if(k.eq.1) then
i=i1-1
if(i.eq.0) nd=1
j=j1
endif
if(k.eq.2) then
i=i1
j=j1+1
if(j.eq.L+1) nd=1
endif
if(k.eq.3) then
i=i1+1
if(i.eq.L+1) nd=1
j=j1
endif
if(k.eq.4) then
i=i1
j=j1-1
if(j.eq.0) nd=1
endif
return
end
F 币
贡献
权杖
窥视卡
发表于 2017-4-3 01:42:32
微信
收藏
点赞
点踩
提升卡
置顶卡
变色卡
千斤顶
显身卡
发表于 2017-4-3 11:00:51
楼主
