#AvekD 4D Hypercube (Now uses equations)

restart:with(plots):with(plottools):

proj4D3D := proc(p, l)
description "Stereographically Project a 4d point into 3d":
Matrix([[1/(l-p[4, 1]), 0, 0, 0], [0, 1/(l-p[4, 1]), 0, 0], [0, 0, 1/(l-p[4, 1]), 0]]).p;
end proc:

proj3D2D := proc(p)
description "Project a 3d point into 2d":

Matrix([[p[1, 1] + cos(Pi/4) * p[3, 1]], [p[2, 1] + sin(Pi/4) * p[3, 1]]]);
end proc:

proj4D2D := proc(p, l)
description "Stereographically Project a 4d point into 3d, then project that point into 2d":

[eval(proj3D2D(proj4D3D(p, l))[1, 1]), eval(proj3D2D(proj4D3D(p, l))[2, 1])]
end proc:

rXZYW := proc(p, theta)
description "Rotate a point by an angle in the xz plane, then rotate it by that same angle in the yw plane":

Matrix([[cos(theta), 0, -sin(theta), 0], [0, cos(theta), 0, -sin(theta)], [sin(theta), 0, cos(theta), 0], [0, sin(theta), 0, cos(theta)]]).p;

end proc:

ltocolmx := proc(l)


Matrix(4,1,l);

end proc:

lines := []:

for a from -1 by 2 to 1 do
 for b from -1 by 2 to 1 do
  for c from -1 by 2 to 1 do
   for d from -1 by 2 to 1 do
    for e from 1 to 4 do
     if ([a, b, c, d][e] = -1) then
      if (e = 1) then
       lines := [op(lines), [ltocolmx([a, b, c, d]), ltocolmx([1, b, c, d])]];
      elif (e = 2) then
       lines := [op(lines), [ltocolmx([a, b, c, d]), ltocolmx([a, 1, c, d])]];
      elif (e = 3) then
       lines := [op(lines), [ltocolmx([a, b, c, d]), ltocolmx([a, b, 1, d])]];
      else
       lines := [op(lines), [ltocolmx([a, b, c, d]), ltocolmx([a, b, c, 1])]];
      end if;
     end if;
    end do;
   end do;
  end do;
 end do;
end do;

frames := []:

plist := []:

redsquares := [
[ltocolmx([-1, -1, -1, -1]), ltocolmx([1, -1, -1, -1]), ltocolmx([1, 1, -1, -1]), ltocolmx([-1, 1, -1, -1])], #XY
[ltocolmx([-1, -1, 1, -1]),  ltocolmx([1, -1, 1, -1]),  ltocolmx([1, 1, 1, -1]),  ltocolmx([-1, 1, 1, -1])], #XY + Z
[ltocolmx([-1, -1, -1, -1]), ltocolmx([1, -1, -1, -1]), ltocolmx([1, -1, 1, -1]), ltocolmx([-1, -1, 1, -1])], #XZ
[ltocolmx([-1, 1, -1, -1]),  ltocolmx([1, 1, -1, -1]),  ltocolmx([1, 1, 1, -1]),  ltocolmx([-1, 1, 1, -1])], #XZ + Y
[ltocolmx([-1, -1, -1, -1]), ltocolmx([-1, -1, 1, -1]), ltocolmx([-1, 1, 1, -1]), ltocolmx([-1, 1, -1, -1])], #ZY
[ltocolmx([1, -1, -1, -1]),  ltocolmx([1, -1, 1, -1]),  ltocolmx([1, 1, 1, -1]),  ltocolmx([1, 1, -1, -1])] #ZY + X
]:

bluesquares := [
[ltocolmx([-1, -1, -1, 1]), ltocolmx([1, -1, -1, 1]), ltocolmx([1, 1, -1, 1]), ltocolmx([-1, 1, -1, 1])], #XY + W
[ltocolmx([-1, -1, 1, 1]),  ltocolmx([1, -1, 1, 1]),  ltocolmx([1, 1, 1, 1]),  ltocolmx([-1, 1, 1, 1])], #XY + Z + W
[ltocolmx([-1, -1, -1, 1]), ltocolmx([1, -1, -1, 1]), ltocolmx([1, -1, 1, 1]), ltocolmx([-1, -1, 1, 1])], #XZ + W
[ltocolmx([-1, 1, -1, 1]),  ltocolmx([1, 1, -1, 1]),  ltocolmx([1, 1, 1, 1]),  ltocolmx([-1, 1, 1, 1])], #XZ + Y + W
[ltocolmx([-1, -1, -1, 1]), ltocolmx([-1, -1, 1, 1]), ltocolmx([-1, 1, 1, 1]), ltocolmx([-1, 1, -1, 1])], #YZ + W
[ltocolmx([1, -1, -1, 1]),  ltocolmx([1, -1, 1, 1]),  ltocolmx([1, 1, 1, 1]),  ltocolmx([1, 1, -1, 1])] #YZ + X + W
]:

nframes := 200:

for i from 1 to nframes do
 for j in lines do
  plist := [plist[], plot([[proj4D2D(rXZYW(j[1], 2*(i-1)*Pi/nframes), 3), proj4D2D(rXZYW(j[2], 2*(i-1)*Pi/nframes), 3)]], color=blue, thickness=3)]:
 end do:
 for k in redsquares do
  plist := [plist[], polygon([proj4D2D(rXZYW(k[1], 2*(i-1)*Pi/nframes), 3), proj4D2D(rXZYW(k[2], 2*(i-1)*Pi/nframes), 3), proj4D2D(rXZYW(k[3], 2*(i-1)*Pi/nframes), 3), proj4D2D(rXZYW(k[4], 2*(i-1)*Pi/nframes), 3)], color=red, filled=true, transparency=0.9)]:
 end do:
 for m in bluesquares do
  plist := [plist[], polygon([proj4D2D(rXZYW(m[1], 2*(i-1)*Pi/nframes), 3), proj4D2D(rXZYW(m[2], 2*(i-1)*Pi/nframes), 3), proj4D2D(rXZYW(m[3], 2*(i-1)*Pi/nframes), 3), proj4D2D(rXZYW(m[4], 2*(i-1)*Pi/nframes), 3)], color=blue, filled=true, transparency=0.9)]:
 end do:
 frames := [frames[], display(plist)]:
 plist := []:
end do:

movie:=display(frames, insequence=true, scaling=constrained, axes=none):

movie;