begin comment JAZ164, R743, Outer Planets;
comment l̲i̲b̲r̲a̲r̲y̲ A0, A1, A4, A5, A12, A15;
integer form1p12e;
integer form1p1e;
integer form7p1;
integer form2p9;
integer k,t; real a,k2,x; boolean fi;
array y,ya,z,za[1:15],m[0:5],e[1:60],d[1:33];
array ownd[1:5,1:5],ownr[1:5];
real procedure f(k); integer k;
begin
integer i,j,i3,j3;
real p;
if k ≠ 1 then goto A;
for i ≔ 1 step 1 until 4 do
begin
i3 ≔ 3×i;
for j ≔ i+1 step 1 until 5 do
begin
j3 ≔ 3×j;
p ≔ (y[i3-2] - y[j3-2]) ⭡ 2 + (y[i3-1] - y[j3-1]) ⭡ 2 + (y[i3] - y[j3]) ⭡ 2;
ownd[i,j] ≔ ownd[j,i] ≔ 1/p/sqrt(p)
end
end ;
for i ≔ 1 step 1 until 5 do
begin
i3 ≔ 3×i;
ownd[i,i] ≔ 0;
p ≔ y[i3-2] ⭡ 2 + y[i3-1] ⭡ 2 + y[i3] ⭡ 2;
ownr[i] ≔ 1/p/sqrt(p)
end ;
A:
i ≔ (k - 1) ÷ 3 + 1;
f ≔ k2 × (- m[0] × y[k] × ownr[i] + SUM(j,1,5,m[j]×((y[3×(j-i)+k]-y[k])×ownd[i,j]-y[3×(j-i)+k]×ownr[j])))
end f;
procedure RK3n(x,a,b,y,ya,z,za,fxyj,j,e,d,fi,n); value b,fi,n;
integer j,n; real x,a,b,fxyj;
boolean fi; array y,ya,z,za,e,d;
begin
integer jj;
real xl,h,hmin,int,hl,absh,fhm,discry,discrz,toly,tolz,mu,mu1,fhy,fhz;
boolean last,first,reject;
array yl,zl,k0,k1,k2,k3,k4,k5[1:n],ee[1:4×n];
if fi
then begin d[3] ≔ a;
for jj ≔ 1 step 1 until n do
begin d[jj+3] ≔ ya[jj]; d[n+jj+3] ≔ za[jj] end
end ;
d[1] ≔ 0; xl ≔ d[3];
for jj ≔ 1 step 1 until n do
begin yl[jj] ≔ d[jj+3]; zl[jj] ≔ d[n+jj+3] end ;
if fi then d[2] ≔ b - d[3];
absh ≔ h ≔ abs(d[2]);
if b - xl < 0 then h ≔ - h;
int ≔ abs(b - xl); hmin ≔ int × e[1] + e[2];
for jj ≔ 2 step 1 until 2×n do
begin hl ≔ int × e[2×jj-1] + e[2×jj];
if hl < hmin then hmin ≔ hl
end ;
for jj ≔ 1 step 1 until 4×n do ee[jj] ≔ e[jj]/int;
first ≔ reject ≔ true ;
if fi
then begin last ≔ true ; goto nstep end ;
test: absh ≔ abs(h);
if absh < hmin
then begin h ≔ if h > 0 then hmin else - hmin;
absh ≔ hmin
end ;
if h > b - xl ≡ h > 0
then begin d[2] ≔ h; last ≔ true ;
h ≔ b - xl; absh ≔ abs(h)
end
else last ≔ false ;
nstep: if reject
then begin x ≔ xl;
for jj ≔ 1 step 1 until n do
y[jj] ≔ yl[jj];
for j ≔ 1 step 1 until n do
k0[j] ≔ fxyj × h
end
else begin fhy ≔ h/hl;
for jj ≔ 1 step 1 until n do
k0[jj] ≔ k5[jj] × fhy
end ;
x ≔ xl + ·27639 32022 50021 × h;
for jj ≔ 1 step 1 until n do
y[jj] ≔ yl[jj] + (zl[jj] × ·27639 32022 50021 +
k0[jj] × ·03819 66011 25011) × h;
for j ≔ 1 step 1 until n do k1[j] ≔ fxyj × h;
x ≔ xl + ·72360 67977 49979 × h;
for jj ≔ 1 step 1 until n do
y[jj] ≔ yl[jj] + (zl[jj] × ·72360 67977 49979 +
k1[jj] × ·26180 33988 74989) × h;
for j ≔ 1 step 1 until n do k2[j] ≔ fxyj × h;
x ≔ xl + h × ·5;
for jj ≔ 1 step 1 until n do
y[jj] ≔ yl[jj] + (zl[jj] × ·5 +
k0[jj] × ·04687 5 +
k1[jj] × ·07982 41558 39840 -
k2[jj] × ·00169 91558 39840) × h;
for j ≔ 1 step 1 until n do k4[j] ≔ fxyj × h;
x ≔ if last then b else xl + h;
for jj ≔ 1 step 1 until n do
y[jj] ≔ yl[jj] + (zl[jj] +
k0[jj] × ·30901 69943 74947 +
k2[jj] × ·19098 30056 25053) × h;
for j ≔ 1 step 1 until n do k3[j] ≔ fxyj × h;
for jj ≔ 1 step 1 until n do
y[jj] ≔ yl[jj] + (zl[jj] +
k0[jj] × ·08333 33333 33333 +
k1[jj] × ·30150 28323 95825 +
k2[jj] × ·11516 38342 70842) × h;
for j ≔ 1 step 1 until n do k5[j] ≔ fxyj × h;
reject ≔ false ; fhm ≔ 0;
for jj ≔ 1 step 1 until n do
begin
discry ≔ abs((- k0[jj] × ·5 + k1[jj] × 1·80901 69943 74947 +
k2[jj] × ·69098 30056 25053 - k4[jj] × 2) × h);
discrz ≔ abs((k0[jj] - k3[jj]) × 2 - (k1[jj] + k2[jj]) × 10 +
k4[jj] × 16 + k5[jj] × 4);
toly ≔ absh × (abs(zl[jj]) × ee[2×jj-1] + ee[2×jj]);
tolz ≔ abs(k0[jj]) × ee[2×(jj+n)-1] + absh × ee[2×(jj+n)];
reject ≔ discry > toly ∨ discrz > tolz ∨ reject;
fhy ≔ discry/toly; fhz ≔ discrz/tolz;
if fhz > fhy then fhy ≔ fhz;
if fhy > fhm then fhm ≔ fhy
end ;
mu ≔ 1/(1 + fhm) + ·45;
if reject
then begin if absh < hmin
then begin d[1] ≔ d[1] + 1;
for jj ≔ 1 step 1 until n do
begin y[jj] ≔ yl[jj];
z[jj] ≔ zl[jj]
end ;
first ≔ true ; goto next
end ;
h ≔ mu × h; goto test
end rej;
if first
then begin first ≔ false ; hl ≔ h; h ≔ mu × h;
goto acc
end ;
fhy ≔ mu × h/hl + mu - mu1; hl ≔ h; h ≔ fhy × h;
acc: mu1 ≔ mu;
for jj ≔ 1 step 1 until n do
z[jj] ≔ zl[jj] + (k0[jj] + k3[jj]) × ·08333 33333 33333 +
(k1[jj] + k2[jj]) × ·41666 66666 66667;
next: if b ≠ x
then begin xl ≔ x;
for jj ≔ 1 step 1 until n do
begin yl[jj] ≔ y[jj]; zl[jj] ≔ z[jj] end ;
goto test
end ;
if ¬ last then d[2] ≔ h;
d[3] ≔ x;
for jj ≔ 1 step 1 until n do
begin d[jj+3] ≔ y[jj]; d[n+jj+3] ≔ z[jj] end
end RK3n;
procedure TYP(x); array x;
begin
integer k;
newline(10, 1);
writetext(10,“T = ”);commentABSFIXT;
write(10,form7p1,t+a);
newline(10, 2);
for k ≔ 1 step 1 until 5 do
begin
if k=1 then writetext(10,“J ”)else
if k=2 then writetext(10,“S ”)else
if k=3 then writetext(10,“U ”)else
if k=4 then writetext(10,“N ”)else
writetext(10,“P ”);
write(10,form2p9,x[3×k-2]);
write(10,form2p9,x[3×k-1]);
write(10,form2p9,x[3×k]);
newline(10, 1)
end
end TYP;
real procedure SUM(i,a,b,xi); value b; integer i,a,b; real xi;
begin
real s;
s ≔ 0;
for i ≔ a step 1 until b do s ≔ s + xi;
SUM ≔ s
end SUM;
form1p12e ≔ format(“s+d.ddddddddddd@+nd”);
form1p1e ≔ format(“+d.d@+nd”);
form7p1 ≔ format(“snnnnnnd.d”);
form2p9 ≔ format(“+nd.ddddddddds”);
open(10);
open(20);
a ≔ read(20);
for k ≔ 1 step 1 until 15 do
begin
ya[k] ≔ read(20); za[k] ≔ read(20);
end ;
for k ≔ 0 step 1 until 5 do
m[k] ≔ read(20);
k2 ≔ read(20); e[1] ≔ read(20);
for k ≔ 2 step 1 until 60 do
e[k] ≔ e[1];
writetext(10,“JAZ164, R743, Outer Planets”);newline(10,2);
for k ≔ 1 step 1 until 15 do
begin
write(10,form1p12e,ya[k]);
write(10,form1p12e,za[k]);
newline(10, 1)
end ;
for k ≔ 0 step 1 until 5 do
begin
newline(10, 1);
write(10,form1p12e,m[k])
end ;
newline(10, 2);
write(10,form1p12e,k2);
newline(10, 2);
writetext(10,“eps = ”);
write(10,form1p1e,e[1]);
newline(10, 1);
t ≔ 0;
TYP(ya);
fi ≔ true ;
for t ≔ 500,1000 do
begin
RK3n(x,0,t,y,ya,z,za,f(k),k,e,d,fi,15);
fi ≔ false ;
TYP(y)
end;
close(20);
close(10);
end