Logistic map & Lorenz attractor
昔のインターフェースを流用して、作図してみます。
横に縞が入るのが、昔の電送写真のようでいいですね。
プログラム
device.f90
module device implicit none type, abstract :: device_t integer :: nsize_x = 640, nsize_y = 480 character(len = 80) :: title = 'Plotter' integer :: width = 1, color = 0 contains procedure (device_on) , deferred, pass :: on procedure (device_off) , deferred, pass :: off procedure (device_show), deferred, pass :: show procedure (device_pen) , deferred, pass :: pen procedure (device_line), deferred, pass :: line procedure (device_move), deferred, pass :: move end type device_t abstract interface subroutine device_on(self) import :: device_t class(device_t), intent(in out) :: self end subroutine device_on subroutine device_off(self, isec) import :: device_t class(device_t), intent(in out) :: self integer, intent(in), optional :: isec end subroutine device_off subroutine device_show(self) import :: device_t class(device_t), intent(in) :: self end subroutine device_show subroutine device_pen(self, iwidth, icolor) import :: device_t class(device_t), intent(in out) :: self integer, intent(in), optional :: iwidth, icolor end subroutine device_pen subroutine device_line(self, ix, iy) import :: device_t class(device_t), intent(in out) :: self integer, intent(in) :: ix, iy end subroutine device_line subroutine device_move(self, ix, iy) import :: device_t class(device_t), intent(in out) :: self integer, intent(in) :: ix, iy end subroutine device_move end interface end module device
unicode.f90
module unicode use device implicit none private public :: unicode_t type, extends(device_t) :: unicode_t private integer :: ix0 = 0, iy0 = 0 integer, allocatable :: array(:, :) contains procedure, pass :: on => unicode_on procedure, pass :: off => unicode_off procedure, pass :: show => unicode_show procedure, pass :: pen => unicode_pen procedure, pass :: line => unicode_line procedure, pass :: move => unicode_move procedure, pass :: point => unicode_point procedure, pass :: line0 end type unicode_t contains subroutine unicode_on(self) class(unicode_t), intent(in out) :: self allocate(self%array(0:(self%nsize_x+1)/2, 0:(self%nsize_y+3)/4), source = 0 ) end subroutine unicode_on subroutine unicode_off(self, isec) class(unicode_t), intent(in out) :: self integer, intent(in), optional :: isec deallocate(self%array) end subroutine unicode_off subroutine unicode_show(self) class(unicode_t), intent(in) :: self integer :: iy do iy = 0, ubound(self%array, 2) print '(*(a))', reverse_endian(shift_code(self%array(:, iy))) end do end subroutine unicode_show subroutine unicode_pen(self, iwidth, icolor) class(unicode_t), intent(in out) :: self integer, intent(in), optional :: iwidth, icolor return end subroutine unicode_pen subroutine unicode_move(self, ix, iy) class(unicode_t), intent(in out) :: self integer, intent(in) :: ix, iy self%ix0 = ix self%iy0 = iy end subroutine unicode_move subroutine unicode_line(self, ix, iy) class(unicode_t), intent(in out) :: self integer, intent(in) :: ix, iy call self%line0(self%ix0, self%iy0, ix, iy) self%ix0 = ix self%iy0 = iy end subroutine unicode_line subroutine unicode_point(self, ix, iy) class(unicode_t), intent(in out) :: self integer, intent(in) :: ix, iy integer :: iax, iay iax = ix / 2 iay = iy / 4 ! clipping if (0 <= iax .and. iax <= ubound(self%array, 1) .and. & 0 <= iay .and. iay <= ubound(self%array, 2) ) then self%array(iax, iay) = ior(self%array(iax, iay), icode(mod(ix, 2), mod(iy, 4))) end if end subroutine unicode_point pure elemental integer function icode(kx, ky) integer, intent(in) :: kx, ky if (ky == 3) then icode = 64 + 64 * kx else ! 0, 1, 2 icode = 2**(ky + 3*kx) end if end function icode subroutine line0(self, ix0, iy0, ix1, iy1) class(unicode_t), intent(in out) :: self integer, intent(in) :: ix0, iy0, ix1, iy1 integer :: i, ix, iy, nx, ny real :: d nx = ix1 - ix0 ny = iy1 - iy0 if (nx == 0 .and. ny == 0) then call self%point(ix0, iy0) else if (abs(nx) < abs(ny)) then d = nx / real(ny) do i = 0, abs(ny) ix = nint(ix0 + d * sign(i, ny)) iy = iy0 + sign(i, ny) call self%point(ix, iy) end do else d = ny / real(nx) do i = 0, abs(nx) iy = nint(iy0 + d * sign(i, nx)) ix = ix0 + sign(i, nx) call self%point(ix, iy) end do end if end subroutine line0 pure elemental integer function shift_code(k) integer, intent(in) :: k integer, parameter :: n0 = Z'E2A080' !14852224 shift_code = n0 + 256 * (k /64) + mod(k, 64) !E2A180, E2A280, E2A380 end function shift_code pure elemental character(len = 4) function reverse_endian(i) integer, intent(in) :: i character:: tmp(4) tmp = transfer(i, ' ', size = 4) reverse_endian = transfer(tmp(4:1:-1), ' ') !array 4 to len 4 end function reverse_endian subroutine line(self, x, y, ipen) class(unicode_t), intent(in out) :: self real, intent(in) :: x, y integer, intent(in) :: ipen integer, save :: ix0 = 0, iy0 = 0 integer :: ix, iy real, parameter :: xn = 80.0, yn = 100.0, fx = 1.0, fy = 0.85 ix = nint( fx * x + xn) iy = nint(-fy * y + yn) if (ipen == 1) call self%line0(ix0, iy0, ix, iy) ix0 = ix iy0 = iy end subroutine line end module unicode
main.f90
program test use unicode implicit none class (unicode_t), allocatable :: fig chaos: block integer :: ix, i, nx = 150, ny = 80 real :: p, x, y allocate(fig, source = unicode_t(nx, ny, 'Chaos')) print *, 'Chaos' call fig%on() do ix = 1, nx p = 0.3 x = ix * (3.0 - 1.5) / nx + 1.5 do i = 1, 50 p = p + x * p * (1.0 - p) end do do i = 51, 100 y = p / 1.5 * ny call fig%move(ix - 1, ny - INT(y)) call fig%line(ix - 0, ny - INT(y)) p = p + x * p * (1 - p) end do end do call fig%show() call fig%off() deallocate(fig) end block chaos Lorenz_attractor: block integer, parameter :: kd = kind(1.0d0) real (kd) :: x, y, z, dx, dy, dz, a, b, c, d integer :: k, mult = 2, nx = 150, ny = 100 allocate(fig, source = unicode_t(nx, ny, 'Lorenz attractor')) print *, 'Lorenz attractor' call fig%on() a = 10.0d0 b = 28.0d0 c = 8.0d0 / 3.0d0 d = 0.01d0 x = 1.0d0 y = 1.0d0 z = 1.0d0 do k = 1, 3000 dx = a * (y - x) dy = x * (b - z) - y dz = x * y - c * z x = x + d * dx y = y + d * dy z = z + d * dz if (k < 100) then call fig%move( INT(mult * (x + 30)), INT(mult * z) ) else call fig%line( INT(mult * (x + 30)), INT(mult * z) ) end if end do call fig%show() deallocate(fig) end block Lorenz_attractor end program test
