#Crib: http://www.math.uni-wuppertal.de/~xsc/literatur/PXSCENGL.pdf # INT n = 4; # x, xd, y, yd # []STRING cmm = ("x ", "xd ", "y ", "yd "); MODE R = LONG LONG REAL; MODE RF = PROC(R)R; RF rsqrt = long long sqrt; RF rsin = long long sin, rcos = long long cos; MODE SV = [n]R; PROC hlf = (SV x)SV: (SV a; FOR j TO n DO a[j] := 0.5*x[j] OD; a); PROC prs = (STRING m, SV x)VOID: (print((m, newline)); FOR j TO n DO print((cmm[j], x[j], newline)) OD); OP * = (R s, SV v)SV: (SV a; FOR j TO n DO a[j] := s*v[j] OD; a); OP * = (REAL s, SV v)SV: (SV a; FOR j TO n DO a[j] := s*v[j] OD; a); OP + = (SV x, y)SV: (SV a; FOR j TO n DO a[j] := x[j] + y[j] OD; a); # This code is to solve the ordinary differential equation dy/dx = f(x, y) # # in particular a Newtonian circular orbit # PROC f = (R tau, SV y)SV: (R d2 = y[1]*y[1] + y[3]*y[3], d = rsqrt(d2), dm3 = 1.0/(d*d2); (SV a; FOR i BY 2 TO n DO a[i] := y[i+1]; a[i+1] := -y[i]*dm3 OD; a)); PROC step = (R dx, x, SV y)SV: ( SV k1 = dx*f(x,y); SV k2 = dx*f(x+dx/2, y + hlf(k1)); SV k3 = dx*f(x+dx/2, y + hlf(k2)); SV k4 = dx*f(x+dx, y+k3); y + 1.0/R(6)*(k1+2.0*k2 + 2.0*k3 + k4)); (SV x := (1, 0, 0, 1); PROC rep = (INT j)VOID: (print((j, newline)); prs("dvs", x); prs("psh", f(0, x))); INT c = 10000; rep(0); FOR k TO c DO x := step (1/c, 0, x); (k=1 OR k=c | rep(k)) OD; print((x[1] - rcos(1), newline, x[3] - rsin(1)))) # +0 dvs x +1.0000000000000000000000000000000000000000000000000000000e +0 xd +0.0000000000000000000000000000000000000000000000000000000e +0 y +0.0000000000000000000000000000000000000000000000000000000e +0 yd +1.0000000000000000000000000000000000000000000000000000000e +0 psh x +0.0000000000000000000000000000000000000000000000000000000e +0 xd -1.0000000000000000000000000000000000000000000000000000000e +0 y +1.0000000000000000000000000000000000000000000000000000000e +0 yd +0.0000000000000000000000000000000000000000000000000000000e +0 +1 dvs x +9.9999999500000000416666665885416677408854154968261735916e -1 xd -9.9999999833333333489583333138020825093587253392537595797e -5 y +9.9999999833333333333333332161458338216145832570393906911e -5 yd +9.9999999500000000416666665885416659830729221700032631048e -1 psh x -9.9999999833333333489583333138020825093587253392537595797e -5 xd -9.9999999500000000416666670312499977213541731296115441865e -1 y +9.9999999500000000416666665885416659830729221700032631048e -1 yd -9.9999999833333333333333336588541652777777796919080979424e -5 +10000 dvs x +5.4030230586813971729713847590873279222936617149931327057e -1 xd -8.4147098480789650778474645471989997252818746298909865834e -1 y +8.4147098480789650586543092183827677525190449567089409277e -1 yd +5.4030230586813971696690300006023594663226673250267114350e -1 psh x -8.4147098480789650778474645471989997252818746298909865834e -1 xd -5.4030230586813971846156577493233868054549194547538784052e -1 y +5.4030230586813971696690300006023594663226673250267114350e -1 yd -8.4147098480789650767891899213668427809323381818154931309e -1 -1.0379813153424381150294424911860895709611314300000000000e -19 -7.8707139979202222437065856512747697290381702000000000000e -19 #