GitHubで公開しているC++用GnuplotライブラリADAPT-GPM2に等高線表示機能を追加した。Gnuplotが持つ仕様上の問題を吸収するために、ちょっと回りくどい実装になってしまった。 ライブラリについては過去記事を参照。
Gnuplotのpm3dにおける最大の悩みは、格子点上のデータを与えた時に、ある四角形の色がその四角形の頂点4つを参照して計算されるところだ。デフォルトではあろうことが4点のz座標の平均値で色を決定するため、Gnuplot初心者を悩ませる原因となっている。
一応、表面的には
set pm3d corners2color { c1 | c2 | c3 | c4 }
とすることで解決する（詳細はググればいくらでも見つかるのでそちらを参照されたい）。しかし場合によっては、これで解決しないこともありうる。例えば、pm3dで描画したデータの上に等高線を描きたいときなどだ。
以下の画像を見て欲しい。
これはADAPT-GPM2の試作機能で作った画像だ。ただしpm3dの問題が発覚したため、現在はこの機能は削除されており、上の画像を出力することは出来ない。その代わり問題を修正した代替機能が追加されている。
この画像は、平面上に二つの点電荷が置かれている場合の、それによって作られる各座標での電位と電場の方向を表している。高校、大学の電磁気学あたりを履修していれば見慣れた図だろう。ポテンシャルの値と電場の向きが食い違っているような気がするが見なかったことにしてほしい。
ところで右図の左上のあたり。本来の等高線であればあるはずのない直線部分があるのが分かるだろうか。これが、と言い切ってしまうと語弊があるが、この部分は先のpm3dの弊害を端的に表した部分である。このグラフは上下対称なものになっているはずだが、等高線が少しだけ左下にずれて表示されているのだ。
この問題は上述のcorners2colorと関連している。私はカラーマップを表示する場合、zの値と色とを厳密に結びつけるために、普段はcorners2color c1と指定し、矩形領域のうち左下の点のz値を参照させている。ただこれは、Gnuplotが認識するxy座標が、本来の位置よりもわずかに左下にずれてしまう、という意味である。例えば(x, y, z) = (0, 0, 0)という値を元に-1<=x<=1、-1<=y<=1の範囲の正方形を塗りつぶす時、corners2color c1であればGnuplotに対して(x, y, z)=(-1, -1, 0)と与えなければならないのだ。
この問題は単にカラーマップを表示したいだけならば無視してよい。ただし等高線を描くとなると話は別だ。本来(0, 0)に位置しているはずの点が、ビン幅/2だけ左下にずれてしまうので、等高線全体が左下にずれる。
もう一つ加えて言うなら、corners2colorを与えた状態でも正常に表示させるためのダミーデータによって、本来存在しないはずの左上の直線が出現してしまっている。
また気持ちの悪さもある。等高線がてんでばらばらの色で表示されているのだ。ある特定の一色のみで統一する方法は見つかったが、折角z座標を反映しているのだからカラーバーの色で表示される方法も用意されるべきではないのか。
というわけで、このあたりの問題を吸収する機能をADAPT-GPM2に実装した。これを使えば以下のような図を作ることが出来る。
ソースコードは次の通り。
#include <ADAPT/GPM2/GPMCanvas.h> #include <thread> using namespace adapt::gpm2; double r(double x, double y) { return std::sqrt(x * x + y * y); }; double potential(double x, double y) { double r1 = r(x - 3, y); double r2 = r(x + 3, y); if (r1 == 0.) r1 = 0.000001; if (r2 == 0.) r2 = 0.000001; double p1 = 1 / r1; double p2 = 3 / r2; return p1 - p2; }; double fieldx(double x, double y) { double f1 = (x - 3) / std::pow(r(x - 3, y), 3); double f2 = 3 * (x + 3) / std::pow(r(x + 3, y), 3); return f1 - f2; } double fieldy(double x, double y) { double f1 = y / std::pow(r(x - 3, y), 3); double f2 = 3 * y / std::pow(r(x + 3, y), 3); return f1 - f2; } int main() { //GPMCanvas::SetGnuplotPath("path to gnuplot.exe"); adapt::Matrix<double> m(100, 100); std::pair<double, double> xrange = { -9.9, 9.9 }; std::pair<double, double> yrange = { -9.9, 9.9 }; for (int iy = -50; iy < 50; ++iy) { double y = iy * 0.2 + 0.1; for (int ix = -50; ix < 50; ++ix) { double x = ix * 0.2 + 0.1; m[ix + 50][iy + 50] = potential(x, y); } } std::vector<double> xfrom(441), yfrom(441), xlen(441), ylen(441); std::vector<double> arrowcolor(441); for (int iy = -10; iy <= 10; ++iy) { for (int ix = -10; ix <= 10; ++ix) { int jx = (ix + 10); int jy = (iy + 10); double xlen_ = fieldx(ix, iy); double ylen_ = fieldy(ix, iy); double len = std::sqrt(xlen_ * xlen_ + ylen_ * ylen_); xlen_ = xlen_ / len * 0.8; ylen_ = ylen_ / len * 0.8; xlen[jy * 21 + jx] = xlen_; ylen[jy * 21 + jx] = ylen_; xfrom[jy * 21 + jx] = ix - xlen_ / 2.; yfrom[jy * 21 + jx] = iy - ylen_ / 2.; arrowcolor[jy * 21 + jx] = potential(ix - xlen_ / 2., iy - ylen_ / 2.); } } GPMMultiPlot multi("example_colormap.png", 1, 2, 1200, 600); //左側カラーマップの作成 GPMCanvasCM g1("example_colormap_tmpfile"); g1.ShowCommands(true); g1.SetTitle("example\\_colormap"); g1.SetPaletteDefined({ {0, "yellow" }, { 4.5, "red" }, { 5., "black" }, { 5.5, "blue"}, { 10, "cyan" } }); g1.SetSizeRatio(-1); g1.SetXLabel("x"); g1.SetYLabel("y"); g1.SetXRange(-10, 10); g1.SetYRange(-10, 10); g1.SetCBRange(-5, 5); g1.PlotColormap(m, xrange, yrange, plot::title = "notitle"). PlotVectors(xfrom, yfrom, xlen, ylen, plot::title = "notitle", plot::color = "white"); std::this_thread::sleep_for(std::chrono::milliseconds(300)); //右側等高線図の作成 GPMCanvasCM g2("example_colormap_tmpfile"); g2.ShowCommands(true); g2.SetTitle("example\\_contour"); g2.SetPaletteDefined({ {0, "yellow" }, { 4.5, "red" }, { 5., "black" }, { 5.5, "blue"}, { 10, "cyan" } }); g2.SetSizeRatio(-1); g2.SetXLabel("x"); g2.SetYLabel("y"); g2.SetXRange(-10, 10); g2.SetYRange(-10, 10); g2.SetCBRange(-5, 5); g2.PlotColormap(m, xrange, yrange, plot::title = "notitle", plot::with_contour, plot::without_surface, plot::variable_cntrcolor, plot::cntrlevels_incremental = { -20., 0.2, 20. }). PlotVectors(xfrom, yfrom, xlen, ylen, plot::title = "notitle", plot::variable_color = arrowcolor); return 0; }
こっちはこっちで描画範囲の端で線が僅かに途切れてしまっているが、さすがに我慢するしかないだろう。Gnuplotの阿呆さにこれ以上付き合ってはいられない。
Gnuplotのpm3dは恐らくそもそもメッシュデータの3Dプロットに対して彩色することを意図して作られたものなので、それを考えれば矩形の頂点4つから色を求めるのは妥当ではある。しかし私はpm3dを主として2次元ヒストグラムを描くために使っているため、この妥当なはずの仕様が却って厄介なものになってしまっている。私はhistepsに相当するような機能がほしいのだ。……いや実を言うと、ヒストグラムを描くのにいちいち小細工を要するsteps系の機能も結構面倒なのだが。
3+1世代のステライルニュートリノを仮定した場合の各物理定数について、90% confidence intervalをグラフに描画するために格子点上のデータに対してちょっと等高線を作ってやる必要が出てきて、こんな機能を用意したわけである。たぶん素粒子宇宙系の人間はCERN ROOTを使ってちょろっと描いてしまうのだろうが、私はやむを得ない状況以外ではROOTを使わないことを心に決めているので、代わりにGnuplotと格闘することになった。ROOTと格闘するよりはきっとマシである。
いい加減、pm3dに代わる優れた機能を用意してくれないものだろうか。あれはちょっと厄介事が多すぎると思う。
