這么簡(jiǎn)單的幾何題居然讓我WA了3次，果然一日不編程即手生矣...
要點(diǎn)： （1） 經(jīng)緯度－>球面坐標(biāo)－>直角坐標(biāo)
（2） 求的是球面距離
（3） 計(jì)算大圓圓心角的時(shí)候，要注意三角形不是以R＝6370為腰的等腰三角形，而是以A，B，O三點(diǎn)各自距離為邊的三角形
（4） 計(jì)算圓心角時(shí)用點(diǎn)乘會(huì)方便很多（還是看到Feli的才發(fā)現(xiàn)，誒）
2007.10.25
Simbaforrest
1 #include<stdio.h>
2 #include<math.h>
3 const double R=6370;
4 const double pi=3.14159265358979323846264338327950;
5 typedef struct tagAngle
6 {
7 double d,m;
8 double ang;
9 char s[20];
10 }Angle;
11 typedef struct tagPoint
12 {
13 Angle a,b;
14 double x,y,z;
15 }Point;
16
17 Point self,dest;
18
19 void change(Point &pt)
20 {
21 pt.a.ang=pt.a.d+pt.a.m/60;
22 pt.b.ang=pt.b.d+pt.b.m/60;
23 //change lat
24 if(pt.a.s[0]=='N')
25 {
26 pt.a.ang=90-pt.a.ang;
27 }
28 else
29 {
30 pt.a.ang+=90;
31 }
32 //change lot;
33 if(pt.b.s[0]=='W')
34 {
35 pt.b.ang=360-pt.b.ang;
36 }
37 pt.z=R*cos(pt.a.ang*pi/180);
38 pt.x=R*sin(pt.a.ang*pi/180)*cos(pt.b.ang*pi/180);
39 pt.y=R*sin(pt.a.ang*pi/180)*sin(pt.b.ang*pi/180);
40 }
41
42 bool init()
43 {
44 if(scanf("%lf",&self.a.d)!=EOF)
45 {
46 scanf("%lf%s%lf%lf%s",&self.a.m,self.a.s,&self.b.d,&self.b.m,self.b.s);
47 change(self);
48 scanf("%lf%lf%s%lf%lf%s",&dest.a.d,&dest.a.m,dest.a.s,&dest.b.d,&dest.b.m,dest.b.s);
49 change(dest);
50 return 1;
51 }
52 return 0;
53 }
54
55 double dist(Point a)
56 {
57 return sqrt( (a.x)*(a.x)+(a.y)*(a.y)+(a.z)*(a.z) );
58 }
59
60 double cross(Point a,Point b)
61 {
62 return (a.x*b.x+a.y*b.y+a.z*b.z);
63 }
64
65 void solve()
66 {
67 double ans,thita;
68 thita=acos( (cross(self,dest))/dist(self)/dist(dest) );
69 printf("%.3lf\n",R*thita);
70 }
71
72 int main()
73 {
74 while(init())
75 {
76 solve();
77 }
78 return 0;
79 }
80
posted on 2007-10-25 17:00
R2 閱讀(170)
評(píng)論(0) 編輯 收藏 引用