こちらの記事が面白かったので
http://d.hatena.ne.jp/rikunora/20100905/p1
昔作った楕円面上の測地線を計算するコードを貼っておきます。コード中の計算式はA4一枚くらい使って手計算して結果だけ書いてあるんで今読むとどう導出したかよく分かりません^^;
筋道としては緯度と経度を時間の関数にしてある時間の間の経路長を積分で表してそれが極小になるようオイラー・ラグランジュ方程式を緯度と経度について出した後、出てきた二つの式を使って緯度を経度の関数と考えた場合の一つの微分方程式に直して解いています。
追記
上の計算はR^nからR^m (楕円面だと緯度経度の二次元面から三次元ユークリッド空間)へ解析関数でマップされる多様体の場合ですが、もっと一般的に三角メッシュで現される曲面について測地線を計算しました。データはフリーのものを拾って改造。これなら手計算不要。CPUパワーに任せて細かく面を分割すれば、どんな面でもOK. まあメッシュ生成自体、結構難しくて奥が深いんですが。アルゴリズムは単純に、三角形上で直線で進み、辺に来たら二つの三角を伸ばして平にして直線を進む。
データはこちらからダウンロードhttp://artist-3d.com/free_3d_models/dnm/model_disp.php?uid=592&count=count
以下楕円面の測地線コード
// Calc geodecics on an ellipsoid // z len = 1 // x and y len = r0 #define r0 (4.0) #include <stdio.h> #include <stdlib.h> #include <math.h> int main(int argc, char **argv) { double th=0.0; //initial longitude double f = 0.0; // initial latitude double dt = 0.00001; double vf = (double)atof(argv[1]); // slope of initial velocity= df/dth double R = r0*r0-1.0; int i=0; while(th<M_PI*2) { double s=sin(f); double c=cos(f); double d0 = (1+R)*c*c; double c0 = 1+R*s*s; double d1 = -2*c*s*(1+R); double c1 = R*2*c*s; double af = (d1-c1*vf*vf)/(2*c0) + vf*vf*d1/d0; // (D'-C'f'f')/2C + D'/D f'f' // V2= C(phi)vphi^2+D(phi)vth^2 f += vf*dt; vf += af*dt; th += dt; if(i%3000==0 || th >= M_PI*2) { double x = cos(th)*r0 * cos(f); double y = sin(th)*r0 * cos(f); double z = sin(f); printf("%f %f %f %f %f\n",th, f, x,y,z); } i++; } printf("\n\n\n"); }
