見出し画像

【Fortranコード】死ぬほどべた書きで書いた「バーンズリーのシダ」のプログラム

物理・プログラミング・Python、Web制作・自己啓発

最近出たこちらの参考書をきっかけにFortranを勉強しています。

思いっきりべた書きで「バーンズリーのシダ」のコードをFortranで書いてみました。あまりきれいで汎用性があるコードではありませんが追々書き換えていくつもりです。

バーンズリーのシダって何？って感じかもしれませんが、下記wikipediaにあるようにフラクタル構造を持った植物です。

これをモデル化してFortranで数値計算させるということです。
理論解説は上記wikipediaに解説があるのでご参照ください。

Fortranコード

program main
    implicit none
    
    integer, parameter :: Nmax=500000, Nx=512, Ny=512
    double precision, parameter :: Lx=6.d0, Ly=6.d0, dx=Lx/dble(Nx), dy=Ly/dble(Ny)
    complex(kind(0d0)), parameter :: uso=(0.d0,1.d0)
    integer :: i, j, n, m, ix, iy
    real(kind=8)::random_value
    double precision:: x, y
    double precision, allocatable, dimension(:) :: x_, y_
    double precision, allocatable, dimension(:,:) :: f

    character(len=40) filename

    call system("mkdir data")

    ! 配列の割付
    allocate(x_(0:Nx-1), y_(0:Ny-1))
    allocate(f(0:Nx-1,0:Ny-1))

    m = 0
    ! x, y座標
    x = 0
    y = 0

    do j=0,Ny-1
        do i=0,Nx-1
       x_(i) = -Lx/2.d0+dx*dble(i)
       y_(j) = -Ly/2.d0+dy*dble(j)
        end do
    end do

    do i=0, Nmax
        call random_seed !乱数の初期値
        call random_number(random_value)

        if (0<= random_value .and. random_value<0.01) then
            x = 0.d0
            y = 0.16d0*y
            ix = int(Nx/2 + x*Nx/10)
            iy = int(y * Ny /12)
            f(ix,iy) = 1.0d0

        else if (0.01d0<= random_value .and. random_value<0.86d0) then
            x = 0.85*x + 0.04d0*y
            y = -0.04d0*x + 0.85d0*y +1.6d0
            ix = int(Nx/2 + x*Nx/10)
            iy = int(y * Ny /12)
            f(ix,iy) = 1.0d0

        else if(0.86d0<= random_value .and. random_value<0.93d0) then
            x = 0.2d0*x - 0.26d0*y
            y = 0.23d0*x + 0.22d0*y +1.6d0
            ix = int(Nx/2 + x*Nx/10)
            iy = int(y * Ny /12)
            f(ix,iy) = 1.0d0

        else
            x =-0.15d0*x + 0.28d0*y
            y = 0.26d0*x + 0.24d0*y +0.44d0
            ix = int(Nx/2 + x*Nx/10)
            iy = int(y * Ny /12)
            f(ix,iy) = 1.0d0
        end if
    end do
!-----------------------------------------------------------------------------------
!output of movie file
    write(filename,"(a,i5.5,a)") "data/barsley_",int(m),".vtk"
    open(10,file=filename)
    write(10,"('# vtk DataFile Version 3.0')")
    write(10,"('test')")
    write(10,"('ASCII ')")
 
    write(10,"('DATASET STRUCTURED_GRID')")
    write(10,"('DIMENSIONS ',3(1x,i3))") Nx, Ny, 1
 
    write(10,"('POINTS ',i9,' float')") Nx*Ny*1
    do j=0,Ny-1
        do i=0,Nx-1
            write(10,"(3(f9.4,1x))") x_(i), y_(j), 0.d0
        end do
    end do
 
    write(10,"('POINT_DATA ',i9)") Nx*Ny*1
 
 !date input
    write(10,"('SCALARS julia float')")
    write(10,"('LOOKUP_TABLE default')")
    do j=0,Ny-1
        do i=0,Nx-1
            write(10,*)  f(i,j)
        end do
    end do
 
    close(10)

    deallocate(x_, y_, f)
end program main

paraviewで可視化するとこのようになります。

可視化にはVTKFortranというライブラリがあるので今後検討してみようと思います。

追ってsubroutineやmodulleといった機能ごとにプログラムを分けて管理しやすいようにして置こうと思います。
そうすることで初期状態、境界条件の設定のみを変えることで他のプログラム部分は意識せずにプログラムを組むことができるということになります。

Twitter➡@t_kun_kamakiri
Instagram➡kamakiri1225
youtube➡https://www.youtube.com/channel/UCbG6_Q9ZRqqVT6YZOpcjDlQ
ブログ➡宇宙に入ったカマキリ（物理ブログ）
ココナラ➡物理の質問サポートサービス
コミュニティ➡製造業ブロガー

この記事が気に入ったらサポートをしてみませんか？