农场阳光 (simpson)-CSDN博客

计算若干个圆与一个矩形的面积并

simpson公式 ans = ( f[l] + f[r] + 4 * f[mid] ) * (r-l) / 6

 1 uses math;
 2 type arr=record x,y:double; end;
 3 const eps=1e-8;
 4 var a,b,n,i:longint;
 5     sl,sr,sm:double;
 6     T:array[0..100] of record x,y,z,r:double; end;
 7     f:array[0..100] of arr;
 8     g:double;
 9 function sum(sl,sr,sm,l,r:double):double;
10 begin
11     exit((sl+sr+4*sm)*(r-l)/6);
12 end; 
13 procedure sort(l,r:longint);
14 var i,j:longint; x,y:double;
15     temp:arr;
16 begin
17     i:=l;j:=r;x:=f[(i+j) div 2].x; y:=f[(i+j) div 2].y;
18     while i<=j do
19     begin
20         while (f[i].x<x) or (f[i].x=x) and (f[i].y<y) do inc(i);
21         while (x<f[j].x) or (x=f[j].x) and (y<f[j].y) do dec(j);
22         if i<=j then
23         begin
24             temp:=f[i]; f[i]:=f[j]; f[j]:=temp;
25             inc(i); dec(j);
26         end;
27     end;
28     if l<j then sort(l,j);
29     if i<r then sort(i,r);
30 end;
31 function calc(xx:double):double;
32 var i,j,k:longint;
33     cnt,yy:double;
34 begin
35     k:=0;
36     for i:=1 to n do 
37         if abs(T[i].x-xx)+eps<T[i].r then
38         begin
39             yy:=sqrt(sqr(T[i].r)-sqr(T[i].x-xx));
40             inc(k);
41             f[k].x:=T[i].y-yy;
42             f[k].y:=T[i].y+yy;
43         end;
44     sort(1,k);
45     i:=1; cnt:=0;
46     //if xx=a/2 then writeln(f[1].x:4:2,f[1].y:4:2);
47     while i<=k do
48     begin
49         j:=i;
50         while (j<k) and (f[j+1].x<=f[i].y) do
51         begin
52             inc(j);
53             if f[j].y>f[i].y then f[i].y:=f[j].y;
54         end;
55         if not ((f[i].y<0) or (f[i].x>b)) then
56         begin
57             if (f[i].y<=b) and (f[i].x>=0) then cnt:=cnt+f[i].y-f[i].x else
58             if (f[i].y>b)  and (f[i].x<0)  then cnt:=cnt+b else
59             if (f[i].y>b) then cnt:=cnt+b-f[i].x else
60             if (f[i].x<0) then cnt:=cnt+f[i].y;
61         end;        
62         i:=j+1;
63     end;
64     exit(cnt);
65 end;
66 function simpson(l,r,mid,sl,sr,sm,w:double):double;
67 var sm1,sm2,m1,m2,w1,w2:double;
68 begin
69     m1:=(l+mid)/2;
70     m2:=(r+mid)/2;
71     sm1:=calc(m1);
72     sm2:=calc(m2);
73     w1:=sum(sl,sm,sm1,l,mid);
74     w2:=sum(sm,sr,sm2,mid,r);
75     if abs(w1+w2-w)<eps then exit(w);
76     exit(simpson(l,mid,m1,sl,sm,sm1,w1)+simpson(mid,r,m2,sm,sr,sm2,w2));
77 end;
78 begin
79     assign(input,'sun.in');reset(input);
80     assign(output,'sun.out');rewrite(output);
81     readln(a,b);
82     readln(g);
83     readln(n);
84     for i:=1 to n do 
85     begin
86         readln(T[i].x,T[i].y,T[i].z,T[i].r);
87         T[i].x:=T[i].x+cos(g/180*pi)/sin(g/180*pi)*T[i].z;
88     end;
89     sl:=calc(0); sr:=calc(a); sm:=calc(a/2);
90     //writeln(sl:4:2,' ',sr:4:2,' ',sm:4:2);
91     writeln(a*b-simpson(0,a,a/2,sl,sr,sm,sum(sl,sr,sm,0,a)):0:2);
92     close(input);
93     close(output);
94 end.