题目链接：http://61.187.179.132/JudgeOnline/problem.php?id=1502
题意：给出如下一棵分层的树，给出每层的高度和每个面的半径。光线是平行的，与地面夹角alpha。求树在地面上投影的面积。
首先，做这题需要知道一点：一个圆从任意一个角度投影都永远是一个圆。
我们可以画出一个简图如下：
如图，这棵树倒影之后，有图中两个圆心p1,p2,他们的横坐标即为这颗树上他们原先的高度乘以cotΘ，而他们的半价却不会变化，因为月光是平行光，所以在圆面与地面平行时，两点间距离不会变化。
如图，倒影最终是圆和他们之间的公切线构成的图形，最右边的点可以看做是半径为eps的圆。之后，可以利用simpson积分公式计算，simpson(l,r)=(f(l)+f(r)+4*f(mid))*(r-l)/6,若是精度差距大可以继续递归，注意：本题的eps要1e-6以下才能过。
1 #include<algorithm> 2 #include<cstdio> 3 #include<cmath> 4 #include<cstring> 5 #include<iostream> 6 const double eps=1e-6; 7 struct Point{ 8 double x,y; 9 Point(){}; 10 Point(double x0,double y0):x(x0),y(y0){}; 11 }e[200005],s[200005],a[200005],b[200005]; 12 int n; 13 double sqr(double x){ 14 return x*x; 15 } 16 void cal(Point &s,Point &e,Point a,Point b){ 17 if (std::fabs(a.y-b.y)<eps){ 18 s=a; 19 e=b; 20 return; 21 } 22 double x0=a.x-a.y*(b.x-a.x)/(b.y-a.y); 23 double k=a.y/(a.x-x0); 24 s.x=a.x-k*a.y; 25 s.y=sqrt(sqr(a.y)-sqr(a.x-s.x)); 26 e.x=b.x-k*b.y; 27 e.y=sqrt(sqr(b.y)-sqr(b.x-e.x)); 28 } 29 double f(double x){ 30 double y=0; 31 for (int i=1;i<=n+1;i++) 32 if (std::fabs(x-a[i].x)<=a[i].y) 33 y=std::max(y,sqrt(sqr(a[i].y)-sqr(std::fabs(x-a[i].x)))); 34 for (int i=1;i<=n;i++) 35 if (a[i+1].x-a[i].x-std::fabs(a[i].y-a[i+1].y)>eps&&x>=s[i].x&&x<=e[i].x){ 36 y=std::max(y,s[i].y+(x-s[i].x)*(e[i].y-s[i].y)/(e[i].x-s[i].x)); 37 } 38 return y; 39 } 40 double work(double L,double R){ 41 double mid=(L+R)/2; 42 return (f(L)+f(R)+f(mid)*4)*(R-L)/6; 43 } 44 double simpson(double L,double R){ 45 double mid=(L+R)/2; 46 double x1=work(L,mid),x2=work(mid,R),x3=work(L,R); 47 if (std::fabs(x1+x2-x3)<eps) return x1+x2; 48 else return simpson(L,mid)+simpson(mid,R); 49 } 50 int main(){ 51 double theta; 52 scanf("%d%lf",&n,&theta); 53 theta=1/(tan(theta)); 54 double h=0; 55 for (int i=1;i<=n+1;i++){ 56 scanf("%lf",&a[i].x); 57 h+=a[i].x; 58 a[i].x=h*theta; 59 } 60 for (int i=1;i<=n;i++) 61 scanf("%lf",&a[i].y); 62 a[n+1].y=a[n+2].y=0; 63 double L=a[1].x,R=a[n+1].x; 64 for (int i=1;i<=n;i++){ 65 L=std::min(a[i].x-a[i].y,L); 66 R=std::max(a[i].x+a[i].y,R); 67 if (a[i+1].x-a[i].x-std::fabs(a[i+1].y-a[i].y)>eps){ 68 cal(s[i],e[i],a[i],a[i+1]); 69 } 70 } 71 printf("%.2f ",2*simpson(L,R)); 72 }