格子点で指定する投影法 (変換番号51) の概要と使用法
● 関連パッケージ
GRPH1
SGPACK 一部追加
STPACK 一部追加
MATH1
GT2DLIB 新設： 2次元格子の座標変換
● 概要
DCLは、U 座標から V 座標への投影法として、各種直交(曲線)座標や地図投影
などをサポートしている。本投影法は、その一つとして追加された、離散的な
格子点間の対応で定義する投影法である(一方、他の変換は解析的な関数によ
る)。変換番号は 51 である。格子点の間は双線形補間 (bilinear
interpolation) による補間を行う。このため、新設された MATH1 の GT2DLIB
を用いる。変換の設定はサブルーチン G2SCTR で行う。実際の変換は下位のルー
チンが G2FCTR を呼んで行う(ユーザーは呼ぶ必要なし)。
他の変換と同様、本変換も U 座標において論理的に矩形な2次元の領域を対象
とする。独立変数をξ, ηとするU座標系における格子点 [ξ_i, η_j]
(i=1,..,Nx, j=1,..,Ny) が、基準となる直交座標系(独立変数x,y)において
[x_ij, y_ij] (i=1,..,Nx, j=1,..,Ny) である場合に、前者をユーザーにとっ
ての描画領域とし、後者の座標系で画面上に表示する。ここで、ξ, ηはそれ
ぞれi,jのみに依存するのに対し、x, y は両方に依存することに注意せよ。即
ち、この格子はU座標において「矩形」である。
この変換は、地球流体力学的には、格子生成された格子における流体シミュレー
ションの結果を可視化するときなどに、有用である。また、地形に沿った座標
系での計算結果の可視化にも役立つ。本変換は、特別な場合として、ξ＝x ま
たは η＝y の場合に対応できるので、U座標でξ＝xとして水平方向をとり、
ηに地形に沿った座標、yに高度を取ればよい。
本変換では、U 座標と V 座標を直接対応させるのでなく、C 座標と呼ぶ中間
的な直交座標を介させる。C 座標は Cartesian を語源とし、V 座標に線形に
対応する直交座標である（原点移動と縦・横軸それぞれの伸縮）。C座標は、
曲線座標定義の基準となる直行座標を具現化したものである。基準となる直交
座標系(上記のx,y座標系)は、V座標でなくC座標に対応させ、C 座標からV座標
への変換は、SGPACKの相似変換パラメター CXMIN, CXMAX, CYMIN, CYMAX で指
定する。一方, U 座標から C 座標への変換は、GT2DLIB で行う。
座標軸は C 座標で U[XY]PACK を使って書くことができる。その場合、変換番
号1に切り替えU座標としてC座標を用いればよい。
現在のところ、本変換では線分はV座標において線分となるように引く。この
ため、例えば座標軸に沿った線を一本の線分として書いても、格子の曲がりが
反映されない。従って、ある程度細かな間隔で線分を描きつなげることが必要
である。これは、変換番号2-7でも同様である。一方、地図投影では線分は大
円に沿って適当に区切って補間されることになっており、例えば経線を線分と
して描画しても、地図に沿った曲線(折れ線)として表示される。将来的には、
地図投影以外の変換に関して、U座標における補間を適宜行うようにすること
を検討している。
● 正規変換パラメタ
・相似変換パラメタ
投影変換の基準となる直交座標系と V 座標の間の相似変換を指定するパラ
メターには下記の2系統がある。DCL 5.2 までは2番目のものだけが存在し、
地図投影などに用いられているが、5.3 ではより一般性のある1番目のもの
が導入され、変換番号51などで用いられる。両系統の設定ルーチンは、それ
ぞれ他方も設定するようになっているので、実際にはどちらを設定しても良
いようになっている。
( CXMIN, CXMAX, CYMIN, CYMAX )
投影変換 (変換番号 5 以上. 地図投影や格子点で指定する変換を含む)
において、投影の基準となる直交座標系を V 座標に対応させるためのパ
ラメター。ウィンドーの四隅のC座標を指定する。CXMIN, CXMAX, CYMIN,
CYMAX が VXMIN, VXMAX, VYMIN, VYMAX に一致するとき、C座標とV座標は
一致する。なお、相似変換パラメターは次に述べる SIMFAC, VXOFF,
VYOFF によっても指定できる。こちらは縦軸と横軸の縮尺が同じであるこ
とを補償する。
(CXMIN, CXMAX, CYMIN, CYMAX) はサブルーチンSGSCWD / SGQCWD で設定
/ 参照できる。SGSCWD は、内部で (SIMFAC, VXOFF, VYOFF) を設
定するが、その際、SIMFAC は横軸のスケールを元に決める。
( SIMFAC, VXOFF, VYOFF )
SGSSIMは内部で (CXMIN, CXMAX, CYMIN, CYMAX) を設定する。
● U座標系の設定
・格子点で指定する座標系 (変換番号51)
CALL SGSTRN(51)
CALL SGSVPT(VXMIN, VXMAX, VYMIN, VYMAX)
CALL SGSCWD(CXMIN, CXMAX, CYMIN, CYMAX)
CALL SGSWND(UXMIN, UXMAX, UYMIN, UYMAX)
CALL G2SCTR(NX, NY, UX, UY, CX, CY)
CALL SGSTRF
ここで、UX, UY はそれぞれ長さ NX, NY の一次元配列で、U座標における
矩形に区切られた格子点を表す(概要説明におけるξ_i, η_jに対応)。ま
た, CY, CY はともにサイズ (NX, NY) の2次元配列とする(概要説明おけ
る x_ij, y_ij に対応)。UX, UYは単調増加か単調減少でなければならな
い。G2SCTR は MATH1のGT2DLIBにおける座標変換関数 G2FCTR の初期化関
数である。G2FCTRは描画時に自動的に呼ばれる。なお、UX, UYの格子がぴっ
たり入る U 座標をとるには、GRSWNDで
CALL SGSWND(UX(1),UX(NX),UY(1),UY(NY))
などとすればよい。
なお GRPH2 では上記のサブルーチンの接頭辞SG を GR に変えたサブルー
チンが用意されている。GR* を用いて変換を設定する場合は、(CXMIN,
CXMAX, CYMIN, CYMAX) の指定は省略できる。省略した場合すべての格子
点が入る最小の範囲が選ばれる。例：
CALL GRSTRN(51)
CALL GRSVPT(VXMIN, VXMAX, VYMIN, VYMAX)
CALL GRSWND(UXMIN, UXMAX, UYMIN, UYMAX)
CALL G2SCTR(NX, NY, UX, UY, CX, CY)
CALL GRSTRF
● SGPACKレファレンス追加分
・SGpGET/SGpSET(SGpSTX)
[変換関数に関するパラメター]
'CXMIN', 'CXMAX', 'CYMIN', 'CYMAX' 相似変換パラメタ. 初期値は
0.0, 1.0, 0.0, 1.0.
・SGSCWD
1. 機能
相似変換を設定する。
2. 呼出し方法
CALL SGSCWD(CXMIN,CXMAX,CYMIN,CYMAX)
3. パラメターの説明
CXMIN,CXMAX,CYMIN,CYMAX (R) C座標の四隅を表す SGPACK の
内部変数。1.4節を参照のこと。
4. 備考
(a) この4パラメターは地図投影等で用いられる相似変換パラメター
SIMFAC, VXOFF, VYOFF の異なる表現法と言える。ただし、一つ多い
ことから分かるように、こちらのほうがより一般的である。後者で
は、縦・横の伸縮パラメターに同一値(SIMFAC)を取ることで、パラ
メターの数が減っているのである。
(b) 主に変換番号 51 で用いるために導入されたパラメターであるが、
実際には、SGSCWDを呼ぶと SGSSIM で扱うパラメターも設定され、
また SGSSIM を呼ぶと SGSCWD で扱うパラメターも設定されるよう
になっている。
● GRPACKレファレンス追加分
・GRSCWD
1. 機能
相似変換を設定する。GRPH1/SGPACKのSGSCWDに引数を渡して呼ぶ
だけである。
2. 呼出し方法
CALL GRSCWD(CXMIN,CXMAX,CYMIN,CYMAX)
3. パラメターの説明
GRPH1/SGPACKのSGSCWDを参照のこと.
4. 備考
GRPH1/SGPACKのSGSCWDを参照のこと.
● チュートリアル
DCLでは、直行曲線座標の表示や地図投影など、数多くの投影変換による描画
が行えます。これらはいずれも解析関数より指定される変換ですが、それ以外
に格子点間の対応で変換を指定する投影法もあり、変換番号は51となっていま
す。ここではその使い方を説明します。なお、格子点の間の点の表示には双線
形補間(bilinear interpolation)が用いられます。補間は MATH1 の GT2DLIB
により行われますので、アルゴリズムの詳細については、そのマニュアルを参
照してください。
・例１：
下記のプログラムを実行すると、図 g2pk01.png が現れます。
プログラムg2pk01.f
--------------- プログラム g2pk01.f -------------------
1: *--------------------------------------------------------------------
2: * Copyright (C) 2000-2004 GFD Dennou Club. All rights reserved.
3: *--------------------------------------------------------------------
4: PROGRAM G2PK01
5:
6: PARAMETER(NX=15,NY=15)
7:
8: REAL UX(NX), UY(NY)
9: REAL UYW(NX), UXW(NY)
10: REAL CX(NX,NY), CY(NX,NY)
11: REAL Z(NX,NY)
12:
13:
14: * / SET PARAMETERS /
15:
16: DO 10 I=1,NX
17: UX(I)=(I-1.0)/(NX-1.0)
18: 10 CONTINUE
19: DO 15 J=1,NY
20: UY(J)=(J-1.0)/(NY-1.0)
21: 15 CONTINUE
22:
23: DO 25 J=1,NY
24: DO 20 I=1,NX
25: CX(I,J) = UX(I) + 0.1*UY(J)
26: CY(I,J) = 0.2*UX(I) + UY(J)
27: 20 CONTINUE
28: 25 CONTINUE
29:
30: CXMIN = 0.0
31: CXMAX = 1.1
32: CYMIN = 0.0
33: CYMAX = 1.1
34:
35: * / GRAPHIC /
36:
37: WRITE(*,*) ' WORKSTATION ID (I) ? ;'
38: CALL SGPWSN
39: READ(*,*) IWS
40:
41: CALL GROPN(IWS)
42: CALL GRFRM
43: CALL GRSVPT(0.15,0.85,0.15,0.85)
44: CALL GRSWND(UX(1),UX(NX),UY(1),UY(NY))
45: CALL G2SCTR(NX,NY,UX,UY,CX,CY)
46: CALL GRSTRN(51)
47: CALL SGSCWD(CXMIN,CXMAX,CYMIN,CYMAX)
48: CALL GRSTRF
49:
50: CALL SGLSET('LCLIP',.TRUE.)
51:
52: * / TONE /
53:
54: DO 35 J=1,NY
55: DO 30 I=1,NX
56: Z(I,J) = UX(I) + UY(J)
57: 30 CONTINUE
58: 35 CONTINUE
59:
60: CALL UELSET('LTONE',.TRUE.)
61: CALL UWSGXA(UX,NX)
62: CALL UWSGYA(UY,NY)
63: CALL UETONE(Z, NX, NX, NY)
64:
65: * / GRID LINES /
66:
67: DO 45 J=1,NY
68: DO 40 I=1,NX
69: UYW(I) = UY(J)
70: 40 CONTINUE
71: CALL SGPLU(NX,UX,UYW)
72: 45 CONTINUE
73:
74: DO 55 I=1,NX
75: DO 50 J=1,NY
76: UXW(J) = UX(I)
77: 50 CONTINUE
78: CALL SGPLU(NY,UXW,UY)
79: 55 CONTINUE
80:
81: * / AXES (Switch to ITR==1) /
82:
83: CALL GRSWND(CXMIN,CXMAX,CYMIN,CYMAX)
84: CALL GRSTRN(1)
85: CALL GRSTRF
86: CALL USDAXS
87: CALL UXSTTL('T','BOTH X & Y TRANSFORMED',0.0)
88:
89: CALL GRCLS
90:
91: END
--------------- プログラム g2pk01.f -------------------
デモのため、格子点はごく簡単な解析関数を用いて定義しました。 U 座標に
おいて「長方形」に並んだ格子点の座標が [ UX(I), UY(J) ] (I=1..NX,
J=1..NY) で定義され、それぞれが基準となる直交座標において [ CX(I,J),
CY(I,J) ] の座標値を持ちます。
42-48行目で、変換を設定しています。U座標における描画範囲は, UX, UY の
カバーする領域ぴったりにすべく、
CALL GRSWND(UX(1),UX(NX),UY(1),UY(NY))
と設定しました(44行目)。一方、ビューポートの四隅のC座標値を指定する
CXMIN,CXMAX,CYMIN,CYMAX には適当な値を設定しています(30-33 行目)。
GRFRM, GRSTR等、GRPH2のGRPACKを使って初期化する場合、C[XY](MIN|MAX)の
設定は省略できます(GRFRMはCXMIN等を不定にします)。省略時には、格子をぴっ
たり収めるべく、CX, CYの最大値を用いて内部的に設定します。
U座標で色塗り(54-63行目)をしたあと、格子を表示します(67-79行目)。現在
のところ、線分はV座標でそのまま線分として書くようになってますので、こ
のように各格子点を結ぶ折れ線として書く必要があります。いずれにしろ、双
線形補間においていは、このようにすると厳密な格子の表示となります。
座標軸は C 座標で書きました(83-87行目)。より正確に言えば、変換を定義し
なおして、それまでC座標であった座標系をU座標として定義した上で(変換番
号は1)、座標軸を書きました。今のところ、変換番号51における直接の差表示
区描画はサポートされていません。
・例２：
下記のプログラムを実行すると、図 g2pk02.png が現れます。
プログラムg2pk02.f
--------------- プログラム g2pk02.f -------------------
1: *--------------------------------------------------------------------
2: * Copyright (C) 2000-2004 GFD Dennou Club. All rights reserved.
3: *--------------------------------------------------------------------
4: PROGRAM G2PK02
5:
6: PARAMETER(NX=15,NY=15)
7:
8: REAL UX(NX), UY(NY)
9: REAL UYW(NX), UXW(NY)
10: REAL CX(NX,NY), CY(NX,NY)
11: REAL Z(NX,NY)
12: REAL TERRAIN(NX)
13:
14:
15: * / SET PARAMETERS /
16:
17: CALL GLRGET('RUNDEF',RUNDEF)
18:
19: DO 10 I=1,NX
20: UX(I)=(I-1.0)/(NX-1.0) - 0.5
21: TERRAIN(I) = 0.1 * EXP(-24*UX(I)**2)
22: 10 CONTINUE
23: DO 15 J=1,NY
24: UY(J)=(J-1.0)/(NY-1.0)
25: 15 CONTINUE
26:
27: CX(1,1) = RUNDEF
28: DO 25 J=1,NY
29: DO 20 I=1,NX
30: CY(I,J) = UY(J)*(1.0-TERRAIN(I)) + TERRAIN(I)
31: 20 CONTINUE
32: 25 CONTINUE
33:
34: * / GRAPHIC /
35:
36: WRITE(*,*) ' WORKSTATION ID (I) ? ;'
37: CALL SGPWSN
38: READ(*,*) IWS
39:
40: CALL GROPN(IWS)
41: CALL GRFRM
42: CALL GRSVPT(0.15,0.85,0.15,0.85)
43: CALL GRSWND(UX(1),UX(NX),UY(1),UY(NY))
44: CALL GRSTRN(51)
45: CALL G2SCTR(NX, NY, UX,UY, CX,CY)
46: CALL GRSTRF
47:
48: * / TONE /
49:
50: DO 35 J=1,NY
51: DO 30 I=1,NX
52: Z(I,J) = UX(I) * (1-UY(J))
53: 30 CONTINUE
54: 35 CONTINUE
55:
56: CALL UELSET('LTONE',.TRUE.)
57: CALL UWSGXA(UX,NX)
58: CALL UWSGYA(UY,NY)
59: CALL UETONE(Z, NX, NX, NY)
60:
61: * / GRID LINES /
62:
63: DO 45 J=1,NY
64: DO 40 I=1,NX
65: UYW(I) = UY(J)
66: 40 CONTINUE
67: CALL SGPLU(NX,UX,UYW)
68: 45 CONTINUE
69:
70: DO 55 I=1,NX
71: DO 50 J=1,NY
72: UXW(J) = UX(I)
73: 50 CONTINUE
74: CALL SGPLU(NY,UXW,UY)
75: 55 CONTINUE
76:
77: * / AXES (Switch to ITR==1) /
78:
79: CALL G2QCTM(CXMIN, CXMAX, CYMIN, CYMAX)
80: CALL GRSWND(CXMIN,CXMAX,CYMIN,CYMAX)
81: CALL GRSTRN(1)
82: CALL GRSTRF
83: CALL USDAXS
84: CALL UXSTTL('T','TERRAIN FOLLOWING',0.0)
85:
86: CALL GRCLS
87:
88: END
--------------- プログラム g2pk02.f -------------------
この例では, CX＝UXとし、Y軸に沿ってのみ座標変換しています。UXは1次元配
列なのに対しCXは2次元ですから、そのような配列を別途用意しなければなら
ないとなると少し面倒です。そこで、CXの最初の要素の値が GLPACK のパラメ
ター 'RUNDEF' の値に等しい場合は、CX＝UXであると解釈されることになって
います。それ以降の値は読み込みませんので、CXのために2次元配列を用意す
る必要もありません（例では用意していますが）。同様なことは Y 軸に沿っ
ても出来ます。ただし、CX＝UXかつCY＝UYならこの変換の意味がないので、エ
ラーになるようになっています。
座標変換の設定(41-46行目)において例１と違うのは、C[XY](MIN|MAX) を設定
していないことです。こうすると、内部的に格子点がぴったり収まるような領
域が確保されますので、図のような結果となります。ただしそのためには 例
のように GRPH2 の設定ルーチン GRFRM 等を用いる必要があります。なお、
内部的に設定されたこれらのパラメターを参照することで、C座標への切り替
えと軸描画もできます（79-84行目)。