Simplifying code that displays a square rolling inside a circle

EDIT 1. Fixed kernel crashing. The problem was WhenEvent choking, because it was continually checking whether the current pivot (which lies exactly on the bounding circle) goes outside the bounding circle.

EDIT 2. Here is a port to GNU Octave: https://github.com/yawnoc/tumbling-polygon

---

System of ODEs in the complex plane

For a completely different take, observe that rigid rotation of the $n$th vertex $z_n$ around a pivot $p$ can be expressed as the differential equation

$$\frac{\mathrm{d}z_n}{\mathrm{d}t} = \mathrm{i} \,(z_n - p),$$

where $\mathrm{i} = \sqrt{-1}$, and we have scaled out the angular velocity. This is because the instantaneous velocity is (1) at right angles to the displacement vector $z_n - p$, and (2) proportional to the distance $|z_n - p|$. (Remember multiplication by $\mathrm{i}$ effects a 90-degree rotation.) If the vertex $z_n$ happens to coincide with the pivot location $p$, this nicely reduces to $\mathrm{d}z_n / \mathrm{d}t = 0$.

So we solve a system of such ODEs in $(z_1, z_2, z_3, z_4, p)$, where we handle the collision with the bounding circle $|z| = \sqrt{2}$ using WhenEvent; if $|z_n| > \sqrt{2}$, we set $p$ to $z_n$.

ClearAll["Global`*"];

vertexInitialPositions = Complex @@@ N @ {
  {-Sqrt[2]/2, -Sqrt[6]/2},
  {Sqrt[2]/2, -Sqrt[6]/2},
  {Sqrt[2]/2, Sqrt[2] - Sqrt[6]/2},
  {-Sqrt[2]/2, Sqrt[2] - Sqrt[6]/2},
  Nothing
};
vertexCount = Length[vertexInitialPositions];

vertexRotationEquations =
  Table[z[i]'[t] == I (z[i][t] - p[t]), {i, vertexCount}];

pivotFixtureEquation = p'[t] == 0;

vertexInitialConditions =
  Table[z[i][0] == vertexInitialPositions[[i]], {i, vertexCount}];

pivotInitialCondition = p[0] == First[vertexInitialPositions];

boundingRadius = Sqrt[2];
collisionHandling =
  Table[
    whenEventDummy[
      (* If this vertex isn't pivot AND goes outside bounding circle *)
      z[i][t] != p[t] && Abs @ z[i][t] > boundingRadius,
      (* Set the pivot to this vertex *)
      p[t] -> z[i][t]
    ]
    , {i, vertexCount}
  ] /. {whenEventDummy -> WhenEvent};
dependentVariables = Table[z[i], {i, vertexCount}] // Append[p];

tMax = Pi;
trajectories =
  NDSolveValue[
    Flatten @ {
      vertexRotationEquations, pivotFixtureEquation,
      vertexInitialConditions, pivotInitialCondition,
      collisionHandling
    },
    dependentVariables,
    {t, 0, tMax}
  ];

vertexRealTrajectories[t_] :=
  Most[trajectories][t] // Through // ReIm // Evaluate;

Manipulate[
  Show[
    Graphics @ Circle[{0, 0}, boundingRadius],
    Graphics @ Polygon @ vertexRealTrajectories[tCurrent],
    If[tCurrent == 0, {},
      ParametricPlot[
        vertexRealTrajectories[t]
        , {t, 0, tCurrent}
      ]
    ]
    , ImageSize -> 240
  ]
  , {tCurrent, 0, tMax}
]

The nice thing about this method is that it generalises to any polygon. Also the initial pivot need not start on the bounding circle (and in fact need not be any of the vertices).