読者です 読者をやめる 読者になる 読者になる

fortran66のブログ

fortran について書きます。

Mandelbrot BMP 改2

    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  
      !
      type :: t_bmp
        type(t_bmp_file_header) :: file_header
        type(t_bmp_info_header) :: info_header
        type(t_rgb), allocatable :: rgb(:, :)
      contains
        procedure :: wr
      end type t_bmp    
    contains   
      subroutine wr(bmp, fn)
        class(t_bmp), intent(in out) :: bmp
        character(len = *) :: fn
        integer :: nx, ny
        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
        open(9, file = trim(fn)//'.bmp', form = 'binary', status = 'unknown')
        write(9) bmp%file_header, bmp%info_header, bmp%rgb
        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_bmp) :: bmp
      c = xypos(nx = 1280, xmin = -2.0_kd, xmax = 2.0_kd,    ny = 1280, ymin = -2.0_kd, ymax = 2.00_kd)
      bmp%rgb = to_rgb( mandel(c) ) 
      call bmp%wr('mandel')
      c = xypos(nx = 1280, xmin = 0.15_kd, xmax = 0.25_kd,   ny = 1280, ymin = 0.50_kd, ymax = 0.60_kd)
      bmp%rgb = to_rgb( mandel(c) ) 
      call bmp%wr('mandel2')
      stop
    contains
      pure function xypos(nx, ny, xmin, ymin, xmax, ymax)
        integer , intent(in) :: nx, ny
        real(kd), intent(in) :: xmin, ymin, xmax, ymax
        complex(kd), allocatable :: xypos(:, :)
        integer :: ix, iy
        allocate(xypos(nx, ny))
        forall (ix = 1:nx, iy = 1:ny) xypos(ix, iy) = cmplx(pos(ix, nx, xmin, xmax), pos(iy, ny, ymin, ymax)) 
      end function xypos  
      !
      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