Fortranから使うPLplot入門 #2

Author: 雨崎しのぶ

Twitter: @amasaki203

Posted on: 2024-12-16 JST

1. Fortranから使うPLplot入門 #1 - インストールとHello World
2. Fortranから使うPLplot入門 #2 - 直線と曲線、散布図を描く（この記事）

はじめに

Qiitaに挙げた初回の記事に引き続きPLplotの使い方について解説します。本稿では直線や曲線と、散布図の描き方について紹介します。

前回は、以下のメインプログラムでHello Worldを描きました。このコードのframeブロックとplotブロックに注目して、線分や曲線を引く方法についてみていきましょう。

program main
   use plplot
   implicit none

   initialize: block
      call plinit
   end block initialize

   frame: block
      real(PLFLT), parameter :: xmin = 0d0, xmax = 1d0
      real(PLFLT), parameter :: ymin = -1d0, ymax = 1d0
      integer :: just, axis
      
      just = 0
      axis = 0
      call plenv( xmin, xmax, ymin, ymax, just, axis)
      call pllab("x", "y", "Title")
   end block frame

   plot: block
      real(PLFLT) :: cx, cy, dx, dy
      
      call plschr(12._plflt, 1._plflt)

      cx = 0.5d0
      cy = 0d0
      dx = 1d0
      dy = 0.5d0
      call plptex(cx, cy, dx, dy, 0.5_plflt, "Hello World!")
   end block plot

   finalize: block
      call plend
   end block finalize
 
end program main

目次

• はじめに
• 線を引く
  • 線分を引く
  • 曲線を描く
    • サインカーブ
    • コサインを追加
• 散布図を描く
• まとめ
• 参考文献

線を引く

線分を引く

線を引くためには、plotブロックにplline手続を呼ぶ行を追加します。

pllineは、x-y平面上の複数の点を線分で繋いだものを書き出します。引数には実数型ランク1の配列を2つとります。配列の添字が同じになる2つの値を、点のx座標とy座標の値とみなして、2つの点を結びます。例えばcall plline(x(:), y(:))とすると、最初の点の座標は(x(1), y(1))、2番目の点は(x(2), y(2))というようになります。

plot: block
   real(PLFLT) :: cx, cy, dx, dy
   call plschr(12._plflt, 1._plflt)

   cx = 0.5d0
   cy = 0d0
   dx = 1d0
   dy = 0.5d0
   call plptex(cx, cy, dx, dy, 0.5_plflt, "Hello World!")
   
   call plline([0d0,dx], [0d0,dy])  
    ! ここでは線分を引くためにサイズ2の配列を与えている。
end block plot

このコードを実行すると、次のような図が出力されます。

この図では原点と点(1, 0.5)に線が引かれていますが、このままでは線分であることが明確ではないのでframeブロックを少し変更します。

 frame: block
    real(PLFLT), parameter :: xmin = -0.1d0, xmax = 1.1d0 ! xminとxmaxを変更
    real(PLFLT), parameter :: ymin = -1d0, ymax = 1d0
    integer :: just, axis

    just = 0
    axis = 0
    call plenv( xmin, xmax, ymin, ymax, just, axis)
    call pllab("x", "y", "Title")
 end block frame

これを反映したプログラムを実行すると、次のようになります。

これで線分を引くことができました。これに加えて、前回の記事で 「plptex手続では、文字列を置く座標をcx = 0.5d0; cy = 0d0とし、その傾きを原点と点dx = 1d0; dy = 0.5d0のなす直線に並行になるように指定しています」 と書きましたが、この図はまさに、文字列の傾きが原点と点(dx, dy)を結ぶ直線に並行であることを示しています。

曲線を描く

サインカーブ

前節ではpllineを使って線分を引きました。この手続は、与える配列のサイズを十分大きいものにすると、曲線を引くことができます。

まずは簡単なSin曲線を引いてみましょう。まずはframeブロックのxmin, xmax, ymin, ymaxをグローバル定数に移し、xmaxには2πを代入します。

program main
   use plplot
   implcit none
   real(PLFLT), parameter :: pi = acos(-1d0)
   real(PLFLT), parameter :: xmin = 0d0, xmax = 2*pi
   real(PLFLT), parameter :: ymin = -1d0, ymax = 1d0
   ...

   frame: block
      integer :: just, axis

      just = 0
      axis = 0
      call plenv( xmin, xmax, ymin, ymax, just, axis)
      call pllab("x", "y", "Title")
   end block frame
   ...

次にplotブロックを次のように変更します。

plot: block
   integer, parameter :: n = 101
   integer, :: i
	 real(PLFLT) :: x(n), y_sine(n)
	 
	 do i = 1, n
      x(i) = (i-1)*(xmax - xmin)/(n-1)
      y_sine(i) = sin(x(i))
   end do
   call plline(x, y_sine)
end block plot

配列xには、x軸方向の区間をn-1個に分割したものの境界（n個）の座標値を代入しています。配列yの第i要素には、x座標上の位置x(i)のSin関数の値を代入しています。

以上の変更を行い、プログラムを実行すると、次のような図が出力されるでしょう。

コサインを追加

次にCos関数の曲線を追加しますが、ここでは線の色を変えてみましょう。plotブロックを以下のように変更します。

 plot: block
    integer, parameter :: n = 101
    integer :: i
    real(PLFLT) :: x(n), y_sine(n), y_cosine(n)

    do i = 1, n
       x(i) = (i-1)*(xmax - xmin)/(n-1)
       y_sine(i) = sin(x(i))
       y_cosine(i) = cos(x(i))
    end do
    call plline(x, y_sine)

    call plcol0(2)
    call plline(x, y_cosine)
 end block plot

新しいこのコードでは、配列y_cosineを追加しています。8行目でこの配列にCos関数の値を代入しています。次にplcol0(2)を呼び出して、その後に描画する色を変更しています。最後にplline(x, y_cosine)を呼び出して曲線を描画しています。これを実行すると以下のようなプロットを得られます。

プロットする色を変更するplcol0手続について解説します。このサブルーチンは0以上の整数値を引数に取り、その手続を呼び出して以降に描画するものの色を番号で対応づけられたものに変更します。番号と色の対応はcmap0_default.palというファイルにHEXカラーコードで定義されています。以下のコマンドは、macOS上でMacPortsを使ってPLplotをインストールしている場合の例です。

% ls /opt/local/share/plplot5.15.0 | grep .pal
cmap0_alternate.pal
cmap0_black_on_white.pal
cmap0_default.pal
cmap0_white_bg.pal
cmap1_blue_red.pal
cmap1_blue_yellow.pal
cmap1_default.pal
cmap1_gray.pal
cmap1_highfreq.pal
cmap1_lowfreq.pal
cmap1_radar.pal

% cat /opt/local/share/plplot5.15.0/cmap0_default.pal
16
#000000
#ff0000
#ffff00
#00ff00
#7fffd4
#ffc0cb
#f5deb3
#bebebe
#a52a2a
#0000ff
#8a2be2
#00ffff
#40e0d0
#ff00ff
#fa8072
#ffffff

call plcol0(2)の2に対応するのは、cmap0_default.palにリストされている値の3番目（0-basedのため）、つまり#ffff00なので黄色に変更しています。

散布図を描く

前節までに、文字列を書くことと線を引くことができるようになりました。この節では、点を打って散布図を描く方法について解説します。ここでは、次のような図をプロットすることを目標にします。

プロットするデータの例として、国際天文学連合のMinor Planet CenterからMOBCORB.DATのデータをダウンロードして使用します。このデータファイルは小天体の軌道に関するデータなどが含まれています。小惑星は楕円軌道で運行しているので、軌道を特徴づけるパラメーターの一部である軌道長半径aと軌道離心率eの値を取り出します。

ダウンロードできたら、まずはこのファイルをFortranで読み取りやすいように加工します（完全なデータフォーマットはExport Format for Minor-Planet Orbits - IAU/MPCに書かれています）。ファイルには説明文とヘッダが含まれているので、以下のコマンドを実行してこれを取り除きます（もしくはテキストエディタでデータ部分以外を切り取ってもかまいません）。

% cat MPCORB.DAT | awk 'NR>= 44' > MPCORB_modified.DAT

次にデータを読み込んでプロットを描くコードについて説明していきます。プログラムの先頭は以下のように記述します。

program main
   use plplot
   implicit none
   real(PLFLT), parameter :: xmin = 0d0, xmax = 6d0
   real(PLFLT), parameter :: ymin = 0d0, ymax = 1d0
   integer, parameter :: n = 1500000
   real(PLFLT) :: x(n), y(n)

   ! 値をプロット範囲外の値で初期化する。
   x(:) = -1d0
   y(:) = -1d0

そして、次のようなread_dataブロックを追加して、MPCORB_modified.DATからデータを読み込みます。

 read_data: block
    character(7) :: id
    character(5) :: epoch
    real(PLFLT) :: h, g, perihelion, epoch_ano, lon_asc, inclination, ecc, mean_motion, a
    	! 使用するのはeccとaのみで、他はダミー変数
    integer :: uni, i, ierr

    open(newunit=uni, file='MPCORB_modified.DAT', status='old')

    i = 1
    ierr = 0
    do while (.true.)
       if (ierr /= 0) then
          exit
       end if
       read(uni, *, iostat=ierr)  &
          id, h, g, epoch, epoch_ano, perihelion, &
          lon_asc, inclination, ecc, mean_motion, a
          ! 1レコードには、第9列に軌道離心率eが、第11列に軌道長半径aが配置されているので、
          ! これらを変数eccとaに読み込む

       x(i) = a
       y(i) = ecc
       i = i + 1
    end do
    close(uni)
 end block read_data

次にinitialize, frame, plotブロックを次のように記述します。

initialize: block
   call plspal0('cmap0_white_bg.pal')  ! このパレットファイルを使用して、背景色を白にする
   call plspage(0._plflt, 0._plflt, 1440, 1080, 0, 0)
      ! 画像の大きさを1440x1080ピクセルに設定する
   call plinit
end block initialize

frame: block
   integer :: just, axis

   just = 0
   axis = 0
   call plcol0(15) ! 色を黒に切り替える（座標軸とラベルのため）
   call plenv( xmin, xmax, ymin, ymax, just, axis)
   call pllab("Semimajor Axis (AU)", "Eccentricity", "Axis-Eccentricity Plot")
end block frame

plot: block
   integer :: code

   code = -1
   call plcol0(9) ! 色を青に切り替える（プロットする点のため）
   call plpoin(x, y, code)
end block plot

plotブロックでは、plpoin手続を使用して点をプロットする処理を実行しています。plpoinに配列x, yを与えるのはpllineと同様ですが、引数としてもう一つ、点として打つグリフ（字体）を指定する値を第3引数に渡す必要があります（ここではcode）。変数codeは-1以上127以下の値を指定します。上のコードではcode = -1としていますが、これは単に点を描画するだけを意味します。そのほかにプロット可能なグリフは公式サイトのExample 6で見ることができます。

なお、plspal0, plspage の手続に関しては次回以降の記事で解説したいと思います。

以下にプログラム全体を示します。

program main
   use plplot
   implicit none
   real(PLFLT), parameter :: xmin = 0d0, xmax = 6d0
   real(PLFLT), parameter :: ymin = 0d0, ymax = 1d0

   integer, parameter :: n = 1500000
   real(PLFLT) :: x(n), y(n)

   x(:) = -1d0
   y(:) = -1d0

   read_data: block
      character(7) :: id
      character(5) :: epoch
      real(PLFLT) :: h, g, perihelion, epoch_ano, lon_asc, inclination, ecc, mean_motion, a
      integer :: uni, i, ierr

      open(newunit=uni, file='MPCORB_modified.DAT', status='old')
   
      i = 1
      ierr = 0
      do while (.true.)
         if (ierr /= 0) then
            exit
         end if
         read(uni, *, iostat=ierr)  &
            id, h, g, epoch, epoch_ano, perihelion, &
            lon_asc, inclination, ecc, mean_motion, a

         x(i) = a
         y(i) = ecc
         i = i + 1
      end do
      close(uni)
   end block read_data

   initialize: block
      call plspal0('cmap0_white_bg.pal')  ! このパレットファイルを使用して、背景色を白にする
      call plspage(0._plflt, 0._plflt, 1440, 1080, 0, 0) 
         ! 画像の大きさを1440x1080ピクセルに設定する
      call plinit
   end block initialize

   frame: block
      integer :: just, axis
   
      just = 0
      axis = 0
      call plcol0(15) ! 色を黒に切り替える（座標軸とラベルのため）
      call plenv( xmin, xmax, ymin, ymax, just, axis)
      call pllab("Semimajor Axis (AU)", "Eccentricity", "Axis-Eccentricity Plot")
   end block frame

   plot: block
      integer :: code
      
      code = -1
      call plcol0(9) ! 色を青に切り替える（プロットする点のため）
      call plpoin(x, y, code)
   end block plot

   finalize: block
      call plend
   end block finalize

end program main

このコードを実行すると、この節の冒頭で示したような画像が出力されます。これで散布図を描くことができるようになりました。

まとめ

今回はPLplotを使って、直線と曲線、散布図の書き方について説明しました。線を引くにはpllineを使い、点を打つにはplpoinを使いますが、x軸とy軸の値を配列としてまとめて引数として渡すことは共通しています。座標軸や背景色等を変更する方法の詳細にはほとんど触れていませんが、次回以降の記事で解説したいと思います。

参考文献

1. IAU/Minor Planet Center
2. Export Format for Minor-Planet Orbits - IAU/MPC
3. PLplot - Examples: Example 6
4. Documentation of the PLplot plotting software - v5.15.0: plline, plpoin, plspal0, plspage