1 {$N+}
2 uses crt,graph;
3 const
4 eps: double=1e-8{cm/s};
5 del: double=1e-5{cm};
6 alp: double=0.975;
7 height: double=8{cm};
8 g: double=980{cm/s^2};
9 d: double=0.1{cm};
10 w: double=200{s-1};
11 A: double=-0.2{cm};
12 n=69;
13 Nimpact: longint=500000;
14 type
15 mas1=array[1..100] of double;
16 mas2=array[1..100] of integer;
17 var
18 xn,vn,AddXn,Addvn2 : mas1;
19 Xnav,Vn2av : mas1;
20 t,delta01,delta : double;
21 ygr1,ygr2,ygr3,ygr4 : mas2;
22 k,rad : integer;
23 l : longint;
24 maxx,maxy,ygrp1,ygrp2 : integer;
25 title1 : string;
26 ps0,ps1,pc0,pc1 : pointer;
27 {---------------------------------------------------------------------}
x0null28 function x0(t : double) : double;
29 begin
30 x0:=A*sin(w*t);
31 end;
32
v0null33 function v0(t : double) : double;
34 begin
35 v0:=A*w*cos(w*t);
36 end;
37 {---------------------------------------------------------------------}
38 procedure init;
39 var
40 i : integer;
41 begin
42 title1:=' n <Xn> <Vn^2> Vn Xn';
43 t:=0;
44 for i:=1 to n do
45 begin
46 xn[i]:=d/2+2*del+2*(d/2+2*del)*(i-1);
47 vn[i]:=0;
48 AddXn[i]:=0;
49 AddVn2[i]:=0;
50 end;
51 end;
52 {-----------------------------------------------------------------}
53 procedure initgr;
54 var
55 gd,gm : integer;
56 error : integer;
57 begin
58 gd:=detect;
59 initgraph(gd,gm,'D:\BP\BGI');
60 error:=graphresult;
61 if error<>grOK then
62 begin
63 writeln(grapherrormsg(error));
64 readln;
65 halt;
66 end;
67 cleardevice;
68 setbkcolor(blue);
69 end{initgr};
70 {-------------------------------------------------------------------}
71 procedure newton(tt1 : double; var t2 : double);
72 var
73 f,derf : double;
74 p : integer;
75 begin
76 for p:=1 to 10 do
77 begin
78 f:=-x0(t+tt1)+xn[1]+vn[1]*tt1-g*sqr(tt1)/2-d/2;
79 derf:=-v0(t+tt1)+vn[1]-g*tt1;
80 tt1:=tt1-f/derf;
81 end;
82 t2:=tt1;
83 end{newton};
84 {--------------------------------------------------------------------}
85 procedure predict1;
86 var
87 tf,t1,dt : double;
88 dx0,dx0init : double;
89 begin
90 t1:=0;
91 dx0init:=-x0(t)+xn[1]-d/2;
92 dt:=1/(1000*w);
93 repeat
94 begin
95 repeat
96 t1:=t1+dt;
97 dx0:=-x0(t+t1)+xn[1]+vn[1]*t1-g*sqr(t1)/2-d/2;
98 until dx0<0;
99 newton(t1,delta01);
100 end;
101 until delta01>0;
102 end{predict1};
103 {---------------------------------------------------------------------}
104 procedure change;
105 var
106 vnold : double;
107 begin
108 if k=0 then
109 begin
110 vn[1]:=-alp*vn[1]+2*v0(t);
111 xn[1]:=x0(t)+d/2+del;
112 end
113 else
114 begin
115 vnold:=vn[k+1];
116 vn[k+1]:=-1/2*((alp-1)*vn[k+1]+vn[k]*(-1-alp));
117 vn[k]:=1/2*((alp+1)*vnold+vn[k]*(1-alp));
118 end;
119 end{change};
120 {---------------------------------------------------------------------}
121 procedure re_count;
122 var
123 i : integer;
124 begin
125 if delta>0 then
126 begin
127 for i:=1 to n do
128 begin
129 xn[i]:=xn[i]+vn[i]*delta-g*sqr(delta)/2;
130 vn[i]:=vn[i]-g*delta;
131 end;
132 for i:=2 to n do
133 begin
134 if xn[i]<(xn[i-1]+d) then xn[i]:=xn[i-1]+d+del;
135 end;
136 t:=t+delta;
137 change;
138 if (k=1) or (k=0) then
139 predict1
140 else
141 delta01:=abs(delta01-delta)
142 end;
143 end{re-count};
144 {---------------------------------------------------------------------}
145 procedure predict;
146 var
147 delta1 : double;
148 i : integer;
149 begin
150 delta:=1e25;
151 for i:=1 to n-1 do
152 begin
153 if vn[i]<>vn[i+1] then
154 delta1:=(abs(xn[i+1]-xn[i]-d))/(vn[i]-vn[i+1])
155 else
156 delta1:=1e25;
157 if (delta1>0) and (delta1<delta) then
158 begin
159 delta:=delta1;
160 k:=i;
161 end;
162 end;
163 if (k>0) and (abs(vn[k+1]-vn[k])<eps) and (abs(xn[k+1]-xn[k]-d)<del/10) then
164 begin
165 vn[k+1]:=vn[k];
166 delta:=0;
167 end;
168 if delta>delta01 then
169 begin
170 delta:=delta01;
171 k:=0;
172 end;
173 end{predict};
174 {----------------------------------------------------------------}
175 procedure mean;
176 var
177 i : integer;
178 begin
179 for i:=1 to n do
180 begin
181 AddXn[i]:=AddXn[i]+Xn[i]*delta;
182 AddVn2[i]:=AddVn2[i]+sqr(Vn[i])*delta;
183 Xnav[i]:=AddXn[i]/t;
184 Vn2av[i]:=AddVn2[i]/t;
185 end;
186 end;{mean}
187 {-----------------------------------------------------------------}
188 procedure drawax;
189 var
190 i,j : integer;
191 size : integer;
192 s1 : string;
193 begin
194 setcolor(white);
195 maxx:=getmaxx;
196 maxy:=getmaxy;
197 size:=Imagesize(maxx-220,maxy-50,maxx-130,maxy-40);
198 getmem(ps0,size);
199 getimage(maxx-220,maxy-50,maxx-130,maxy-40,ps0^);
200 setlinestyle(solidln,$AAA,normwidth);
201 rectangle(6,6,maxx-6,maxy-6);
202 rectangle(3,3,maxx-3,maxy-3);
203 settextstyle(defaultfont,horizdir,1);
204 line(maxx-120,50,maxx-120,maxy-50);
205 line(maxx-220,maxy-50,maxx-20,maxy-50);
206 for j:=0 to 10 do
207 line(maxx-112,round((maxy-100)*j/10)+50,maxx-128,round((maxy-100)*j/10)+50);
208 for j:=1 to 20 do
209 line(round(10*j)+maxx-220,maxy-50,round(10*j)+maxx-230,maxy-40);
210 str(height:8,s1);
211 s1:='H='+s1+' cm';
212 settextstyle(defaultfont,vertdir,1);
213 outtextxy(maxx-95,50,s1);
214 settextstyle(defaultfont,horizdir,1);
215 ygrp1:=maxy-50;
216 getmem(ps1,size);
217 getimage(maxx-220,maxy-50,maxx-130,maxy-40,ps1^);
218
219 rad:=round(d/(2*height)*(maxy-100));
220 size:=Imagesize(maxx-200-rad,maxy-100-rad,maxx-200+rad,maxy-100+rad);
221 getmem(pc0,size);
222 getimage(maxx-200-rad,maxy-100-rad,maxx-200+rad,maxy-100+rad,pc0^);
223 setcolor(red);
224 setfillstyle(1,12);
225 fillellipse(maxx-200,maxy-100,rad,rad);
226 getmem(pc1,size);
227 getimage(maxx-200-rad,maxy-100-rad,maxx-200+rad,maxy-100+rad,pc1^);
228 putimage(maxx-200-rad,maxy-100-rad,pc0^,0);
229 for i:=1 to n do
230 begin
231 ygr1[i]:=round(maxy-50-(maxy-100)*Xn[i]/height-rad);
232 ygr2[i]:=round(maxy-50-(maxy-100)*Xn[i]/height-rad);
233 putimage(maxx-180,ygr1[i],pc1^,0);
234 putimage(maxx-70,ygr2[i],pc1^,0);
235 end;
236 end{drawax};
237 {--------------------------------------------------------------------}
238 procedure table( k1 : longint);
239 var
240 i,j : integer;
241 s,s1,s2 : string;
242 begin
243 setlinestyle(solidln,$AAA,normwidth);
244 setfillstyle(1,1);
245 setcolor(15);
246 if n>30 then j:=30 else j:=n;
247 bar3d(30,25,385,25+12*(j+3)+10,10,TopOn);
248 str(k1,s);
249 s:='Collision number = '+s;
250 str(t:12:5,s1);
251 outtextxy(40,30,s);
252 str(k:3,s2);
253 s:='Time='+s1+'s K = '+s2;
254 outtextxy(40,42,s);
255 outtextxy(40,54,title1);
256 for i:=1 to j do
257 begin
258 str(i:3,s);
259 str(Xnav[i]:8:5,s1);
260 str(Vn2av[i]:8:3,s2);
261 s:=s+' '+s1+' '+s2;
262 str(Vn[i]:8:3,s1);
263 str(xn[i]:8:5,s2);
264 s:=s+' '+s1+' '+s2;
265 outtextxy(40,54+(i)*12,s);
266 end;
267 end;{table}
268 {----------------------------------------------------------------}
269 procedure drawgr;
270 var
271 j,i : integer;
272 s,s1,s2 : string;
273 begin
274 setlinestyle(solidln,$AAA,normwidth);
275 setfillstyle(1,1);
276 setcolor(15);
277 for i:=1 to n do
278 begin
279 if abs(xnav[i])<(1.1*height) then
280 begin
281 ygr3[i]:=round(maxy-50-(maxy-100)*Xn[i]/height-rad);
282 ygr4[i]:=round(maxy-50-(maxy-100)*Xnav[i]/height-rad);
283 end;
284 if ygr3[i]<>ygr1[i] then
285 begin
286 putimage(maxx-180,ygr1[i],pc0^,0);
287 putimage(maxx-180,ygr3[i],pc1^,0);
288 ygr1[i]:=ygr3[i];
289 end
290 end;
291 for i:=1 to n do
292 begin
293 if ygr4[i]<>ygr2[i] then
294 begin
295 putimage(maxx-70,ygr2[i],pc0^,0);
296 putimage(maxx-70,ygr4[i],pc1^,0);
297 ygr2[i]:=ygr4[i];
298 end;
299 end;
300 ygrp2:=round(maxy-50-(maxy-100)*X0(t)/height);
301 if ygrp2<>ygrp1 then
302 begin
303 putimage(maxx-220,ygrp1,ps0^,0);
304 putimage(maxx-220,ygrp2,ps1^,0);
305 ygrp1:=ygrp2;
306 end;
307 end;{drawgr}
308 {----------------------------------------------------------------}
309 BEGIN
310 init;
311 initgr;
312 drawax;
313 predict1;
314 l:=1;
315 repeat
316 begin
317 predict;
318 re_count;
319 mean;
320 if l mod 2 =0 then
321 drawgr;
322 if l mod 1000 =0 then
323 table(l);
324 l:=l+1;
325 end;
326 until KeyPressed;
327 closegraph;
328 END.
329