先是用经纬度算，调了很久，实在不会如何用经纬度求交点，百度谷歌也找不到。用三分法求了交点，有问题。。。
不得已各种三角函数反三角函数坐标旋转神马的都用上了，硬是求出三维坐标系的各种点再算，精度都不知损到哪里去了，竟然能AC了。。。
1、算临界线的纬度
2、判断行程所在平面的“斜率”是否能与临界线相交，不相交则100.000
3、把行程平面的x方向旋转到x轴正方向算出与临界线交点，再在临界面旋转到行程面与临界面交点
4、用球面距离和z坐标结合判断计算区间。
1 #include<stdio.h> 2 #include<string.h> 3 #include<stdlib.h> 4 #include<math.h> 5 #include<vector> 6 #include<list> 7 #include<algorithm> 8 #include<iostream> 9 using namespace std; 10 const double eps = 1e-8; 11 const double pi = acos(-1.0); 12 inline int dcmp(double x){return (x > eps) - (x < -eps);} 13 inline double pz(double x) {return dcmp(x) ? x : 0;} 14 inline double Sqr(double x) {return x * x;} 15 inline double pcs(double x) {return x > 1 ? 1 : (x < -1 ? -1 : x);} 16 struct Point3 17 { 18 double x, y, z; 19 Point3(){x = y = z = 0;} 20 Point3(double a, double b, double c){x = a, y = b, z = c;} 21 Point3 cross(Point3 p){return 22 Point3(y * p.z - p.y * z, z * p.x - x * p.z, x * p.y - y * p.x);} 23 double dot(Point3 p){return x * p.x + y * p.y + z * p.z;} 24 Point3 operator-(const Point3 &p)const{return Point3(x - p.x, y - p.y, z - p.z);} 25 Point3 operator-()const{return Point3(-x, -y, -z);} 26 Point3 operator+(const Point3 &p)const{return Point3(x + p.x, y + p.y, z + p.z);} 27 Point3 operator*(const double &b)const{return Point3(x * b, y * b, z * b);} 28 Point3 operator/(const double &b)const{return Point3(x / b, y / b, z / b);} 29 Point3 fxl(Point3 b, Point3 c){return (*this - b).cross(b - c);} 30 double Dis(Point3 b){return sqrt((*this - b).dot(*this - b));} 31 double Rdis(Point3 b, double R){return R * asin(pcs((*this).Dis(b) * 0.5 / R)) * 2;} 32 double vlen(){return sqrt(dot(*this));} 33 Point3 RotePoint(const Point3 &p, double ang) 34 { 35 return Point3((p.x - x) * cos(ang) - (p.y - y) * sin(ang) + x, 36 (p.x - x) * sin(ang) + (p.y - y) * cos(ang) + y, p.z); 37 } 38 }; 39 double R, H, lat1, lng1, lat2, lng2, lat; 40 inline double p2cross(const Point3 &a, const Point3 &b, const Point3 &c) 41 {return (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y);} 42 Point3 r(0, 0, 0); 43 int main() 44 { 45 while(scanf("%lf%lf%lf%lf%lf%lf", &R, &H, &lat1, &lng1, &lat2, &lng2) != EOF) 46 { 47 double R_ = R; lat = acos(R / (R + H)); R += H; 48 lat1 += pi * 0.5, lat2 += pi * 0.5; 49 Point3 s(R * sin(lat1) * cos(lng1), R * sin(lat1) * sin(lng1), R * cos(lat1)); 50 Point3 e(R * sin(lat2) * cos(lng2), R * sin(lat2) * sin(lng2), R * cos(lat2)); 51 Point3 f = r.fxl(s, e); 52 double as = asin(pcs(fabs(f.cross(Point3(0, 0, 100)).vlen() / f.vlen() / 100))); 53 if(as < lat + eps) {printf("100.000\n"); continue;} 54 double len = R * sin(lat) / tan(as); 55 Point3 jd1(len, sqrt(Sqr(R_) - Sqr(len)), -R * sin(lat)); 56 Point3 jd2(len, -sqrt(Sqr(R_) - Sqr(len)), -R * sin(lat)); 57 if(f.z < -eps) f = -f; 58 double ang = atan2(f.y, f.x); 59 jd1 = r.RotePoint(jd1, ang), jd2 = r.RotePoint(jd2, ang); 60 double s1d = s.Rdis(jd1, R), s2d = s.Rdis(jd2, R); 61 double e1d = e.Rdis(jd1, R), e2d = e.Rdis(jd2, R); 62 double sed = s.Rdis(e, R), jd = jd1.Rdis(jd2, R); 63 if(!dcmp(s1d + e1d - sed) && !dcmp(s2d + e2d - sed)) 64 printf("%.3f\n", (sed - jd) / sed * 100); 65 else if(!dcmp(s1d + e1d - sed)) 66 { 67 if(e.z < s.z) printf("%.3f\n", s1d / sed * 100); 68 else printf("%.3f\n", e1d / sed * 100); 69 } 70 else if(!dcmp(s2d + e2d - sed)) 71 { 72 if(e.z < s.z) printf("%.3f\n", s2d / sed * 100); 73 else printf("%.3f\n", e2d / sed * 100); 74 } 75 else printf(e.z < jd1.z ? "0.000\n" : "100.000\n"); 76 } 77 return 0; 78 }