후배가 좀 도와달래서 짜준 프로그램. 출력물을 gnuplot으로 점 찍어보면 뭔가 재미난게 나온다. 뭐하는 프로그램일까요? Giordano의 전산물리 교재에 보면 정답이 있다고 합니다. 저는 그냥 알고리즘만 이해하고 절반정도 도와줬지요.
참고로 포트란 프로그램입니다.
more..
dimension I2(0:101,0:101)
do ii=0,101
do jj=0,101
I2(ii,jj)=0
end do
end do
I2(50,50)=1
DO k=1,100
100 i=2
j=2
if (I2(3,2).EQ.1) then
goto 300
else if (I2(1,2).EQ.1) then
goto 300
else if (I2(2,3).EQ.1) then
goto 300
else if (I2(2,1).EQ.1) then
goto 300
else
101 if (I2(i+1,j).EQ.1) then
goto 100
else if (I2(i-1,j).EQ.1) then
goto 100
else if (I2(i,j+1).EQ.1) then
goto 100
else if (I2(i,j-1).EQ.1) then
goto 100
else
call random_number(a)
c print *,a
if ((0.0.LT.a).and.(a.LT.0.25)) then
if (i.EQ.100) then
I2(1,j)=1
I2(i,j)=0
i=1
else
I2(i+1,j)=1
I2(i,j)=0
i=i+1
end if
else if ((0.25.LT.a).and.(a.LT.0.5)) then
if (i.EQ.1) then
I2(100,j)=1
I2(i,j)=0
i=100
else
I2(i-1,j)=1
I2(i,j)=0
c print *,i
i=i-1
end if
else if ((0.5.LT.a).and.(a.LT.0.75)) then
if (j.EQ.100) then
I2(i,1)=1
I2(i,j)=0
j=1
else
I2(i,j+1)=1
I2(i,j)=0
j=j+1
end if
else if ((0.75.LT.a).and.(a.LT.1.0)) then
if (j.EQ.1) then
I2(i,100)=1
I2(i,j)=0
j=100
else
I2(i,j-1)=1
I2(i,j)=0
j=j-1
end if
end if
end if
goto 101
end if
300 i=i
end do
200 open(1,file=’bbb.dat’)
do i=1,100
do j=1,100
if (I2(i,j).EQ.1) then
write(1,999) i,j
end if
end do
end do
close(1)
999 format (2x,I10,2x,I10)
cccccccc mass versus radius cccccccccccccc
ll=0
open(2,file=’density.dat’)
do n=1,49
do i=50-n,50+n
j=50+n
ll=ll+I2(i,j)
end do
do j=50-n,49+n
i=50+n
ll=ll+I2(i,j)
end do
do i=50-n,49+n
j=50-n
ll=ll+I2(i,j)
end do
do j=49-n,49+n
i=50-n
ll=ll+I2(i,j)
end do
write(2,999) n,ll
end do
close(2)
end
댓글 남기기