fortran66のブログ

fortran について書きます。

Mandelbrot BMP 改

ISO_FORTRAN_ENV で KIND を指定

拡大倍率をあげられるように実数は倍精度にしました。

f:id:fortran66:20140209015126j:plain

ソース・プログラム

戯れに RETURN 文を省略してみました。本来は付ける派閥なんですが・・・

    module m_bmp
      use, intrinsic :: iso_fortran_env
      implicit none
      type :: t_bmp_file_header
        sequence  
        integer(int16) :: bfType = transfer('BM', 0_int16) ! BitMap
        integer(int32) :: bfSize          ! file size in bytes
        integer(int16) :: bfReserved1 = 0 ! always 0
        integer(int16) :: bfReserved2 = 0 ! always 0
        integer(int32) :: bfOffBits
      end type t_bmp_file_header
      !
      type :: t_bmp_info_header
        sequence
        integer(int32) :: biSize     = Z'28' ! size of bmp_info_header ; 40bytes 
        integer(int32) :: biWidth
        integer(int32) :: biHeight
        integer(int16) :: biPlanes   = 1 ! always 1
        integer(int16) :: biBitCount
        integer(int32) :: biCompression = 0 ! 0:nocompression, 1:8bitRLE, 2:4bitRLE, 3:bitfield
        integer(int32) :: biSizeImage
        integer(int32) :: biXPelsPerMeter = 3780 ! 96dpi
        integer(int32) :: biYPelsPerMeter = 3780 ! 96dpi 
        integer(int32) :: biClrUsed       = 0
        integer(int32) :: biClrImportant  = 0 
      end type t_bmp_info_header
      !
      type :: t_rgb
        sequence
        integer(int8) :: ib, ig, ir
      end type t_rgb  
    contains   
      subroutine wr(bmp)
        type(t_rgb), intent(in) :: bmp(:, :)
        type(t_bmp_file_header) :: bmp_file_header
        type(t_bmp_info_header) :: bmp_info_header
        integer :: nx, ny
        nx = size(bmp, 1)
        ny = size(bmp, 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
        open(9, file = 'mandel.bmp', form = 'binary', status = 'unknown')
        write(9) bmp_file_header, bmp_info_header, bmp
        close(9)
      end subroutine wr  
    end module m_bmp
    
    module m_mandel
      use, intrinsic :: iso_fortran_env
      implicit none
      integer, parameter :: kd = real64
      integer, parameter :: maxiter = 255
    contains
      pure elemental integer function mandel(c)
        complex(kd), intent(in) :: c
        complex(kd) :: z
        z = c
        do mandel = maxiter, 1, -1
          if (abs(z) > 2.0_kd) exit
          z = z**2 + c
        end do 
      end function mandel
    end module m_mandel
    
    program BMP1
      use m_bmp 
      use m_mandel
      implicit none
      complex(kd), allocatable :: c(:, :)
      type(t_rgb), allocatable :: bmp(:, :)
      integer  :: ix, iy, nx, ny
      real(kd) :: xmin, ymin, xmax, ymax
      xmin = 0.15_kd
      xmax = 0.25_kd
      ymin = 0.50_kd
      ymax = 0.60_kd
      nx = 1280
      ny = 1280
      allocate( c(nx, ny) )
      forall (ix = 1:nx, iy = 1:ny) c(ix, iy) = cmplx(pos(ix, nx, xmin, xmax), pos(iy, ny, ymin, ymax)) 
      bmp = to_rgb( mandel(c) ) 
      call wr(bmp)
      stop
    contains
      pure real(kd) function pos(i, n, rmin, rmax) 
        integer , intent(in) :: i, n
        real(kd), intent(in) :: rmin, rmax
        pos = (rmax - rmin) / (n - 1) * i + rmin
      end function pos
      !
      pure elemental type(t_rgb) function to_rgb(k)
        integer, intent(in) :: k
        to_rgb = t_rgb(10 * k,  k, 5 * k)
      end function to_rgb
    end program BMP1