/*************************************************************************
２地点間の距離を求める
入力：２地点の緯度、経度
出力：直距離、地表面距離、半射程角
19999.08.08
By T.Kubota
*************************************************************************/
#include <stdio.h>
#include <math.h>
#define PI 3.141592653589793238462 /* 円周率 */
int main(void) {
double f1,f2; /* ２地点の緯度(°) */
double fr1,fr2; /* ２地点の緯度(rad) */
double g1,g2; /* ２地点の経度（東経基準）(°) */
double gr1,gr2; /* ２地点の経度（東経基準）(rad) */
double h1,h2; /* 標高(m) */
double a=6378136.0; /* 赤道半径(m) */
double e2=0.006694470; /* 地球の離心率の自乗 */
double x1,y1,z1,x2,y2,z2; /* ２地点の直交座標値(m) */
double r; /* ２地点間の直距離(m) */
double s; /* ２地点間の地表面距離(m) */
double w; /* ２地点間の半射程角(°) (中心角の１／２) */
double wr; /* ２地点間の半射程角(rad) */
double rad; /* 度→ラジアン変換係数 */
double N1,N2; /* 緯度補正した地球の半径(m) */
rad=PI/180.0;
h1=h2=0.0; /* ここでは、標高を無視 */
printf("第１地点の緯度、経度(°,西経は負符号)を入力して下さい。--->");
scanf("%lf %lf",&f1,&g1);
if(g1<0) g1=360.0+g1;
fr1=f1*rad; gr1=g1*rad;
printf("第２地点の緯度、経度(°,西経は負符号)を入力して下さい。--->");
scanf("%lf %lf",&f2,&g2);
if(g2<0) g2=360.0+g2;
fr2=f2*rad; gr2=g2*rad;
N1=a/(sqrt(1.0-e2*sin(fr1)*sin(fr1)));
x1=(N1+h1)*cos(fr1)*cos(gr1);
y1=(N1+h1)*cos(fr1)*sin(gr1);
z1=(N1*(1.0-e2)+h1)*sin(fr1);
N2=a/(sqrt(1.0-e2*sin(fr2)*sin(fr2)));
x2=(N2+h2)*cos(fr2)*cos(gr2);
y2=(N2+h2)*cos(fr2)*sin(gr2);
z2=(N2*(1.0-e2)+h2)*sin(fr2);
r=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2)); /* 直距離 */
wr=asin(r/2/a); /* 半射程角(rad) */
w=wr/rad; /* 半射程角(°) */
s=a*2*wr; /* 地表面距離 */
putchar('\n');
printf("第１地点の緯度：%lf 経度：%lf\n",f1,g1);
printf("第２地点の緯度：%lf 経度：%lf\n",f2,g2);
printf("直 距 離 ：%9.3lf km\n",r/1000);
printf("地表面距離 ：%9.3lf km\n",s/1000);
printf("半 射 程 角：%8.4lf °\n",w);
getchar(); getchar();
return 0;
}