code 33013 ;
procedure RK2N(X, A, B, Y, YA, Z, ZA, FXYZJ, J, E, D,
FI, N); value B, FI, N; integer J, N; real X, A, B, FXYZJ;
boolean FI; array Y, YA, Z, ZA, E, D;
begin integer JJ;
real XL, H, INT, HMIN, 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:= true; if FI then
begin LAST:= true; goto STEP end;
TEST: ABSH:= ABS(H); if ABSH < HMIN then
begin H:= if H > 0 then HMIN else - HMIN;
ABSH:= ABS(H)
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;
STEP: X:= XL;
for JJ:= 1 step 1 until N do
begin Y[JJ]:= YL[JJ]; Z[JJ]:= ZL[JJ] end;
for J:= 1 step 1 until N do K0[J]:= FXYZJ * H;
X:= XL + H / 4.5;
for JJ:= 1 step 1 until N do
begin Y[JJ]:= YL[JJ] + (ZL[JJ] * 18 + K0[JJ] * 2) /
81 * H; Z[JJ]:= ZL[JJ] + K0[JJ] / 4.5;
for J:= 1 step 1 until N do K1[J]:= FXYZJ * H;
X:= XL + H / 3;
for JJ:= 1 step 1 until N do
begin Y[JJ]:= YL[JJ] + (ZL[JJ] * 6 + K0[JJ]) / 18 * H;
Z[JJ]:= ZL[JJ] + (K0[JJ] + K1[JJ] * 3) / 12
end;
for J:= 1 step 1 until N do K2[J]:= FXYZJ * H;
X:= XL + H * .5;
for JJ:= 1 step 1 until N do
begin Y[JJ]:= YL[JJ] + (ZL[JJ] * 8 + K0[JJ] + K2[JJ])
/ 16 * H;
Z[JJ]:= ZL[JJ] + (K0[JJ] + K2[JJ] * 3) / 8
end;
for J:= 1 step 1 until N do K3[J]:= FXYZJ * H;
X:= XL + H * .8;
for JJ:= 1 step 1 until N do
begin Y[JJ]:= YL[JJ] + (ZL[JJ] * 100 + K0[JJ] * 12 +
K3[JJ] * 28) / 125 * H;
Z[JJ]:= ZL[JJ] + (K0[JJ] * 53 - K1[JJ] * 135 +
K2[JJ] * 126 + K3[JJ] * 56) / 125
end;
for J:= 1 step 1 until N do K4[J]:= FXYZJ * H;
X:= if LAST then B else XL + H;
for JJ:= 1 step 1 until N do
begin Y[JJ]:= YL[JJ] + (ZL[JJ] * 336 + K0[JJ] * 21 +
K2[JJ] * 92 + K4[JJ] * 55) / 336 * H;
Z[JJ]:= ZL[JJ] + (K0[JJ] * 133 - K1[JJ] * 378 +
K2[JJ] * 276 + K3[JJ] * 112 + K4[JJ] * 25) / 168
end;
for J:= 1 step 1 until N do K5[J]:= FXYZJ * H;
REJECT:= false; FHM:= 0;
for JJ:= 1 step 1 until N do
begin DISCRY:= ABS(( - K0[JJ] * 21 + K2[JJ] * 108 -
K3[JJ] * 112 + K4[JJ] * 25) / 56 * H);
DISCRZ:= ABS(K0[JJ] * 21 - K2[JJ] * 162 + K3[JJ]
* 224 - K4[JJ] * 125 + K5[JJ] * 42) / 14;
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
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;
if FIRST then
begin FIRST:= false; HL:= H; H:= MU * H; goto ACC
end;
FHM:= MU * H / HL + MU - MU1; HL:= H; H:= FHM * H;
ACC: MU1:= MU;
for JJ:= 1 step 1 until N do
begin Y[JJ]:= YL[JJ] + (ZL[JJ] * 56 + K0[JJ] * 7 +
K2[JJ] * 36 - K4[JJ] * 15) / 56 * HL;
Z[JJ]:= ZL[JJ] + ( - K0[JJ] * 63 + K1[JJ] * 189
- K2[JJ] * 36 - K3[JJ] * 112 + K4[JJ] * 50) / 28
end;
for J:= 1 step 1 until N do K5[J]:= FXYZJ * HL;
for JJ:= 1 step 1 until N do
begin Y[JJ]:= YL[JJ] + (ZL[JJ] * 336 + K0[JJ] * 35 +
K2[JJ] * 108 + K4[JJ] * 25) / 336 * HL;
Z[JJ]:= ZL[JJ] + (K0[JJ] * 35 + K2[JJ] * 162 +
K4[JJ] * 125 + K5[JJ] * 14) / 336
end;
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 RK2N;
eop