gnuplot でメルカトルな平面に円を描きたい

# 目的
球面を平面に投影する方法としてもっともポピュラーなものはメルカトル図法ではないかと思います．緯度方向と経度方向が直行するので図示するのはとても簡単です．一方で，極に近い部分では形状が大きく歪むという点もメルカトル図法の特徴です．

メルカトル図法で表現された平面上で，ある点を中心とした角半径 &theta; の円を描くことを考えます．このとき，緯度経度平面上で角半径 &theta; の円を描いてはいけません．半径 &theta; とはあくまで球面上で定義されるものだからです．円を構成する各点の座標を何とかして求める必要があります．

<ins title="ふええ(´；ω；｀)">今更気づきましたがこれメルカトル図法じゃないですね……正しくは Equirectangular projection です． Mercator projection との違いについては末尾の Wikipedia へのリンクを参照してください．</ins>

# 解決方法
方針は次のとおりです．まずは，円を最も描きやすい場所で考えます．緯度 &alpha; - 経度 &delta; ととった系で円を描くには極が最も簡単です．まずは極を中心として半径 &theta; の円を描きます．その後，円の中心が指定された位置に来るように座標系を回転させます．座標系の回転には回転行列を用います．最終的に次のようなスクリプトを作成しました．

```gnuplot:plot.gp
#!/usr/local/bin/gnuplot
#
set mxtics 5
set mytics 4
set grid front
set grid x y mx my
set grid lt 0 lw 1 lc rgb 'gray40', lt 0 lc rgb 'gray60'
set key samplen 2
set xtics format "%.0f" offset 0,0.2
set xlabel "Lon (degree)" offset 0,0.1
set xrange [ 0.0000 : 360.0000 ] reverse
set ytics format "%+4.0f" offset 0.2,0
set ylabel "Lat (degree)" offset 0.8,0
set yrange [ -90.0000 : 90.00000 ] noreverse

rad2deg(r) = 180.*(r/(4.*atan(1.)));
deg2rad(d) = (d/180.)*4.*atan(1.);
ad2x(a,d) = cos(a)*cos(d)
ad2y(a,d) = sin(a)*cos(d)
ad2z(a,d) = sin(d)
xyz2a(x,y,z) = atan2(x,y)+pi
xyz2d(x,y,z) = asin(z)
euler1a(p,a,d) = \
xyz2a(cos(p)*ad2x(a,d)-sin(p)*ad2y(a,d),\
      sin(p)*ad2x(a,d)+cos(p)*ad2y(a,d),\
      ad2z(a,d))
euler1d(p,a,d) = \
xyz2d(cos(p)*ad2x(a,d)-sin(p)*ad2y(a,d),\
      sin(p)*ad2x(a,d)+cos(p)*ad2y(a,d),\
      ad2z(a,d))
euler2a(p,a,d) = \
xyz2a(ad2x(a,d),\
      cos(p)*ad2y(a,d)-sin(p)*ad2z(a,d),\
      sin(p)*ad2y(a,d)+cos(p)*ad2z(a,d))
euler2d(p,a,d) = \
xyz2d(ad2x(a,d),\
      cos(p)*ad2y(a,d)-sin(p)*ad2z(a,d),\
      sin(p)*ad2y(a,d)+cos(p)*ad2z(a,d))
euler3a(p,a,d) = \
xyz2a(cos(p)*ad2x(a,d)-sin(p)*ad2y(a,d),\
      sin(p)*ad2x(a,d)+cos(p)*ad2y(a,d),\
      ad2z(a,d))
euler3d(p,a,d) = \
xyz2d(cos(p)*ad2x(a,d)-sin(p)*ad2y(a,d),\
      sin(p)*ad2x(a,d)+cos(p)*ad2y(a,d),\
      ad2z(a,d))
###
alpha1(x,y) = euler1a(atan2(x,y),0,pi/2-sqrt(x**2+y**2))
delta1(x,y) = euler1d(atan2(x,y),0,pi/2-sqrt(x**2+y**2))
alpha2(x,y,d) = euler2a(d-pi/2.0,alpha1(x,y),delta1(x,y))
delta2(x,y,d) = euler2d(d-pi/2.0,alpha1(x,y),delta1(x,y))
alpha3(x,y,a,d) = euler3a(-a+pi/2.0,alpha2(x,y,d),delta2(x,y,d))
delta3(x,y,a,d) = euler3d(-a+pi/2.0,alpha2(x,y,d),delta2(x,y,d))
alpha(x,y,a,d) = rad2deg(alpha3(deg2rad(x),deg2rad(y),deg2rad(a),deg2rad(d)))
delta(x,y,a,d) = rad2deg(delta3(deg2rad(x),deg2rad(y),deg2rad(a),deg2rad(d)))

set parametric
set tr [0:2*pi]

set sample 600
x(t) = 15*sin(t)
y(t) = 15*cos(t)

plot \
     for [n=1:9] \
     alpha(x(t),y(t),30+30*n,-100+20.*n),\
     delta(x(t),y(t),30+30*n,-100+20.*n) \
     w p notitle ps 0.2 lc rgb hsv(160+5*n,196,165+10*n)
#    EOF
```

このスクリプトでは指定の位置に半径 15 度の円を描きます．中心の緯度は -80 から +80 度まで変化します．実際に wxt terminal で作図した結果を以下に表示します．極に近い部分では円の形状が大きくゆがんでいる様子が確認できます．

![mercator_c.png](https://qiita-image-store.s3.amazonaws.com/0/46396/fdef7c82-c855-d390-1773-b69c5b36f4b1.png "半径 15 度の円をメルカトル図法で")

`rad2deg`, `deg2rad` は度数法と弧度方の変換に用いた関数です．`ad2x`, `ad2y`, `ad2z` という 3 つの関数は経度-緯度の形式 (&alpha;,&delta;) から (x,y,z) の表示に変換するためのものです．`xyz2a`,`xyz2d` は逆に (x,y,z) の形式から (&alpha;,&delta;) の表示に変換するための関数です．`euler1a`, `euler1d` という関数はオイラー角の形式で (&alpha;,&delta;) のペアを回転させる関数です．オイラー角に習って 3 ペアの関数がありますが `euler1a`, `euler1d` と `euler3a`, `euler3d` は本質的に同じ演算をしています．オイラー角の詳細については Wikipedia の [オイラー角][euler] の項を参照してください．また，色に指定に使用している `hsv` 関数については[僕の過去の記事][color.gp]を参照してください．

今回用いた手法は円以外の図形にも用いることができます．関数 `x`, `y` を以下のように変更すると，下記のようなプロットになります．

```gnuplot
x(t) = 10*sin(3*t)
y(t) = 10*cos(5*t)
```

![mercator_r.png](https://qiita-image-store.s3.amazonaws.com/0/46396/40dd2754-a6f8-d2b1-0235-45a9ef29a899.png "リサジュー図形をメルカトル図法で")

今回はなんとか目的のプロットを得ることができました．しかし，かなり力技で何とかした感じが否めないので，もしもっとエレガントな方法を知っている方がいらっしゃいましたらコメントしていただけるとうれしいです．また，今回作成したスクリプトはあまり詳細なテストをしていないのでバグがあるかもしれません．試してみて変な挙動を見せたら報告してくださるとありがたいです．


# 参考資料
- [Wikipedia - オイラー角][euler]
- [gnuplot でプロットに好きな色を使う][color.gp]
- [Wikipedia - Equirectangular projection][wiki-eqtan]
- [Wikipedia - Mercator projection][wiki-merca]

[euler]: http://ja.wikipedia.org/wiki/%E3%82%AA%E3%82%A4%E3%83%A9%E3%83%BC%E8%A7%92
[color.gp]: http://qiita.com/xr0038/items/80e30db5fa7bcd391b18
[wiki-eqtan]: http://en.wikipedia.org/wiki/Equirectangular_projection
[wiki-merca]: http://en.wikipedia.org/wiki/Mercator_projection