2次元のスカラー場を手早くコンタリングして見たいというときには, 等高線 描画パッケージ UDPACK にあるメソッド(元サブルーチン) udcntr を1つ呼ぶだけで 十分です. 例題として, 球面調和メソッド(元関数)の重ね合わせ の実部を描いてみましょう(quick3).
# quick3.rb require "narray" require "numru/dcl" include NumRu include Math nx = 37 ny = 37 xmin = 0 xmax = 360 ymin = -90 ymax = 90 drad = PI/180 p = NArray.sfloat(nx, ny) #-- data ---- for j in 0..ny-1 for i in 0..nx-1 alon = (xmin + (xmax-xmin)*i/(nx-1))*drad alat = (ymin + (ymax-ymin)*j/(ny-1))*drad slat = sin(alat) p[i,j] = 3*sqrt(1-slat**2)*slat*cos(alon) - 0.5*(3*slat**2-1) end end #-- graph ---- iws = (ARGV[0] || (puts ' WORKSTATION ID (I) ? ;'; DCL::sgpwsn; gets)).to_i DCL::gropn iws DCL::grfrm DCL::grswnd(xmin, xmax, ymin, ymax) DCL::grsvpt(0.2, 0.8, 0.2, 0.8) DCL::grstrn(1) DCL::grstrf DCL::usdaxs DCL::udcntr(p) DCL::grclsprogram quick3
udcntr は等高線図を描くだけで, 枠や座標軸を描いてくれませんから,
この例ではまず, いくつかの「おまじない」(grswnd, grsvpt,
grstrn, grstrf)で座標系を決め, USPACKのルーチン usdaxs を使っておまかせの座標軸を描画したのちに, udcntr を呼ん
で等高線図を描いています.
udcntr の引数は, 最初のp が作画しようとする場を与える2次元
配列, 2番目のnx が実際に宣言した配列の第1次元寸法, 3, 4番目の
nx, ny がそれぞれ, 作画に使う配列の第1次元および第2次元寸
法です. この例のように, 宣言した配列をフルに用いて作画するときは, 2番
目と3番目の引数は同じになりますが, 第1次元目の寸法をこのように2つの引
数で指定するといいことがあります. その利点は第9.3節で説
明します.
この例のように udcntr だけを呼ぶ場合には, 等間隔の格子点を設定し
てコンタリングを行ないます. つまり, 座標軸の目盛に関係なく, P(1,1) が左下隅, P(NX,NY) が右上隅にくるような等間隔の格子点座
標が設定されます. また, コンターレベルも自動的に決定されて, 図の下にそ
のコンター間隔が表示されます.
GRPH1 の ユーザーインターフェイス SGPACK には, 折れ線を描いたり文字列
を描いたりする機能の他にも, 多角形閉領域のぬりつぶしを行なう機能があり
ます. UEPACK はこれに対応するトーンぬりつぶしの上位ルーチンですが,
quick3 の例で uetone メソッド(元サブルーチン)を1つ呼ぶだけで, とりあえ
ず負の領域に斜線のハッチをつけることができます(quick4). ここで
uetone を呼ぶのが31行めで, 座標軸やコンターを描くよりも前である
ことに注意しておきましょう. これらの順序を変えると, 出力装置によっては
結果が異なります. それまでに描いた部分がぬりつぶされ消されてしまうこと
があるのです. この点については, 第3.5節, 第
10.1節で詳しく説明します.
UDPACK や UEPACK, UWPACK の諸機能を使うことによって, 高度な等高線図が 描けるようになります. カラーグラフィクスが利用できる環境では, このよう な等高線図を色の塗り分け(カラー グラデーション)で表現することも可能に なります. これらは, 第9, 10, 12章で説 明することにしましょう.
# quick4.rb require "narray" require "numru/dcl" include NumRu include Math nx = 37 ny = 37 xmin = 0 xmax = 360 ymin = -90 ymax = 90 drad = PI/180 p = NArray.sfloat(nx, ny) #-- data ---- for j in 0..ny-1 for i in 0..nx-1 alon = (xmin + (xmax-xmin)*i/(nx-1))*drad alat = (ymin + (ymax-ymin)*j/(ny-1))*drad slat = sin(alat) p[i,j] = 3*sqrt(1-slat**2)*slat*cos(alon) - 0.5*(3*slat**2-1) end end #-- graph ---- iws = (ARGV[0] || (puts ' WORKSTATION ID (I) ? ;'; DCL::sgpwsn; gets)).to_i DCL::gropn iws DCL::grfrm DCL::grswnd(xmin, xmax, ymin, ymax) DCL::grsvpt(0.2, 0.8, 0.2, 0.8) DCL::grstrn(1) DCL::grstrf DCL::uetone(p) DCL::usdaxs DCL::udcntr(p) DCL::grclsprogram quick4