|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 12288 (0x3000)
Types: TextFile
Names: »tpbev1«
└─⟦2c55ea56f⟧ Bits:30001844 SW-save af projekt 1000, Alarm-system
└─⟦6b41451d2⟧
└─⟦this⟧ »tpbev1«
o pbevlist
mode list.yes
(pbev1=algol list.yes bossline.yes
if ok.yes
scope user pbev1
o c
convert pbevlist
lookup pbev1 tpbev1
finis)
\f
begin
procedure intgr(t,c,f,n,h,t1,t2,eps,print); value t,n,h,t1,t2,eps;
real t,h,t1,t2,eps; array c; procedure f,print; integer n;
begin
array d,d1,d2,c1,c2,b,b1,b2,a,a1,a2,p,p1,p2(1:n);
real norm; integer i; boolean dbl;
procedure l1(a,a1,a2,b,b1,b2,t); array a,a1,a2,b,b1,b2; real t;
begin integer i; real x;
for i:=1 step 1 until n do a(i):=b1(i)+b2(i)*h+b(i);
f(a2,a,t+h); norm:=0;
for i:=1 step 1 until n do begin
a(i):=x:=b1(i)+(b2(i)+a2(i))*h/2+b(i);
a1(i):=b(i)-x+(b2(i)+a2(i))*h/2+b1(i);
if abs(a2(i)-b2(i))/2 > norm then norm:=abs(a2(i)-b2(i))/2;
end;
end l1;
procedure l2(a,a1,a2,b,b1,b2,c,c1,c2,t);
array a,a1,a2,b,b1,b2,c,c1,c2; real t;
begin integer i; real x;
for i:=1 step 1 until n do a(i):=c1(i)+(9/4*b2(i)+3/4*c2(i))*h+c(i);
f(a2,a,t+h); norm:=0;
for i:=1 step 1 until n do begin
x:=27/32*(b(i)-c(i)+b1(i)-c1(i))
-(45/32*b2(i)+21/32*c2(i)-3/8*a2(i))*h;
a(i):=a(i)+x;
a1(i):=c(i)-a(i)+(9/4*b2(i)+3/4*c2(i))*h+x+c1(i);
if abs x > norm then norm:=abs x;
end;
norm:=norm/h;
end l2;
if t >= t1 then print(t,c,t1,t2);
if h=0 or t+3*h > t2 then h:=((t1+t2)/2-t)*1.5 else h:=h*2;
for i:=1 step 1 until n do c1(i):=0;
f(c2,c,t);
repeat
repeat
h:=h/2; l1(b,b1,b2,c,c1,c2,t);
until norm <= eps;
l1(a,a1,a2,b,b1,b2,t+h);
until norm <= eps;
dbl:=false; t:=t+2*h;
repeat
l2(p,p1,p2,a,a1,a2,c,c1,c2,t);
if norm > eps then begin
for i:=1 step 1 until n do begin
d(i):=c(i); d1(i):=c1(i); d2(i):=c2(i);
c(i):=b(i); c1(i):=b1(i); c2(i):=b2(i);
end;
h:=h/2; l2(b,b1,b2,c,c1,c2,d,d1,d2,t-2*h); dbl:=false;
end else begin
t:=t+h; if t >= t1 then print(t,p,t1,t2);
if t < t2 then begin
while t+h > t2 do begin
for i:=1 step 1 until n do begin
d(i):=b(i); d1(i):=b1(i); d2(i):=b2(i);
b(i):=a(i); b1(i):=a1(i); b2(i):=a2(i);
a(i):=c(i); a1(i):=c1(i); a2(i):=c2(i);
end;
h:=h/2;
l2(c,c1,c2,d,d1,d2,a,a1,a2,t-4*h);
l2(a,a1,a2,b,b1,b2,d,d1,d2,t-2*h);
norm:=eps;
end t+h > t2;
if norm < eps/16 and t+2*h <= t2 and dbl then begin
for i:=1 step 1 until n do begin
a(i):=p(i); a1(i):=p1(i); a2(i):=p2(i);
c(i):=d(i); c1(i):=d1(i); c2(i):=d2(i);
end;
h:=2*h; dbl:=false;
end else begin
for i:=1 step 1 until n do begin
d(i):=c(i); d1(i):=c1(i); d2(i):=c2(i);
c(i):=b(i); c1(i):=b1(i); c2(i):=b2(i);
b(i):=a(i); b1(i):=a1(i); b2(i):=a2(i);
a(i):=p(i); a1(i):=p1(i); a2(i):=p2(i);
end;
dbl:=true;
end eps/16 < norm < eps;
end t<t2;
end norm < eps;
until t>= t2;
end intgr;
procedure decompose(a,b,n); value n;
array a; integer array b; integer n;
begin integer r,s,i,j,p;
real ai,ar,np,nr,ars;
np:=0; i:=j:=0;
repeat
if np=0 then begin
if j>0 then b(j):=0;
j:=j+1;
if j<n then begin
for r:=i+1 step 1 until n do begin
nr:=0;
for s:=j step 1 until n do nr:=nr+abs a(r,s);
if abs a(r,j)>nr*np then begin
p:=r; np:=abs a(r,j)/nr;
end;
end;
end;
end else begin
i:=i+1;
for s:=n step -1 until j do begin
ai:=a(p,s); a(p,s):=a(i,s); a(i,s):=ai;
end;
b(j):=p; np:=0; j:=j+1;
for r:=i+1 step 1 until n do begin
ar:=a(r,j-1):=a(r,j-1)/ai; nr:=0;
for s:=n step -1 until j do begin
ars:=a(r,s):=a(r,s)-a(i,s)*ar;
nr:=nr+abs ars;
end s;
if abs ars>nr*np then begin
p:=r; if nr>0 then np:=abs ars/nr;
end;
end r;
end np;
until j>n;
end;
procedure solve(a,b,n,x); value n;
array a,x; integer array b; integer n;
begin integer r,s,i; real v;
i:=1;
for s:=1 step 1 until n do if b(s)<>0 then begin
v:=x(b(s)); x(b(s)):=x(i); x(i):=v; i:=i+1;
for r:=i step 1 until n do x(r):=x(r)-v*a(r,s);
end;
for r:=n step -1 until 1 do if b(r)<>0 then begin
i:=i-1; v:=x(i);
for s:=r+1 step 1 until n do v:=v-x(s)*a(i,s);
x(r):=v/a(i,r);
end else x(r):=0;
end;
<*test*******************************************************>
procedure o(n,v); value n,v; integer n; real v;
if n<test then
write(out,<:<10>:>,<<ddd>,n,<< -.ddd ddd ddd ddd'+zdd>,v);<************>
integer n,n1,n2,n3,n6,n16,n26,i,j,k,i6,no,test;
boolean prt;
real g,alfa,eps,delta,a,pi2,u0,u1,u2,t0,tx,to,tr,dts,dto,dt;
zone zo(128,1,stderror);
<* mainprogram **********************************************************>
read(in,test,g,alfa,dts,dt,eps,delta,n,n2,n1,n3,a);
open(zo,4,<:pbout:>,0); outrec6(zo,12*4);
for i:=1 step 1 until 12 do zo(i):=case i of
(test,g,alfa,dts,dt,eps,delta,n,n2,n1,n3,a);
n16:=-n1*6; n6:=n*6; n2:=n2+n; n26:=n2*6;
begin array m(-n1:n2),p,p0(1:n2,1:8),x(n16+1:n26+6),y(1:n2,1:7,0:7),
x1(1:3,0:7),p1(1:8),
p2(if n2<1 then -1 else -n2:-1,1:if n3<1 then 1 else n3);
procedure grv(x1,x,t); array x1,x; real t;
begin
integer i,j,k; real r; array a(1:3);
for i:=n16 step 6 until -6 do orb1(t,x,i);
for i:=0 step 6 until n26 do begin
for k:=4,5,6 do begin x1(i+k-3):=x(i+k); x1(i+k):=0 end;
for j:=n16 step 6 until n6 do if i<>j then begin
r:=0;
for k:=1,2,3 do begin
a(k):=x(j+k)-x(i+k); r:=r+a(k)**2;
end;
r:=g*m(j//6)/r/sqrt(r);
for k:=4,5,6 do x1(i+k):=x1(i+k)+r*a(k-3);
end;
end;
end;
procedure orb(t,x,p); value t; real t; array x<*1:3,0:7*>,p<*1:8*>;
<*global real delta,pi2*>
begin integer j;
real s,y,e,e1,sy,cy,r,q,q1,sv,cv,ex,ey,a;
array r1,y1(1:4);
procedure dif(i); integer i;
begin
x(i,7):=q/a; q:=q/r;
for j:=1 step 1 until 4 do x(i,j):=q*r1(j)+q1*y1(j);
end;
q:=p(1); j:=t/q; e:=q+q*j-q*j; q1:=pi2/q;
ex:=p(3); ey:=p(4); r:=arg(ex,ey);
y:=(t-e*j-(q-e)*j-p(8)*j)*q1+p(2)-r;
s:=(round(y/pi2+.5)-.5)*pi2; cy:=sign(y-s);
e:=sqrt(ex**2+ey**2); e1:=sqrt(1-e**2);
repeat
sy:=s;
s:=((sin(s)-s*cos(s))*e+y)/(1-e*cos(s));
until cy*(s-sy)<=delta;
q1:=-q1*t/q; sy:=sin(s); cy:=cos(s);
sv:=sin(s+r); cv:=cos(s+r);
y:=arg(cy-e,e1*sy)+r;
r:=1-e*cy; r1(2):=q:=e*sy/r; r1(1):=q*q1; r1(3):=(ex-cv)/r;
r1(4):=(ey-sv)/r;
y1(2):=q:=e1/r**2; y1(1):=q*q1; e:=sy*(cy-e)/e1/r**2;
y1(3):=(2*sv-ey*cy**2)*q-ey/(1+e1)-ex*e;
y1(4):=ex/(1+e1)-(2*cv-ex*cy**2)*q-ey*e;
a:=p(7); e1:=a*r; sy:=sin(y); cy:=cos(y);
ex:=p(5); ey:=p(6); e:=(ex**2+ey**2+1)/2;
s:=e1*(ex*cy+ey*sy)/e; sv:=e1*(ex*sy-ey*cy)/e;
x(1,0):=q:=e1*cy-ex*s; q1:=-e1*sy+ex*sv;
x(3,5):=cv:=-q/e; x(1,5):=cv*ex-s; x(2,5):=cv*ey; dif(1);
x(2,0):=q:=e1*sy-ey*s; q1:=e1*cy+ey*sv;
x(3,6):=cv:=-q/e; x(1,6):=cv*ex; x(2,6):=cv*ey-s; dif(2);
x(3,0):=q:=-s; q1:=sv; dif(3);
end;
procedure orbinit(t,x,p,i); real t; array x,p; integer i;
begin array y(1:3,0:7); integer j;
orb(t,y,p);
for j:=1,2,3 do begin
x(i+j):=y(j,0); x(i+j+3):=y(j,2)*pi2/p(1);
end;
end;
procedure orb1(t,x,i);
real t; array x; integer i;
begin end;
procedure orbinit1(t,x,i);
real t; array x; integer i;
begin end;
procedure pbevprt(t,x,t1,t2); real t,t1,t2; array x;
begin integer pn,pn6,i,j; real g,g1,g2;
array p1(1:8),x1(1:3,0:7),z(1:7,1:7),z1(1:7);
integer array z2(1:7);
if prt then begin
outrec6(zo,4); zo(1):=t;
end;
g2:=(t0-t)/alfa/4; g1:=exp(g2*4); g:=1-g1;
for pn:=1 step 1 until n2 do begin
if prt then outrec6(zo,13*4);
pn6:=pn*6; p1(1):=p0(pn,1);
for i:=2 step 1 until 8 do p1(i):=p(pn,i)+p0(pn,i);
orb(t,x1,p1);
<*test********************************************************>
for i:=1 step 1 until 3 do begin
o(pn6+i,x(pn6+i)); o(pn6+3+i,x(pn6+3+i)); o(i,x(i)); o(3+i,x(3+i));
for j:=0 step 1 until 7 do o(10*i+j,x1(i,j));
end;
for i:=1 step 1 until 8 do o(i,p1(i));
for i:=1 step 1 until 3 do begin
if prt then zo(i*2-1):=x1(i,0);
x1(i,0):=x1(i,0)-x(pn6+i)+x(i);
if prt then zo(i*2):=x1(i,0);
end;
for i:=1 step 1 until 7 do begin
for j:=0 step 1 until i do begin
y(pn,i,j):=y(pn,i,j)*g1+(x1(1,i)*x1(1,j)+x1(2,i)*x1(2,j)
+x1(3,i)*x1(3,j))*g;
<*test****************************************************************>
o(100+10*j+i,y(pn,i,j));
if j=0 then z1(i):=y(pn,i,0) else
z(i,j):=z(j,i):=y(pn,i,j);
end;
end;
decompose(z,z2,7); solve(z,z2,7,z1);
<*test*******************************************************************>
for i:=1 step 1 until 7 do begin o(i,z1(i)); o(200+i,z2(i));
for j:=1 step 1 until 7 do o(200+i*10+j,z(i,j));
end;
p(pn,8):=p(pn,8)+g2*z1(1);
if prt then zo(7):=p(pn,8);
for i:=2 step 1 until 7 do begin
p(pn,i):=p(pn,i)+g2*z1(i);
if prt then zo(i+6):=p(pn,i);
end;
end;
t0:=t;
if prt then begin
if no=0 then begin
read(in,test,no,tr); dto:=(tr-tx)/no;
end;
no:=no-1; to:=tr-no*dto;
end;
tx:=tx+dts; prt:= tx>=to;
if prt then begin
tx:=to; t1:=tx-dt; t2:=tx+dt;
end else begin
t1:=tx-dts/10; t2:=tx+dts/10;
end;
end;
<*initialize *******************************************************>
pi2:=2*3.1415926536;
for j:= 1 step 1 until 6 do x(j):=0;
m(0):=a;
for i:=1 step 1 until n2 do begin
read(in,m(i),u0,u1,p(i,8));
u2:=p0(i,1):=u0+u1; p0(i,8):=u0-u2+u1;
for j:=2 step 1 until 7 do read(in,p0(i,j),p(i,j));
a:=a+m(i);
outrec6(zo,16*4); zo(1):=m(i); zo(2):=p0(i,0);
for j:= 2 step 1 until 8 do begin
zo(j+1):=p0(i,j); zo(j+8):=p(i,j);
end;
end;
for i:=-n1 step 1 until -1 do begin
read(in,m(i));
for j:=1 step 1 until n3 do read(in,p2(i,j));
a:=a+m(i);
outrec6(zo,(1+n3)*4); zo(1):=m(i);
for j:=1 step 1 until n3 do zo(j+1):=p2(i,j);
end;
outrec6(zo,0);
read(in,tx); prt:=true; no:=0;
for i:=1 step 1 until n2 do begin
u0:=p1(1):=p0(i,8)+p(i,8)+p0(i,1);
p1(8):=p0(i,1)-u0+p(i,8)+p0(i,8);
for j:= 2 step 1 until 7 do p1(j):=p0(i,j)+p(i,j);
i6:=i*6; orbinit(tx,x,p1,i6);
for j:= 1 step 1 until 6 do x(j):=x(j)-m(i)*x(i6+j)/a;
end;
for i:=-n1 step 1 until -1 do begin
i6:=i*6; orbinit1(tx,x,i);
for j:=1 step 1 until 6 do x(j):=x(j)-m(i)*x(i6+j)/a;
end;
for i:=1 step 1 until n2 do begin
i6:=i*6;
for j:=1 step 1 until 6 do x(i6+j):=x(i6+j)+x(j);
for j:=1 step 1 until 7 do begin
p(i,j+1):=0;
for k:=0 step 1 until j do y(i,j,k):=0;
end;
end;
t0:=tx-alfa; to:=t0+dts;
repeat
u1:=exp((t0-to)/alfa); u2:=1-u1;
for k:=1 step 1 until n2 do begin
for i:=1 step 1 until 8 do p1(i):=p0(k,i);
orb(to,x1,p1);
for i:=1 step 1 until 7 do begin
for j:=1 step 1 until i do y(k,i,j):=y(k,i,j)*u1
+(x1(1,i)*x1(1,j)+x1(2,i)*x1(2,j)+x1(3,i)*x1(3,j))*u2;
end;
end;
t0:=to; to:=to+dts;
until to>=tx;
intgr(tx,x,grv,n26+6,dts,tx-dt,tx+dt,eps,pbevprt);
outrec6(zo,4); zo(1):=tr;
close(zo,false);
end;
end;
▶EOF◀