# This does the 3 body 3-4-5 configuration. It reports ke, pe & te. This includes a variable dt. This has logic to compare two nearby paths thru phase space. The Runge Kutta scheme. Switching to d dimensions instead of 2. This code has the dubious ability to track masses in 4D thru an inverse square law field. The three lines below with "chnge" must modified together # PROC i2s = (INT j)STRING: (j<0 | "-" + i2s(-j) |: j < 10 | STRING("0123456789"[j]) | i2s(ENTIER (j/10)) + i2s(j%*10)); MODE FL = # LONG LONG chnge # REAL; # prefix LONG for long real types # PROC (FL)FL srt = # long long chnge # sqrt; # prefix long for long real types # PROC f2r = (FL x)REAL: # SHORTEN SHORTEN chnge # x; INT case = 1; # which problem? # INT d = 2, # d (number of dimensions) # n = (case|3, 2, 2); # the mass count. # PROC sq = (FL x)FL: x*x; MODE DV = STRUCT(FL x, xd); MODE PT = STRUCT([d]DV dv, FL m); # all about a point # MODE SS = STRUCT([n]PT pos, FL t, dt); # a system state inorder to compare states # PROC spin = (REF [] SS gurg, FL del, BOOL draw)VOID: ( print(("Beginning "+(case | "345 configuration", "Circle", "Parabola"), newline)); INT gx := 0; FILE out; INT nm := 1; PROC op = (INT j)VOID: open (out, "trail" + "ABCDEFGHIJK"[j], standback channel); (draw|op(1)); INT ccx := 0; REAL x := 0; PROC putc = (INT c)VOID: (draw| x := x*256 + c; ccx +:= 1; (ccx = 4 | putbin(out, ROUND(x>2147483647.0|x-4294967296.0|x)); ccx := 0; x := 0)); FL dt := .000001, t := 0; INT c := 0; # c was n # PROC clse = VOID: (draw | WHILE ccx > 0 DO putc(0) OD; close(out)); [n]PT z := (case | ((((del,0),(4,0)),3),(((3,0),(0,0)),4),(((0,0),(0,0)),5)), # 345 # ((((-1+del,0),(0,-.5)),1),(((1,0),(0,.5)),1)), # Circle # ((((0,srt(2)),(del,0)),1),(((0,-srt(2)),(0.5,0)),1)) # parabolas #); PROC ri = VOID: (FL tm := 0; [d]FL cg, m; [d,d]FL l; FOR j TO d DO cg[j] := m[j] := 0; FOR k TO d DO l[j,k] := 0 OD OD; FOR j TO n DO tm +:= m OF z[j]; FOR k TO d DO cg[k] +:= x OF (dv OF z[j])[k]; m[k] +:= xd OF (dv OF z[j])[k] OD; FOR k TO d DO FOR m TO d DO l[k,m] +:= m OF z[j] * ( x OF (dv OF z[j])[k] * xd OF (dv OF z[j])[m] - x OF (dv OF z[j])[m] * xd OF (dv OF z[j])[k]) OD OD OD; PROC pdex = ([]FL tp, FL f, STRING s) VOID: ( FOR j TO UPB tp DO print(tp[j]*f) OD; print(("<-", s, newline))); pdex(cg, 1/tm, "cg"); pdex(m, 1, "mass"); FOR j TO d DO pdex(l[j,], 1, "L") OD); PROC e = (REF PT j, k)FL:(FL s := 0; FOR m TO d DO s +:= sq (x OF (dv OF j)[m] - x OF (dv OF k)[m]) OD; - m OF j*m OF k/srt(s)); PROC ke = (REF PT j)FL: (FL s := 0; FOR m TO d DO s +:= sq(xd OF (dv OF j)[m]) OD; m OF j*s/2); PROC pr = (STRING mes)VOID: (FL ek:=0, ep:=0; FOR j TO n DO FOR k TO j-1 DO ep +:= e(z[j], z[k]) OD; ek +:= ke(z[j]) OD; print((mes+", ", "t=", t, " dt=", dt, newline, ek+ep)); print((" cy=", c, " pe=", ep, " ke=", ek, newline)); FOR k TO n DO FOR n TO d DO print(((dv OF z[k])[n])) OD; print(newline) OD; (case | SKIP, print(("disc = ", long cos(t/2) - (x OF (dv OF z[2])[1]), long sin(t/2) - (x OF (dv OF z[2])[2]), newline)), (FL x = (x OF (dv OF z[1])[1]), y = (x OF (dv OF z[1])[2]); print(("disc = ", x*x - y, " area=", x*(3 + x*x*4) - 3*srt(2)*t, newline))))); PROC step = VOID:(MODE KT = [n,d]DV; BOOL slow := 0<0, notfast := 0<0; PROC f = (REF KT y, BOOL dtc, REF KT nv)VOID: ( FOR n TO n DO FOR da TO d DO x OF nv[n,da] := xd OF y[n,da]; xd OF nv[n,da] := 0 OD OD; FOR j TO n DO FOR k TO j-1 DO [d]FL dx; FL s2 := 0; FOR d TO d DO s2 +:= sq(dx[d] := x OF y[j,d] - x OF y[k,d]) OD; FL ds = m OF z[j] * m OF z[k]/(s2 * srt(s2)); FOR dd TO d DO FL f = dx[dd]*ds; xd OF nv[j,dd] -:= f/m OF z[j]; xd OF nv[k,dd] +:= f/m OF z[k]; IF dtc THEN FL b := 0; FOR dd TO d DO b +:= sq(xd OF (dv OF z[j])[dd] - xd OF (dv OF z[k])[dd]) OD; b *:= sq(dt); (s2 < 600000.0 *b | slow := 0=0 | (s2 < 1500000.0 *b | notfast := 0=0)) FI OD OD OD); PROC lc = (FL sc, REF KT x, y, REF KT rt)VOID: (FOR j TO n DO FOR dd TO d DO rt[j,dd] := (sc*x OF x[j,dd] + x OF y[j,dd], sc*xd OF x[j,dd] + xd OF y[j,dd]) OD OD); KT k0; FOR j TO n DO FOR dd TO d DO k0[j,dd] := (dv OF z[j])[dd] OD OD; KT k1, k2, k3, k4, k5, kx, ky; f(k0, 0=0, k1); lc(dt/2, k1, k0, kx); f(kx, 0<0, k2); lc(dt/2, k2, k0, kx); f(kx, 0<0, k3); lc(dt, k3, k0, kx); f(kx, 0<0, k4); lc(2, k2, k1, ky); lc(2, k3, k4, kx); lc(1, ky, kx, k5); FOR j TO n DO FOR dd TO d DO x OF (dv OF z[j])[dd] +:= dt*x OF k5[j,dd]/6; xd OF (dv OF z[j])[dd] +:= dt*xd OF k5[j,dd]/6 OD OD; t +:= dt; (slow| dt *:= 0.5 |: NOT notfast | dt *:= 2)); PROC qt = (INT q, BOOL r)VOID: (clse; print(("qt ", q)); pr("fin"); ri; (r|stop)); IF UPB z /= n THEN qt(4, 0=0) FI; pr("init"); ri; MODE PX = STRUCT(INT x, y); INT nct := -1; FL tretic := 0; [n]PX v; INT ol := 200000; FL fst := 8, cpnt := 1; PROC pc = (FL y)INT: (REAL x = f2r(y); x<-2.5 | 0 |: x < 17.5 |ENTIER (100*x + 250)|999); PROC cx = (PX w, INT j)BOOL: (PX a = v[j]; v[j] := w; x OF a /= x OF w ORF y OF a /= y OF w); OP & = (INT x, y)INT: ABS (BIN x & BIN y); DO (c>nct | pr(">"); nct +:= 1000); (draw| PROC ph = (INT s)VOID: (putc(ENTIER (s/256)); putc(255 & s)); (t > fst | clse; nm +:= 1; op(nm); fst +:= 8); FOR k TO n DO (cx((pc(x OF (dv OF z[k])[1]), pc(x OF (dv OF z[k])[2])), k) | putc(k-1); ph(x OF v[k]); ph(y OF v[k]); (ol -:= 1; ol<0 |qt(2, 0=0))) OD; (t > tretic | FOR k TO n DO putc(3); ph(x OF v[k]); ph(y OF v[k]) OD; tretic+:=1)); ((gurg ISNT REF [] SS(NIL)) ANDF t > cpnt | PROC grab = SS: (SS s; FOR j TO n DO (pos OF s)[j] := z[j] OD; t OF s := t; dt OF s := dt; s); SS tmp = grab; dt := cpnt - t; step; gurg[gx +:= 1] := grab; cpnt +:= 1.0; FOR j TO n DO z[j] := (pos OF tmp)[j] OD; t := t OF tmp; dt := dt OF tmp; (gx = UPB gurg | GO TO end spin)); step; c +:= 1 OD EXIT end spin: SKIP); (INT last = 62; PROC mabs = (FL a, b)FL: (FL x = (a<0 | -a | a); (x