begin comment JAZ164, R743, Outer Planets;
library 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]) !up 2 + (y[i3-1] - y[j3-1]) !up 2 + (y[i3] - y[j3]) !up 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] !up 2 + y[i3-1] !up 2 + y[i3] !up 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 eqv 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 or discrz > tolz or 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 not 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 _ = _ }); comment ABSFIXT;
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
|
2430000.5,
+0.342947415189º+1,
-0.557160570446º-2,
+0.335386959711º+1,
+0.505696783289º-2,
+0.135494901715º+1,
+0.230578543901º-2,
+0.664145542550º+1,
-0.415570776342º-2,
+0.597156957878º+1,
+0.365682722812º-2,
+0.218231499728º+1,
+0.169143213293º-2,
+0.112630437207º+2,
-0.325325669158º-2,
+0.146952576794º+2,
+0.189706021964º-2,
+0.627960525067º+1,
+0.877265322780º-3,
-0.301552268759º+2,
-0.240476254170º-3,
+0.165699966404º+1,
-0.287659532608º-2,
+0.143785752721º+1,
-0.117219543175º-2,
-0.211238353380º+2,
-0.176860753121º-2,
+0.284465098142º+2,
-0.216393453025º-2,
+0.153882659679º+2,
-0.148647893090º-3,
+0.100000597682º+1,
+0.954786104043º-3,
+0.285583733151º-3,
+0.437273164546º-4,
+0.517759138449º-4,
+0.277777777778º-5,
+0.295912208286º-3,
+0.10º-3;
|