MODE V = [3]REAL;
MODE TEN = [3,3]REAL;
PROC tv = (TEN t)V: (FOR i TO 3 DO FOR j TO 3 DO
   (t[i,j] + t[j,i] NE 0 | print(("Not asym", t, newline))) OD OD;
   (t[2,3], t[3,1], t[1,2]));
PROC tt = (V x)TEN: ((0, x[3], -x[2]), (-x[3], 0, x[1]), (x[2], -x[1], 0));
PROC rx = REAL: 2*next random-1;
PROC pv = (V x)VOID:
    (FOR j TO 3 DO print(fixed(x[j], 12, 5)) OD; print(newline));
PROC pnv = (STRING n, V  x)VOID: (print((n, newline)); pv(x));
PROC pnt = (STRING n, TEN x)VOID: (print((n, newline));
         FOR j TO 3 DO pv(x[j,]) OD);
PROC ip = (V a, b)REAL: (REAL s:=0; FOR j TO 3 DO s +:= a[j]*b[j] OD; s);
first random(43);
[20]V ma;
FOR j TO UPB ma DO ma[j] := (rx, rx, rx) OD;
PROC sum = (PROC (INT)REAL pr)REAL: (REAL sm := 0;
  FOR j TO 3 DO sm +:= pr(j) OD; sm);
#print(sum((INT k)REAL: ma[19][k]));#
TEN omet;
FOR i TO 3 DO FOR j TO 3 DO omet[i,j] := (i<j|rx|:i=j|0|-omet[j,i]) OD OD;
V omev = tv(omet);
pnt("Omega tensor", omet); pnv("Omega axial vec", omev);
PROC v = (V m)[]REAL: ([3]REAL s;
  FOR j TO 3 DO s[j] := sum((INT q)REAL: omet[j,q]*m[q]) OD; s);
  
TEN momen; #From first principles#
FOR i TO 3 DO FOR j TO 3 DO momen[i,j] := 0 OD OD;
FOR k TO UPB ma DO V v; FOR i TO 3 DO v[i] :=
   (0=0|ip(omet[,i],ma[k][])| sum((INT a)REAL: (omet[a,i]*ma[k][a]))) OD;
  FOR i TO 3 DO FOR j TO 3
    DO momen[i,j] +:= ma[k][i]*v[j] - ma[k][j]*v[i] OD OD OD;
pnt("L =", momen);

#Now calculate Inertia Tensor.#
#but first the second moments:#
TEN m2;
FOR i TO 3 DO FOR j TO 3 DO m2[i,j] := 0 OD OD;
FOR k TO UPB ma DO
  FOR i TO 3 DO FOR j TO 3 DO m2[i,j] +:= ma[k][i]*ma[k][j] OD OD OD;
pnt("moments", m2);
TEN ti;
FOR i TO 3 DO FOR j TO 3 DO ti[i,j] :=
  (i=j|sum((INT a)REAL: m2[a,a])|0) - m2[i,j] OD OD;
pnt("Inertia Tensor", ti); 
(#And now Angular momentum from TI:#
TEN lti;
FOR i TO 3 DO FOR j TO 3 DO
  lti[i,j] := sum((INT a)REAL: ti[i,a]*omet[a,j]) OD OD;
pnt("L", lti);
#This is neither antisymetric nor obviously related
to the previous calculation of L.#
(TEN x; FOR i TO 3 DO FOR j TO 3 DO x[i,j] := lti[i,j] - lti[j,i] OD OD;
   pnt("antisymmetrized version of L", x)));
(V ccl; #Conventionally computed L summed over points#
PRIO X = 5;
OP X = (V a, b)V:
  (a[2]*b[3] - a[3]*b[2], a[3]*b[1] - a[1]*b[3], a[1]*b[2] - a[2]*b[1]);
OP + = (V a, b)V: (V x; FOR j TO 3 DO x[j] := a[j]+b[j] OD; x);
FOR j TO 3 DO ccl[j] := 0 OD;
FOR j TO UPB ma DO ccl := ccl + (ma[j] X (omev X ma[j])) OD;
pnv("Conventional L summed over points", ccl));
(V cxl; #Convention I omega calculation#
FOR i TO 3 DO cxl[i] := ip(ti[i,], omev[]) OD;
pnv("I omega version of L", cxl));
([3,3,3,3]REAL it2;
FOR i TO 3 DO FOR v TO 3 DO FOR j TO 3 DO FOR k TO 3 DO
  it2[i,v,j,k] := (v=k|1|0)*m2[j,i] + (j=i|1|0)*m2(k,v) OD OD OD OD;
(TEN l4; 
FOR j TO 3 DO FOR k TO 3 DO l4[j,k] :=
(0<0|# This clause awaits MK5! #
   sum((INT i)REAL: sum((INT v)REAL: it2[i,v,j,k]*omet[i,v]))
  |sum((INT i)REAL: ip(it2[i,,j,k],omet[i,])))
OD OD;
pnt("my L", l4); pnv("and axial vec thereof", tv(l4)));
FOR i TO 3 DO FOR j TO 3 DO pnt((j=1|"---"|""), it2[i,,j,]) OD OD)