program diffusion
implicit none
real,parameter :: x0=-1.0,x1=1.0
real :: dx,dt
integer,parameter :: kosuux=128
integer,parameter :: kosuut=200
integer,parameter :: nint=10
real,dimension(0:kosuux) :: x,ondo
real,dimension(0:kosuux) :: ondo2d
integer :: i,k,n
!
dx = (x1-x0)/real(kosuux) ! x 方向の「小区間」の幅
dt = 0.25*dx*dx ! 時間刻み (x方向の刻み幅の2乗に比例させる必要がある)
do i = 0,kosuux
x(i) = x0+dx*real(i)
end do
do i = 0,kosuux ! 初期条件を設定
if (abs(x(i)) <= 1.0/32.0 ) then ! |x| <= 1/32 の範囲だけ
ondo(i) = 1.0 ! ondo = 1.0 としておく
else
ondo(i) = 0.0 ! その範囲外では ondo = 0.0 にしておく
end if
end do
k = 0 ! お絵描き用サブルーチンを呼び
call oekaki(k,x,ondo) ! 初期条件での ondo の分布を書き出させる
do n = 1,kosuut ! 時間発展を計算するループの始まり
do i = 1,kosuux-1
ondo2d(i) = & ! 現在の時刻での「温度」の
((ondo(i+1)-ondo(i))/dx & ! 2階微分をとる
-(ondo(i)-ondo(i-1))/dx )/dx ! (1階微分をさらに微分する)
end do
do i = 1,kosuux-1 ! 次の時点での「温度」の値を計算し、
ondo(i) = ondo(i)+dt*ondo2d(i) ! 配列 ondo の中身を更新
end do ! (更新の手順はオイラー法に準じる)
ondo(0) = 0.0 ! この2点の値は常に変化しないものとする
ondo(kosuux) = 0.0 ! (要するに「境界条件」)
if (mod(n,nint) == 0) then ! 時間発展ループが nint 回動くたびに
k = n/nint ! お絵描き用サブルーチンを呼び
call oekaki(k,x,ondo) ! その時刻での ondo の分布を書き出させる
end if ! ちなみに mod は「あまり」を求める関数
end do ! 時間発展を計算するループの終わり
stop
contains
subroutine oekaki(k,x,ondo)
integer,intent(in) :: k
real,dimension(0:kosuux),intent(in) :: ondo,x
character(len=2) :: ck2 ! 長さ2の文字型変数
integer :: i
! 出力用ディレクトリがなければ作成
call system('if [ ! -d diffusion ]; then mkdir -p diffusion; fi')
write(ck2,'(i2.2)')k ! 整数 k を文字列に変換 (11 -> 11 とか、3 -> 03 とか)
open(10,file="diffusion/"//ck2//".txt") ! 書き出すためのファイルを開く
do i = 0,kosuux ! ちなみに // は文字列を連結して新しい文字列を作る命令
write(10,'(12f8.4)')x(i),ondo(i)
end do
close(10) ! 書き出しが完了したらファイルを閉じる
end subroutine oekaki
end program diffusion