MODE V = [3]REAL; #Type for vector#
MODE TEN = [3,3]REAL; #Type for tensor#
PROC pv = (V x)VOID: #Print vector#
    (FOR j TO 3 DO print(fixed(x[j], 12, 5)) OD; print(newline));

# Print named tensor:#
PROC pnt = (STRING n, TEN x)VOID: (print((n, newline));
         FOR j TO 3 DO pv(x[j,]) OD);

#Convert anti-symmetric tensor to vector.#
PROC tv = (TEN t)V: (FOR i TO 3 DO FOR j TO 3 DO
   (t[i,j] + t[j,i] NE 0 | pnt("Not anti-symmetric!!", t)) OD OD;
   (t[2,3], t[3,1], t[1,2]));

#Convert vector to anti-symmetric tensor.#
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;

#Print named vector#
PROC pnv = (STRING n, V x)VOID: (print((n, newline)); pv(x));
first random(43); #Initialize random generator#

[20]V ma; #20 equal mass points of rigid body#

FOR j TO UPB ma DO ma[j] := (rx, rx, rx) OD; #Initial locations randomly.#
TEN omet; #Angular velocity: randommize it below:#
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); #Axial Vector form of Angular velocity#
pnt("Omega tensor", omet); pnv("Omega axial vec", omev);

PROC sum = (PROC (INT)REAL pr)REAL: (REAL sm := 0; #general summer#
  FOR j TO 3 DO sm +:= pr(j) OD; sm);

PROC ip = (V a, b)REAL: sum((INT j)REAL: a[j]*b[j]); #inner product#
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] := ip(omet[,i],ma[k][]) 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);
V lv = tv(momen); pnv("Vector L", lv);
V ld;
FOR i TO 3 DO ld[i] := ip(ti[i,], omev[]) OD;
pnv ("Derived angular momentum", ld)

# Yields:
Omega tensor
    +0.00000    +0.62349    +0.36719
    -0.62349    +0.00000    -0.21410
    -0.36719    +0.21410    +0.00000
Omega axial vec
    -0.21410    -0.36719    +0.62349
L =
    +0.00000    +8.11694    +3.06260
    -8.11694    +0.00000    -1.20996
    -3.06260    +1.20996    +0.00000
moments
    +6.95937    +1.64894    -2.11501
    +1.64894    +8.20448    -2.40957
    -2.11501    -2.40957    +6.43406
Inertia Tensor
   +14.63854    -1.64894    +2.11501
    -1.64894   +13.39343    +2.40957
    +2.11501    +2.40957   +15.16385
Vector L
    -1.20996    -3.06260    +8.11694
Derived angular momentum
    -1.20996    -3.06260    +8.11694
#