#Crib: http://www.math.uni-wuppertal.de/~xsc/literatur/PXSCENGL.pdf # INT n = 3; # alpha, x, y # []STRING cmm = ("alpha ", "x ", "y "); MODE R = LONG LONG REAL; MODE RF = PROC(R)R; RF rsin = long long sin, rcos = long long cos; RF rsqrt = long long sqrt; 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 beam bent so its ends meet # # 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 f = (R tau, SV y)SV: (-y[3], rcos(y[1]), rsin(y[1])); 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 := (0, 0, 1); PROC rep = (REAL j)VOID: (print((newline, j, newline)); prs("dvs", x); (0=1|prs("psh", f(0, x)))); rep(0); R tau := 0; REAL dt = 0.0001, dr=0.1; REAL tr := 0; WHILE x[3] > 0 DO x := step (dt, tau, x); tau := tau + dt; IF tau > tr-dt/2.0 THEN tr := tr + 0.1; rep(SHORTEN SHORTEN tau) FI OD)