平面における回転操作には、複素数が有効である。
z'=z*(cosθ+isinθ) において、z'=x'+y'i 、z=x+yi と置くと、
x'+y'i=(x+yi)(cosθ+isinθ)=xcosθ-ysinθ+(xsinθ+ycosθ)i から、
x'=xcosθ-ysinθ 、y'=xsinθ+ycosθ
よって、
[x'] [cosθ -sinθ][x]
[y']=[sinθ cosθ][y]
と回転行列が現れることからもわかる。
では、空間における回転操作はどうでしょう？これにはあまり学習する機会がないが聞い
たことはある四元数(quaternion)が有効であると知った。
四元数 ｑ は、 ｑ=a+b*i+c*j+d*k （a、b、c、d∈R） で表される数である。
（ただし、i2=j2=k2=i*j*k=-1）
これから、i*j=k、j*k=i、k*i=j、i*j=-j*i、j*k=-k*j、k*i=-i*k 等の計算規則が派生する。
従って、積に関しては非可換となり、２つの四元数ｐ、ｑの積は、一般に、積の順番を入れ
替えると異なるものになる。
これだけの法則から、３次元での一般に回転軸が単位ベクトル(ax，ay，az)の回りを点
（a，b，c)がθ回転（回転軸方向へ右ねじが進む動きがプラスとする。）したときの座標は、
２つの四元数ｐ、ｑを、 四元数p=a*i+b*j+c*k
q=cos(θ/2)+ax*sin(θ/2)*i+ay*sin(θ/2)*j+az*sin(θ/2)*k
とおくことで、q*p*conjugate(q) を計算すればよい。
（共役とは虚数部分の符号を全部qと反対にする。）
(計算結果は実数部は必ず0となり、i、j、k の係数にあたるものがそれぞれ回転後のx、y、z
座標を表す。）
そこで、一般に、ベクトル(l，m，n)を回転軸方向としたとき、
（ax=l/(l2+m2+n2)、ay=m/(l2+m2+n2)、az=n/(l2+m2+n2)が対応）
点(a，b，c)がθ°回転した後の座標をひたすら計算してみました。
移動後のｘ座標は、
1/(2(l2+m2+n2)2)*(al2+al4+2blm-am2+2al2m2+am4+2cnl-an2+2al2n2+2am2n2+an4+(-2l(bm+cn)
+a(l4+m2+m4+n2+2m2n2+n2+l2(-1+2m2+2n2)))*cos(πθ/180)+2(cm-bn)(l2+m2+n2)*sin(πθ/180))
移動後のy座標は
1/(2(l2+m2+n2)2)*(bm2+bl4+2cmn-bn2+2bl2m2+bm4+2alm-bl2+2bl2n2+2bm2n2+bn4+
(-2m(al+cn)+b(l4-m2+m4+n2+2m2n2+n24+l2(1+2m2+2n2)))*cos(πθ/180)+
2(an-cl)(l2+m2+n2)*sin(πθ/180))
移動後のz座標は
1/(2(l2+m2+n2)2)*(cn2+cl4+2anl-cl2+2cl2m2+cm4+2bmn-cm2+2cl2n2+2cm2n2+cn4+
(-2n(al+bm)+c(l4+m2+m4-n2+2m2n2+n4+l2(1+2m2+2n2)))*cos(πθ/180)+
2(bl-am)(l2+m2+n2)*sin(πθ/180))
平面の回転に比べ格段に複雑きわまりない式と相成りましたが（タイプしていてサイクリックな部
分はあるのですが微妙にずれる部分がまちがっていないか心配です。）、これを四元数の道具を使わず
に導き出す事は想像できません。
イヤー四元数は高校や中学でも教えるべきですよ。実数での整数や複素数でのガウス整
数にいろいろ面白い性質があるのと同様、四元数での整数はどう定義するか？、素数に相
当するものは何か？、ラグランジュの平方和の定理（整数は必ず高々４つの平方和で表せる。）など
への応用等ハミルトンが思わず石橋に四元数の規則を刻みつけた気持ちが何となく理解で
きます。
複素数で四元数も再定義できるのでしょうが、だからといって現行のようにこの数にまった
く触れない教育制度は改めるべきです。
上の式は一般的にしたから複雑なようにみえますが、x、y、z軸の回りでの回転などならちょ
ちょいのちょいで計算できます。
アニメーション作成のプログラムや宇宙船の姿勢制御などには頻繁に四元数での回転計
算で処理されているらしいです。現代においては、この数は必須のアイテムでしょう。
ＧＡＩさんからの続報です。（平成２７年１月１４日付け）
３次元での回転計算を使って、回転軸を地球の中心と日本の東京（東経140度、北緯36度）
を結ぶ直線とし、地球をその回転軸の回りに60度回転させたとき北極がどの位置になるか計
算してみました。
経度と緯度から、３次元での(x，y，z)座標に変換すると、地球の半径を約6400ｋｍとして計
算すると、(-3957，3310，3738)なるので、これを回転軸ベクトルにとり、点(0，0，6400)の60
度回転による行き先を計算してみると、(0.450436，0.53854，4783.5)になりました。
この点を逆に緯度、経度へ変換すると、
緯 度: -6°25′56.11717″
経 度: +173°10′57.56817″
となり、これを地図で確認しましたら、パプアニューギニアの東約2000ｋｍの太平洋の海上の
場所となりました。まさに、ここは赤道直下であり、北極の氷も白熊も一切なくなります。地図
ではメルカトル的視野に馴れているせいか、緯度が高い地域の遠近感がずれてしまうので
北極が日本からどれ位離れているか実感が涌きませんが、こうして回転させてみると少し感
覚的に理解できます。
これから先も地球の自転軸がずれて、東京などに近づかないことを願うばかりです。
移動先の(x，y，z)座標から緯度・経度へ変換するとき、入力単位をmで入力するところを、
kmで入力していたため改めてその部分をやり直したら（また入力時y座標の値0.53854を
0.053854と入力していた可能性もあります。）
緯 度: +89°59′29.99495″
経 度: +50°05′26.94335″
と大きく食い違っていました。
緯度と経度から移動先の地球上の場所を探すと緯度はほとんど変わらない結果になりま
した。なお座標から緯度と経度への変換計算は「緯度・経度と地心直交座標の相互換算」
のサイトを利用しています。
こんどは逆に、あまりにも変化がなさ過ぎるような気もしてきました？
これは地球を楕円体と見なした正確に数値へ変換する公式を基にする計算なので、私が
地球を球形として求めた回転後の(x,y,z)の数値が、この入力には不適切と考えないといけ
ないのかも・・・
らすかるさんからのコメントです。（平成２７年１月１４日付け）
そうですね。そんなに変わらないはずはありません。普通に考えて、北緯36°から北極ま
では6000kmぐらいであり、60°回転するということは（左回転として）東京から能登半島ぐら
いの方向に6000km行ったところですから、「Googleマップで大圏航路を表示する」のサイトで
手作業で適当に合わせると、タシケント（ウズベキスタン）あたり（北緯41°ぐらい）にならな
いとおかしいと思います。
DD++さんからのコメントです。（平成２７年１月１５日付け）
北極を点P、東京を点T、求める点をA、地球半径を1とすると、球面三角形 PAT は、
PT = AT = 54°, ∠T = 60° の二等辺三角形で、余弦定理より、
cosPA = cosPT cosAT + sinPT sinAT cos∠T = (13-√5)/16
なので PA = 47.72°、つまり点Aは、北緯 42.28° くらいの地点です。
ついでに、再び余弦定理より、
cos∠P = (cosAT - cosPA cosPT) / (sinPA sinPT) = √((35-6√5)/209)
なので ∠P = 71.25°、つまり点Aは東経 68.75° または西経 148.75° ですね。
ＧＡＩさんからのコメントです。（平成２７年１月１５日付け）
ヘーこうやって求められるんですね。初めて球面三角法の使用法の有効性がわかりま
した。上記の計算を追随したんですが、結果の√((35-6√5)/209)に整理するまでが大
変ですね。
sinPA?と思ったんですが、通常に、sin2PA+cos2PA=1 でいいんですよね。
球面三角形PTAの内角の和が60+2*71.25=202.5(度)であるのは面白いです。
PTの6000kmに較べ60度回転にもかかわらずPAの大圏航路は5330km程度であることも
面白い。それだけ地球は曲がっているんですね。球面三角法も高校で教えましょうよ！
らすかるさんからのコメントです。（平成２７年１月１５日付け）
私が求めた「タシケント（ウズベキスタン）」から約100km、ほぼ合っていて、考え方に間
違いがなかったことが確認できました。
りらひいさんからのコメントです。（平成２７年１月１５日付け）
「Googleマップで大圏航路を表示する」のサイト、いいですね。起点（または終点）マーカー
をクリックすると、緯度経度とともに距離と方角も出てくるので６０度または３００度の方角も
きちんとあわせられますね。
わたしは東京中心の正距方位図法の地図を使えばいいんじゃないかと思いましたけど、
大体の位置の目安を見付けることはできても、詳細な地図でないと正確な場所まではわ
かりませんね。目星を付けるだけなら地球儀でいいかもしれませんが。
らすかるさんからのコメントです。（平成２７年１月１５日付け）
その機能を知らず、60°を目分量で見てました・・・
ＧＡＩさんからのコメントです。（平成２７年１月１７日付け）
一般に、ベクトル(l，m，n)を回転軸方向としたとき、（ax=l/(l2+m2+n2)，ay=m/(l2+m2+n2)，
az=n/(l2+m2+n2)が対応）点(a，b，c)がθ°回転した後の座標をひたすら計算してみました。
移動後のｘ座標は、
FX(a，b，c，θ，l，m，n)
=1/(2(l2+m2+n2)2)*(al2+al4+2blm-am2+2al2m2+am4+2cnl-an2+2al2n2+2am2n2+an4+(-2l(bm+cn)
+a(l4+m2+m4+n2+2m2n2+n2+l2(-1+2m2+2n2)))*cos(πθ/180)+2(cm-bn)(l2+m2+n2)*sin(πθ/180))
移動後のｙ座標は、
FＹ(a，b，c，θ，l，m，n)
=1/(2(l2+m2+n2)2)*(bm2+bl4+2cmn-bn2+2bl2m2+bm4+2alm-bl2+2bl2n2+2bm2n2+bn4+
(-2m(al+cn)+b(l4-m2+m4+n2+2m2n2+n24+l2(1+2m2+2n2)))*cos(πθ/180)+
2(an-cl)(l2+m2+n2)*sin(πθ/180))
移動後のｚ座標は、
FＺ(a，b，c，θ，l，m，n)
=1/(2(l2+m2+n2)2)*(cn2+cl4+2anl-cl2+2cl2m2+cm4+2bmn-cm2+2cl2n2+2cm2n2+cn4+
(-2n(al+bm)+c(l4+m2+m4-n2+2m2n2+n4+l2(1+2m2+2n2)))*cos(πθ/180)+
2(bl-am)(l2+m2+n2)*sin(πθ/180))
を折角求めていたので、回転軸ベクトル(l，m，n)を単位ベクトル
l = Cos[36 Degree]*Cos[140 Degree]、m = Cos[36 Degree]*Sin[140 Degree]、n = Sin[36 Degree]
で設定し（地球中心の東京方向のベクトル）上記の関数式で北極点(a,b,c)=(0，0，1)、θ=60
を代入して計算しましたら、
FX(0，0，1，60，l，m，n)=0.268218、FY(0，0，1，60，l，m，n)=0.689545、FZ(0，0，1，60，l，m，n)=0.672746
これから北極点の移動先での地心直交座標(X，Y，Z)は、これを6400*103(m)倍することで、
(X，Y，Z)=(1716595，4413088，4305574)
がとれました。これを例のサイトの緯度、経度へ変換するものへ代入しましたら、
緯 度 +42°28′12.06133″ 、経 度 +68°44′42.18002″
となり、DD++さんが球面三角法を使って出された結果と殆ど一致しました。
これで上記の関係式も正常に機能出来ることが確認でき苦労の甲斐がありました。
（前回の投稿も大体この趣旨で計算させた積もりでしたが、どこでどう間違ったのかプログラ
ムを消したのでわからなくなりました。）
