fpm に c プログラム入れて X window に図
fpm を利用して unicode plot と X-window で同じ図を出すことにします。
そのためには、fpm に c プログラムを認識させ、X11 のライブラリとリンクさせる必要があります。fpm 公式 gihub に例題集があって、それを見ることで何とかなりました。
fpm.toml の [build] に link = "X11" を付け加えるだけで解決します。
[build] auto-executables = true auto-tests = true auto-examples = true link = "X11"
実行結果
プログラム
OO 式に abstract type で抽象クラスを書いておきました。これによって fig オブジェクトを allocate する時に、どのデバイスに出すかを決定しています。そこ以外は共通のプログラムで書けます。
改訂版
下に貼ったソースコードを改訂したものを github に上げました。
git clone https://github.com/f66blog/xplot cd xplot fpm run fpm test
メインプログラム app/main.f90
program main use device implicit none class (device_t), allocatable :: fig real, allocatable :: x(:), y(:) integer :: n = 10**3 allocate(x(n), y(n)) call random_seed() call random_number(x) call random_number(y) uplot: block use uniplot allocate(fig, source = fig_t(100, 100)) call monte_carlo(fig, x, y) deallocate(fig) end block uplot xplot: block use xplot allocate(fig, source = fig_t(450, 450)) call monte_carlo(fig, x, y) deallocate(fig) end block xplot contains subroutine monte_carlo(fig, x, y) class(device_t), intent(in out) :: fig real, intent(in) :: x(:), y(:) integer :: i, ix0, iy0, ix1, iy1, k k = fig%nx print * print *, 'Monte Carlo: estimated pi =', 4.0 * count(x**2 + y**2 < 1.0) / size(x) call fig%init() ! draw box call fig%line0(0, 0, k-1, 0) call fig%line0(0, 0, 0, k-1) call fig%line0(0, k-1, k-1, k-1) call fig%line0(k-1, 0, k-1, k-1) ! draw 1/4 circle ix0 = 0 iy0 = 0 do ix1 = 0, k - 1 iy1 = k - 1 - int(sqrt(real((k-1)**2 - ix1**2))) call fig%line0(ix0, iy0, ix1, iy1) ix0 = ix1 iy0 = iy1 end do do i = 1, n call fig%point(int(k * x(i)), int(k * y(i))) end do call fig%show() end subroutine monte_carlo end program main
抽象クラス定義 src/device.f90
module device implicit none type, abstract :: device_t integer :: nx, ny contains procedure (device_init) , deferred, pass :: init procedure (device_point), deferred, pass :: point procedure (device_show) , deferred, pass :: show procedure (device_line0), deferred, pass :: line0 procedure (device_line) , deferred, pass :: line end type device_t abstract interface subroutine device_init(fig) import :: device_t class(device_t), intent(in out) :: fig end subroutine device_init subroutine device_point(fig, ix, iy) import :: device_t class(device_t), intent(in out) :: fig integer, intent(in) :: ix, iy end subroutine device_point subroutine device_show(fig) import :: device_t class(device_t), intent(in) :: fig end subroutine device_show subroutine device_line0(fig, ix0, iy0, ix1, iy1) import :: device_t class(device_t), intent(in out) :: fig integer, intent(in) :: ix0, iy0, ix1, iy1 end subroutine device_line0 subroutine device_line(fig, x, y, ipen) import :: device_t class(device_t), intent(in out) :: fig real, intent(in) :: x, y integer, intent(in) :: ipen end subroutine device_line end interface end module device
unicode plot 用モジュール src/uniplot.f90
module uniplot use device implicit none private public :: fig_t type, extends(device_t) :: fig_t private integer, allocatable :: array(:, :) contains procedure :: init procedure :: point procedure :: show procedure :: line0 procedure :: line end type fig_t contains subroutine init(fig) class(fig_t), intent(in out) :: fig allocate(fig%array(0:(fig%nx+1)/2, 0:(fig%ny+3)/4) ) end subroutine init subroutine point(fig, ix, iy) class(fig_t), intent(in out) :: fig integer, intent(in) :: ix, iy integer :: iax, iay iax = ix / 2 iay = iy / 4 ! clipping if (0<=ix .and. ix<fig%nx .and. 0<=iy .and. iy<fig%ny) then fig%array(iax, iay) = ior(fig%array(iax, iay), icode(mod(ix, 2), mod(iy, 4))) end if end subroutine 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(fig, ix0, iy0, ix1, iy1) class(fig_t), intent(in out) :: fig 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 fig%point(ix, iy) else if (abs(nx) < abs(ny)) then d = nx / real(ny) do i = 0, ny, sign(1, ny) call fig%point(nint(ix0 + i * d), iy0 + i) end do else d = ny / real(nx) do i = 0, nx, sign(1, nx) call fig%point(ix0 + i, nint(iy0 + i * d)) end do end if end subroutine line0 subroutine show(fig) class(fig_t), intent(in) :: fig integer :: iy do iy = 0, ubound(fig%array, 2) print '(*(a))', reverse_endian(shift_code(fig%array(:, iy))) end do end subroutine show 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(fig, x, y, ipen) class(fig_t), intent(in out) :: fig 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 fig%line0(ix0, iy0, ix, iy) ix0 = ix iy0 = iy end subroutine line end module uniplot
X window 用 plot モジュール src/xplot.f90
module xplot use device implicit none private public :: fig_t interface subroutine Xopen(nx, ny) bind(c, name = 'X_open') integer, value :: nx, ny end subroutine Xopen subroutine Xpoint(ix, iy) bind(c, name = 'X_point') integer, value :: ix, iy end subroutine Xpoint subroutine Xclose() bind(c, name = 'X_close') end subroutine Xclose subroutine Xflush() bind(c, name = 'X_flush') end subroutine Xflush end interface type, extends(device_t) :: fig_t private contains procedure :: init procedure :: point procedure :: show procedure :: line0 procedure :: line final :: x_close end type fig_t contains subroutine init(fig) class(fig_t), intent(in out) :: fig call Xopen(fig%nx, fig%ny) ! call sleep(1) ! non-standard end subroutine init subroutine point(fig, ix, iy) class(fig_t), intent(in out) :: fig integer, intent(in) :: ix, iy call Xpoint(ix, iy) end subroutine point subroutine show(fig) class(fig_t), intent(in) :: fig call XFlush() end subroutine show subroutine line0(fig, ix0, iy0, ix1, iy1) class(fig_t), intent(in out) :: fig integer, intent(in) :: ix0, iy0, ix1, iy1 integer :: nx, ny, i real :: d nx = ix1 - ix0 ny = iy1 - iy0 if (nx == 0 .and. ny == 0) then call Xpoint(ix0, iy0) else if (abs(nx) < abs(ny)) then d = nx / real(ny) do i = 0, ny, sign(1, ny) call Xpoint(nint(ix0 + i * d), iy0 + i) end do else d = ny / real(nx) do i = 0, nx, sign(1, nx) call Xpoint(ix0 + i, nint(iy0 + i * d)) end do end if end subroutine line0 subroutine x_close(fig) type(fig_t), intent(in out) :: fig print *, 'press ENTER to continue' read * call Xclose() end subroutine x_close subroutine line(fig, x, y, ipen) class(fig_t), intent(in out) :: fig real, intent(in) :: x, y integer, intent(in) :: ipen integer, save :: ix0 = 0, iy0 = 0 integer :: ix, iy real, parameter :: xn = 250.0, yn = 600.0, fx = 4.0, fy = 4.0 ix = nint( fx * x + xn) iy = nint(-fy * y + yn) if (ipen == 1) call fig%line0(ix0, iy0, ix, iy) ix0 = ix iy0 = iy end subroutine line end module xplot
X11 c 言語インターフェース src/plot.c
#include <X11/Xlib.h> #include <X11/Xutil.h> static Display* d; static Window w; static GC gc; unsigned long white, black; void X_open(int nx, int ny){ // Open a display d = XOpenDisplay(0); if ( !d ) return; // white = WhitePixel(d, DefaultScreen(d)); black = BlackPixel(d, DefaultScreen(d)); // Create a window w = XCreateSimpleWindow(d, DefaultRootWindow(d), 0, 0, nx, ny, 0, black, white); XMapWindow(d, w); gc = XCreateGC(d, w, 0, 0); XFlush(d); } void X_point(int ix, int iy){ XDrawPoint(d, w, gc, ix, iy); XFlush(d); } void X_flush(){ XFlush(d); } void X_close(void){ XFreeGC(d, gc); XDestroyWindow(d, w); XFlush(d); XCloseDisplay(d); }
