「プログラムでシダを描画する」を Fortran でシダ植物

空前のシダ描画ブーム到来！？(^^;)
あなたも得意なプログラミング言語でシダを描画してみよう！

出遅れてますが、夏休みということで。

Fortran は絵が苦手なので、BMP ファイルを生成して書き込むことにします。利点はライブラリに頼らず、全部自前で書けることです。
###実行結果

BMP は MS-Paint で jpg に直してから upload しました。

###ソースプログラム
注意！RGB の 赤と青が逆転していたので修正しました。(H26.8.20)

    module m_bmp
      implicit none
      type :: t_bmp_file_header
        sequence  
        integer(2) :: bfType = transfer('BM', 0_2, 1) ! BitMap
        integer(4) :: bfSize          ! file size in bytes
        integer(2) :: bfReserved1 = 0 ! always 0
        integer(2) :: bfReserved2 = 0 ! always 0
        integer(4) :: bfOffBits
      end type t_bmp_file_header
      !
      type :: t_bmp_info_header
        sequence
        integer(4) :: biSize     = Z'28' ! size of bmp_info_header ; 40bytes 
        integer(4) :: biWidth
        integer(4) :: biHeight
        integer(2) :: biPlanes   = 1 ! always 1
        integer(2) :: biBitCount
        integer(4) :: biCompression = 0 !0:nocompression,1:8bitRLE,2:4bitRLE,3:bitfield
        integer(4) :: biSizeImage
        integer(4) :: biXPelsPerMeter = 3780 ! 96dpi
        integer(4) :: biYPelsPerMeter = 3780 ! 96dpi 
        integer(4) :: biClrUsed      = 0
        integer(4) :: biClrImportant = 0 
      end type t_bmp_info_header
      !
      type :: t_rgb
        sequence
        character :: b, g, r
      end type t_rgb 
      type :: t_bmp
        type(t_rgb), allocatable :: rgb(:, :)
      contains 
        procedure :: wr => wr_bmp
      end type  
    contains   
      subroutine wr_bmp(bmp, fn)
        class(t_bmp), intent(in) :: bmp
        character(len = *), intent(in) :: fn
        type(t_bmp_file_header) :: bmp_file_header
        type(t_bmp_info_header) :: bmp_info_header
        associate(nx => size(bmp%rgb, 1), ny => size(bmp%rgb, 2))
          bmp_file_header%bfSize      = 14 + 40 + 0 + nx * ny * 3
          bmp_file_header%bfOffBits   = 14 + 40
          bmp_info_header%biWidth     = nx
          bmp_info_header%biHeight    = ny
          bmp_info_header%biBitCount  = 24 
          bmp_info_header%biSizeImage = nx * ny * 3
        end associate
        open(9, file = fn//'.bmp', form = 'binary', status = 'unknown')
        write(9) bmp_file_header
        write(9) bmp_info_header
        write(9) bmp%rgb
        close(9)
        return
      end subroutine wr_bmp  
    end module m_bmp

    program fern
      use m_bmp
      implicit none
      integer, parameter :: n = 20, nx_shida = 500, ny_shida = 500
      type(t_bmp) :: shida
!
!  define statement functions
!      
      real :: w1x, w1y, w2x, w2y, w3x, w3y, w4x, w4y, x, y
      w1x(x, y) =  0.836 * x + 0.044 * y
      w1y(x, y) = -0.044 * x + 0.836 * y + 0.169
      w2x(x, y) = -0.141 * x + 0.302 * y
      w2y(x, y) =  0.302 * x + 0.141 * y + 0.127
      w3x(x, y) =  0.141 * x - 0.302 * y
      w3y(x, y) =  0.302 * x + 0.141 * y + 0.169
      w4x(x, y) =  0.0
      w4y(x, y) =              0.175337 * y
!
      call random_seed()
      allocate(shida%rgb(nx_shida, ny_shida))
      call f(n, 0.0, 0.0) 
      call shida%wr('shida')
      stop
    contains
      subroutine point(x, y)
        real, intent(in) :: x, y
        shida%rgb(nint(x), nint(y) + 1) = t_rgb(achar(0), achar(255), achar(0))
      end subroutine point  
    
      recursive subroutine f(k, x, y)
        integer, intent(in) :: k
        real, intent(in) :: x, y
        real :: r
        if (k > 0) then
          call f(k - 1, w1x(x, y), w1y(x, y))
          call random_number(r)
          if (r < 0.3) call f(k - 1, w2x(x, y), w2y(x, y)) 
          call random_number(r)
          if (r < 0.3) call f(k - 1, w3x(x, y), w3y(x, y)) 
          call random_number(r)
          if (r < 0.3) call f(k - 1, w4x(x, y), w4y(x, y)) 
        else
          call point(x * nx_shida * 0.98 + 0.5 * nx_shida, y * ny_shida * 0.98)
        end if    
      end subroutine f
    end program fern

bmp 形式は y 軸の向きが普通と上下逆になっているので、適宜修正してあります。しかし、その副作用で値０が出て配列をはみ出すので１ドット分ずらしてあります。

廃止予定項目ですが文関数も使ってみました。