グラフ作画ルーチン（仮

本番の清書などの場合、適切な作図ソフトで作図する分には不満は無いのですが、計算結果のチェックやパラメータ振りのような捨て作図のときには、計算プログラム内で簡単にスコッと描いてくれると便利なことが多いです。計算結果をスクリプトなどで処理してもいいのですが、window がポップしたり埋もれて煩わしかったり、作図ソフト側がイマイチ気の利かないこともままあります。

そういうわけで本体のメインルーチンに、デバッグ行的に軽く挿入して作図してくれるルーチンが前々から欲しかったのですが、Fortran2008 で block construct が導入されて、一気にやりやすくなりました。block 内では use 文によるモジュール読み込みや、ローカル変数の利用が可能なので、局所的に完全に独立したプログラム構造が構成可能だからです。

さて作図ルーチンは今のところ気の利かなさという点では事態はむしろ悪化してるんですが、とりあえず備忘録代わりに書きつけておきます。

プログラム

ここでは、エラトステネスのふるいで素数を求めて、素数の逆数の和を求めます。素数の逆数の和は、LogLog に比例して発散するらしいので、それをグラフによって確かめます。整数N以下の素数を求めるとして、横軸をlog10(log10(N)) にとり、縦軸を N 以下の素数の和に取ります。

グラフを見ると N=10^9 まで直線的に伸びているので、この範囲で確かに loglog で発散していっていることが分かります。最初の方で直線の傾きがずれているのでは、最初の点が素数の数が少なすぎるためだと思われます。

10**9 以下の素数を求めて、10 のべき乗毎に、そこまでの素数の逆数の和を求めています。10**9 以下の場合、約五千万個の素数の逆数の和を求めていることになります。

メイン・プログラム

    module m_eratos
      implicit none
    contains
      function eratos(n) result(ires)
        integer, intent(in) :: n
        integer, allocatable :: ires(:)
        logical(1), allocatable :: tab(:)
        integer :: i, k
        allocate(tab(n))
        tab = .true.
        tab(1) = .false.
        do i = 1, int(sqrt(real(n)))
          if (tab(i)) tab(i**2::i) = .false.  
        end do
        allocate(ires(count(tab)))
        k = 0
        do i = 1, n
          if (.not. tab(i)) cycle
          k = k + 1
          ires(k) = i
        end do            
      end function eratos
    end module m_eratos
    
    program Console6
      use m_eratos
      implicit none
      integer, allocatable :: iprim(:)
      integer :: i, np(9)
      real(8) :: tw(9)
      do i = 1, 9
        iprim = eratos(10**i)     
        np(i) = size(iprim)
        tw(i) = sum(1.0d0 / iprim)
        print *, i, np(i), tw(i)
      end do

      block
        use m_graph
        call plot(log10(log10(real(np))), real(tw), &
                  yrange = [0.0, 4.0], &
                  title = ['Log\_1\_0Log\_1\_0(N)', 'Sum of 1/primes', '\gR1/p~log(log(N))'], &
                  fmt = ['f4.1', 'f3.1'])
      end block    
    end program Console6

グラフ作図用モジュール

x 軸と y 軸の作図を深く考えず共通化しようとしてドツボにはまったので、あきらめて同じようなルーチンを並べています。

    module m_axis
      use m_plot
      implicit none
      real, parameter :: pi = 4 * atan(1.0), rad = pi / 180.0 ! degree to radian
      type, extends (t_fig) :: t_graph
        real :: xmin, xmax, ymin, ymax
      contains
        procedure :: graph
        procedure :: x_axis
        procedure :: y_axis
        procedure :: x_mk_axis
        procedure :: y_mk_axis
        procedure :: x_tick
        procedure :: y_tick
      end type t_graph    
    contains  
      subroutine graph(self, xmin, ymin, xmax, ymax)
        class (t_graph), intent(in out) :: self
        real, intent(in) :: xmin, ymin, xmax, ymax
        self%xmin = xmin
        self%ymin = ymin
        self%xmax = xmax
        self%ymax = ymax         
      end subroutine graph
      
      subroutine x_axis(self, fmt, xstep, mtick, d, title, log)
        class (t_graph), intent(in) :: self
        character (len = *), intent(in) :: fmt
        real   , intent(in) :: xstep, d
        integer, intent(in) :: mtick
        character (len = *), intent(in), optional :: title
        logical, intent(in), optional :: log
        character (len = :), allocatable :: title_
        logical :: log_
        real :: d_scaled
        call self%move(self%xmin, self%ymin)
        call self%line(self%xmax, self%ymin)
        title_ = 'x-axis'
        if (present(title)) title_ = trim(title)
        log_ = .false.
        if (present(log)) log_ = log
        d_scaled = d * (self%y1 - self%y0) / 100.0
        call self%x_mk_axis(fmt, self%xmin, self%xmax, xstep, mtick, d_scaled, title_, log_)
      end subroutine x_axis
      
      subroutine y_axis(self, fmt, ystep, mtick, d, title, log)
        class (t_graph), intent(in) :: self
        character (len = *), intent(in) :: fmt
        real   , intent(in) :: ystep, d
        integer, intent(in) :: mtick
        character (len = *), intent(in), optional :: title
        logical, intent(in), optional :: log
        character (len = :), allocatable :: title_
        logical :: log_
        real :: d_scaled
        call self%move(self%xmin, self%ymin)
        call self%line(self%xmin, self%ymax)
        title_ = 'y-axis'
        if (present(title)) title_ = trim(title)
        log_ = .false.
        if (present(log)) log_ = log
        d_scaled = d * (self%x1 - self%x0) / 100.0
        call self%y_mk_axis(fmt, self%ymin, self%ymax, ystep, mtick, d_scaled, title_, log_)
      end subroutine y_axis

      subroutine x_mk_axis(self, fmt, rmin, rmax, rstep, mtick, d, title, log)
        class (t_graph), intent(in) :: self
        character (len = *), intent(in) :: fmt, title
        real   , intent(in) :: rmin, rmax, rstep, d
        integer, intent(in) :: mtick
        real, allocatable :: r(:)
        logical, intent(in), optional :: log
        logical :: log_
        real :: tmin, xt, yt
        character (len = 20), allocatable :: tr(:)
        integer :: i, j, n
        tmin = rstep * int(rmin / rstep) - rstep
        n = (rmax - tmin) / rstep + 1
        ! tick    
        log_ = .false.
        if (present(log)) log_ = log
        allocate(r(n), tr(n))
        call mk_rt(fmt, tmin, rstep, r, tr, log_)
        call self%x_tick(r, rmin, rmax, 2.0 * d, 'center', tr)
        ! subtick
        r = [(( tmin + rstep * (i - 1) + rstep / mtick * j, j = 1, mtick), i = 0, n)]  
        if (log_) r = [((( tmin + rstep * (i - 1) + rstep * log10(real(j))), j = 1, mtick), i = 0, n)]  
        call self%x_tick(r, rmin, rmax, d)
        ! title
        xt = (self%xmax + self%xmin) / 2 
        yt = self%ymin + 8 * d        
        call self%text(xt, yt, title, 4.0, 0.0, 'center')
      end subroutine x_mk_axis
      
      subroutine y_mk_axis(self, fmt, rmin, rmax, rstep, mtick, d, title, log)
        class (t_graph), intent(in) :: self
        character (len = *), intent(in) :: fmt, title
        real   , intent(in) :: rmin, rmax, rstep, d
        integer, intent(in) :: mtick
        logical, intent(in), optional :: log
        logical :: log_
        real, allocatable :: r(:)
        real :: tmin, xt, yt
        character (len = 20), allocatable :: tr(:)
        integer :: i, j, n
        tmin = rstep * int(rmin / rstep) - rstep 
        n = (rmax - tmin) / rstep + 1   
        ! tick     
        log_ = .false.
        if (present(log)) log_ = log
        allocate(r(n), tr(n))
        call mk_rt(fmt, tmin, rstep, r, tr, log_)
        call self%y_tick(r, rmin, rmax, 2.0 * d, 'right', tr)
        ! subtick
        r = [(( tmin + rstep * (i - 1) + rstep / mtick * j, j = 1, mtick), i = 0, n)]  
        if (log_) r = [((( tmin + rstep * (i - 1) + rstep * log10(real(j))), j = 1, mtick), i = 0, n)]  
        call self%y_tick(r, rmin, rmax, d)
        ! title 
        xt =  self%xmin - (5.0 + 0.4 * maxval(len_trim(tr))) * d
        yt = (self%ymax + self%ymin) / 2 
        call self%text(xt, yt, title, 4.0, 90.0, 'center')
      end subroutine y_mk_axis
      
      subroutine mk_rt(fmt, tmin, rstep, r, tr, log)
        character (len = *), intent(in) :: fmt
        real, intent(in ) :: tmin, rstep
        real, intent(out) :: r(:)
        character (len = *), intent(out) :: tr(:)
        logical, intent(in) :: log
        integer :: i
        real :: t
        do i = 1, size(r) 
          r(i) = tmin + rstep * (i - 1)
          t = r(i)
          if (log) t = 10.0**(r(i))
          if (scan(fmt, 'iI') /= 0) then
            write(tr(i), fmt) int(t)  
          else  
            write(tr(i), fmt) t
          end if 
        end do   
      end subroutine mk_rt

      subroutine x_tick(self, r, rmin, rmax, d, adjust, tr)
        class (t_graph), intent(in) :: self
        real, intent(in) :: r(:), rmin, rmax, d
        character (len = *), intent(in), optional :: adjust, tr(:)
        integer :: i
        real :: r0(2), r1(2), rt(2)
        do i = 1, size(r)
          if (rmin <= r(i) .and. r(i) <= rmax) then
            call self%move(r(i), self%ymin    )
            call self%line(r(i), self%ymin + d)
            if (present(tr)) then
              call self%text(r(i) + 0.01 * (self%xmax - self%xmin), self%ymin + 2.2 * d, trim(tr(i)), 2.2, 0.0, trim(adjust))
            end if  
          end if  
        end do    
      end subroutine x_tick

      subroutine y_tick(self, r, rmin, rmax, d, adjust, tr)
        class (t_graph), intent(in) :: self
        real, intent(in) :: r(:), rmin, rmax, d
        character (len = *), intent(in), optional :: adjust, tr(:)
        integer :: i
        real :: r0(2), r1(2), rt(2)
        do i = 1, size(r)
          if (rmin <= r(i) .and. r(i) <= rmax) then
            call self%move(self%xmin    , r(i))
            call self%line(self%xmin - d, r(i))
            if (present(tr)) then
              call self%text(self%xmin - 0.08 * (self%xmax - self%xmin), r(i) - 0.02 * (self%ymax - self%ymin), trim(tr(i)), 2.2, 0.0, trim(adjust))
            end if  
          end if  
        end do    
      end subroutine y_tick
    end module m_axis
    !================================================================
    module m_graph
      use m_axis
      implicit none
    contains
      subroutine plot(x, y, xrange, yrange, title, fmt, logq)
        real, intent(in) :: x(:)
        real, intent(in), optional :: y(:), xrange(:), yrange(:)
        character (len = *), intent(in), optional :: title(3), fmt(2)
        logical, intent(in), optional :: logq(2)
        real, allocatable :: xx(:), yy(:)
        real :: x0, y0, x1, y1, xmin, ymin, xmax, ymax, xstep, ystep, f
        integer :: i, m, mstepx, mstepy
        character (len = :), allocatable :: fmt_(:)  
        type (t_graph) :: fig
        if (present(y)) then
          xx = x
          yy = y
        else
          allocate(xx(0:size(x)))
          forall (i = 0:size(x)) xx(i) = i  
          allocate(yy(0:size(x)))
          yy(0) = 0.0 
          yy(1:) = x
        end if
        
        xmin = minval(xx)
        ymin = minval(yy) 
        xmax = maxval(xx) 
        ymax = maxval(yy) 
        ymin = ymin - 0.001 * abs(ymax - ymin)
        ymax = ymax + 0.001 * abs(ymax - ymin)
        if (present(xrange)) then 
          xmin = xrange(1)
          xmax = xrange(2)
        end if  
        if (present(yrange)) then 
          ymin = yrange(1)
          ymax = yrange(2)
        end if  
        x0 = xmin - 0.4 * abs(xmax - xmin)
        y0 = ymin - 0.4 * abs(ymax - ymin)
        x1 = xmax + 0.2 * abs(xmax - xmin)
        y1 = ymax + 0.2 * abs(ymax - ymin )
        
        xstep = (xmax - xmin) / 4
        m = floor(log10(xstep))
        if (m==0) m = m - 1
        f = 10.0**m
        xstep = int(xstep / f) * f
        ystep = (ymax - ymin) / 2
        m = floor(log10(ystep))
        if (m==0) m = m - 1
        f = 10.0**m
        ystep = int(ystep / f) * f
        
        call fig%on(t_win32(800, 600, title = 'Graph Auto Plot'))
        call fig%window(x0, y0, x1, y1)
        call fig%graph (xmin, ymin, xmax, ymax)       
        call fig%pen(2)
      
        mstepx = 5
        mstepy = 5
        if (present(logq)) then
          if (logq(1)) mstepx = 10
          if (logq(2)) mstepy = 10
        end if  
        fmt_ = ['(F10.1)', '(ES8.1)'] 
        if (present(fmt)) fmt_ = '('//fmt//')'
        call fig%x_axis(fmt_(1), xstep, mstepx, -2.0, title(1), logq(1))
        call fig%y_axis(fmt_(2), ystep, mstepy,  2.0, title(2), logq(2))
 
        call fig%move(fig%xmin, fig%ymin)
        do i = 1, ubound(xx, 1)
          if (i == 1) call fig%move(xx(i), yy(i))
          call fig%line(xx(i), yy(i))
        end do
        if (present(title)) call fig%text((fig%xmax + fig%xmin) / 2, fig%ymax + 0.02 * abs(fig%ymax - fig%ymin), title(3), 4.0, 0.0, 'center')
        call fig%off()
      end subroutine plot
    end module m_graph