2010年05月10日

子午線孤長を緯度から求める

http://www.gsi.go.jp/common/000051178.pdf
にあるソースをAWKからC++に移植しました。jtは3で固定しています。

[訂正]5/12

タイポ修正。

pre環境でもinclude内がタグと認識されてしまうので、関数だけ載せることにしました。



//緯度(degree)から子午線孤長(km)を求める。
//



double Meridian_arc_length(double deg)
{

double a=6378.137 ;//地球半径[km]
double rf=298.257222101 ;//?
double n=0.5/(rf-0.5) ;//_?
double n15; n15=1.5*n ;
double an; an=a/(1.0+n);
double rho=3.14159265358979323846264/180.0 ;


double s[7]={0.0};//sin2l\phi 
int jt=3;
int jt2=6;
double ep=1.0 ;
double e[7];
double phi;
double phi2;
double c2;


for(int k=1; k<=jt; k++)
{
e[k]=n15/k-n ;
e[k+jt]=n15/(k+jt)-n;
ep*=e[k];
}




phi=deg*rho ;//radian

phi2=2.0*phi ;double dc;dc=2.0*cos(phi2) ;

double t[7];
t[1]=sin(phi2);s[1]=t[1];

for(int i=2; i<=jt2; i++) { s[i]=dc*s[i-1]-s[i-2] ; t[i]=s[i]/(double)i; } //

double sum=0.0 ;double cs1=ep;
//ここまで準備
for(int j=jt; j>=1; j--)
{
c2=1.0 ;double cs2=phi;
for(int l=1; l<=j; l++)
{
c2/=e[j-l+1] ;
cs2+=c2*t[2*l-1];
c2*=e[j+l] ;
cs2+=c2*t[2*l];
}

sum+=cs1*cs1*cs2 ;
cs1/=e[j];
}

double S=an*(sum-2.0*n*s[1]/sqrt(n*(n+dc)+1.0)+phi) ;
cout.precision(10);

return S;




}



ラベル:C++
posted by hougi at 16:22| Comment(1) | 仕事 | このブログの読者になる | 更新情報をチェックする
×

この広告は180日以上新しい記事の投稿がないブログに表示されております。