ZOJ3663_Polaris of Pandora 2012 ICPC Changchun Site I 题

先是用经纬度算,调了很久,实在不会如何用经纬度求交点,百度谷歌也找不到。用三分法求了交点,有问题。。。

不得已各种三角函数反三角函数坐标旋转神马的都用上了,硬是求出三维坐标系的各种点再算,精度都不知损到哪里去了,竟然能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  }
posted @ 2012-10-17 17:24  CSGrandeur  阅读(656)  评论(2编辑  收藏  举报