Symbolic model of the locked triangle mesh

This Mathematica file accompanies the article “Locked rigid origami with multiple degrees of freedom” by Zachary Abel, Thomas C. Hull, and Tomohiro Tachi, which is to be published in the book (Origami^6): Proceedings of the Sixth International Meeting of Origami in Science, Mathemaics, and Education (to be published by the American Mathematical Society).

This Mathematica file contains notes for (1) the code for our symbolic model of the locked triangle mesh, using kinematic equations and allowing for the parameter angle θ to be changed, and (2) an attempt at a different configuration space---one that captures the self-intersection of the faces.

Here’s a drawing of the crease pattern:

And here is the vertex O rotated and with the creases labeled as in our model.  (So v2 = OB, v3 = OC, and v4 = OA.)

Graphics[{Line[{{0,0},Drop[v1,-1]}],Line[{{0,0},Drop[v2[0],-1]}],Line[{{0,0},Drop[v3[0.001,67°],-1]}],Line[{{0,0},Drop[v4p[0],-1]}],Line[{{0,0},Drop[v5,-1]}],
Text[Style["v1",Large],Drop[v1,-1]],Text[Style["v2",Large],Drop[v2[0],-1]],Text[Style["v3",Large],Drop[v3[0.001,67°],-1]],Text[Style["v4",Large],Drop[v4p[0],-1]],Text[Style["v5",Large],Drop[v5,-1]]}]

Vector stuff:  We’ll let a be our variable angle of the inner-most equilateral triangle.
v1 = (3/2, -sqrt(3)/2), 
v2 = (2,0) 
v3 = (k*Cos(180-a), k*Sin(180-a)) so that v3-v2 is perp to v2-v1
v4 is cut
v5 = (0,-sqrt(3))

Solve[(k*{Cos[180°-a],Sin[180°-a]}-{2,0}).({2,0}-{3/2,-Sqrt[3]/2})==0,k]

So k = -2/(Cos(a)-sqrt(3)Sin(a)).
Let’s do the dependent angles (p2=180°-a and p3 = a in Devlin’s notation):
(({{x}, {y}, {z}})) = (Subscript[R, v2])(phi1)(Subscript[R, v3])(phi2)(({{Cos(180°)}, {Sin(180°)}, {0}}))  
= (({{1, 0, 0}, {0, Cos(phi1), -Sin(phi1)}, {0, Sin(phi1), Cos(phi1)}}))(({{Cos(180-a), -Sin(180-a), 0}, {Sin(180-a), Cos(180-a), 0}, {0, 0, 1}}))(({{1, 0, 0}, {0, Cos(phi2), -Sin(phi2)}, {0, Sin(phi2), Cos(phi2)}}))(({{Cos(a)}, {Sin(a)}, {0}}))
Let’s see what that is:

RotX[t_]:={{1,0,0},{0,Cos[t],-Sin[t]},{0,Sin[t],Cos[t]}}
RotY[t_]:={{Cos[t],0,-Sin[t]},{0,1,0},{Sin[t],0,Cos[t]}}
RotZ[t_]:={{Cos[t],-Sin[t],0},{Sin[t],Cos[t],0},{0,0,1}}

RotX[phi1].RotZ[180°-a].RotX[phi2].{Cos[a],Sin[a],0}//Simplify

Now the independent angles, which I’m making equal and calling t:
(({{x}, {y}, {z}}))=(Subscript[R, v1])(t)(Subscript[R, v5])(t)(({{Cos[180]}, {Sin[180]}, {0}}))
= RotZ(-30)RotX(-t)RotZ(30)RotZ(-90)RotX(-t)RotZ(90)(({{-1}, {0}, {0}}))

RotZ[-30°].RotX[-t].RotZ[30°].RotZ[-90°].RotX[-t].RotZ[90°].{-1,0,0}//Simplify

So the x-coordinates of these two things will give us phi2:

Solve[(1/8) (1-6 Cos[t]-3 Cos[2 t])==-((Cos[a])^2)-Cos[phi2] ((Sin[a])^2),phi2]

ph2[t_,a_]:=ArcCos[-(1/8) (1+8 ((Cos[a])^2)-6 Cos[t]-3 Cos[2 t]) ((Csc[a])^2)]

Let’s try getting phi1 by doing the angle between vectors thing.

Ang[u_,v_]:=ArcCos[(u.v)/(Norm[u]*Norm[v])]

RotZ[180°-a].RotX[phi2].{Cos[a],Sin[a],0}//Simplify

ph1[t_,a_]:=Ang[{-Cos[a] (-1+Cos[ph2[t,a]]) Sin[a],Sin[a] Sin[ph2[t,a]]},{(1/2) (Sqrt[3]) (2+3 Cos[t]) ((Sin[(t/2)])^2),(3/2) Cos[t] Sin[t]}]

Plot[ph1[t,66°],{t,0,Pi}]

Awesome.  Let’s try it.  I wanna keep v1 and v5 fixed.

v1={3/2,-Sqrt[3]/2,0};
v2[t_]:=RotZ[-30°].RotX[t].RotZ[30°].{2,0,0}
v3[t_,a_]:=RotZ[-30°].RotX[t].RotZ[30°].RotX[-ph1[t,a]].(-(2/(Cos[a]-(Sqrt[3]) Sin[a])){Cos[180°-a],Sin[180°-a],0})
v4[t_,a_]:=RotZ[-30°].RotX[t].RotZ[30°].RotX[-ph1[t,a]].RotZ[180°-a].RotX[ph2[t,a]].RotZ[-(180°-a)].{-1,0,0}
v5={0,-Sqrt[3],0};
v4p[t_]:=RotZ[-90°].RotX[-t].RotZ[90°].{-1,0,0}

Manipulate[
Graphics3D[{
Polygon[{{0,0,0},v1,v2[t]}],
Polygon[{{0,0,0},v2[t],v3[t,a]}],
Polygon[{{0,0,0},v3[t,a],v4[t,a]}],
Polygon[{{0,0,0},v1,v5}],
Polygon[{{0,0,0},v4p[t],v5}],
Polygon/@{#+v5&/@(RotZ[120°].#&/@{{0,0,0},v2[t],v3[t,a]})},
Polygon/@{#+v5&/@(RotZ[120°].#&/@{{0,0,0},v3[t,a],v4[t,a]})},
Polygon/@{#+v5&/@(RotZ[120°].#&/@{{0,0,0},v4p[t],v5})},
Polygon/@{#+v1&/@(RotZ[240°].#&/@{{0,0,0},v2[t],v3[t,a]})},
Polygon/@{#+v1&/@(RotZ[240°].#&/@{{0,0,0},v3[t,a],v4[t,a]})}
},PlotRange->{{-3,3},{-3,2},{-2,2}},Axes->True],
{t,0.001,Pi},{a,60°, 70°}]

So!  Look at that!  You can clearly see where the paper starts to self-intersect, and that happens right about when the kinematic equations break down.  This corresponds to the discontinuity in the configuration space that Zach’s analysis predicted (when we were at Barbados 2013).  

Note that we only want the “a” angle to go up to 66.7 degrees or so.  And this crease pattern does NOT fold flat, so we shouldn’t take the “t” folding angle all the way to 180 degrees.

Below is a single vertex , to look at in isolation.

Manipulate[
Graphics3D[{
Polygon[{{0,0,0},v1,v2[t]}],
Polygon[{{0,0,0},v2[t],v3[t,65°]}],
Polygon[{{0,0,0},v3[t,65°],v4[t,65°]}],
Polygon[{{0,0,0},v1,v5}],
Polygon[{{0,0,0},v4p[t],v5}]
},PlotRange->{{-3,3},{-3,2},{-2,2}},Axes->True],
{t,0.001,Pi}]

Below is animation code that can be exported to Quicktime.

Animate[
Graphics3D[{
Polygon[{{0,0,0},v1,v2[t]}],
Polygon[{{0,0,0},v2[t],v3[t,70°]}],
Polygon[{{0,0,0},v3[t,70°],v4[t,70°]}],
Polygon[{{0,0,0},v1,v5}],
Polygon[{{0,0,0},v4p[t],v5}],
Polygon/@{#+v5&/@(RotZ[120°].#&/@{{0,0,0},v2[t],v3[t,70°]})},
Polygon/@{#+v5&/@(RotZ[120°].#&/@{{0,0,0},v3[t,70°],v4[t,70°]})},
Polygon/@{#+v5&/@(RotZ[120°].#&/@{{0,0,0},v4p[t],v5})},
Polygon/@{#+v1&/@(RotZ[240°].#&/@{{0,0,0},v2[t],v3[t,70°]})},
Polygon/@{#+v1&/@(RotZ[240°].#&/@{{0,0,0},v3[t,70°],v4[t,70°]})}
},PlotRange->{{-3,3},{-3,2},{-2,2}},Axes->True,ViewPoint->{0,0,1}],
{t,0.001,2.6}]

Export["triunit70",%,"QuickTime"]

Below is a drawing of the parameter space when the folding angles α, β, and γ are all equal.  I tried to fit a quadratic curve (in red dash) to the curve near the critical value.

Sam=RegionPlot[(-1/2 Cos[a °] Cos[a °]-3/2 Cos[a °]+Sin[a °] Sin[a °])/2≤Cos[180 °-2θ °],{θ,60,70},{a,0,180},FrameTicks->{Automatic,Range[0,360,45]},
FrameLabel->{"Parameter angle θ","Folding Angle α=β=γ"},PlotLegends->{"Rigid-foldable region"}]
Bill=ContourPlot[θ °-1.164==-.4(a °-2.1)^2,
{θ,60,70},{a,0,180},ContourStyle->{Red,Thick,Dashed},PlotLegends->{"Quadratic θ=-0.4(a-2.1)^2+1.164"}]

Show[Sam,Bill]

Self-intersection configuration space

Now I want to look at self-intersections.  My claim is that the locked triangle mesh will have pieces of the paper self-intersecting when the crease v4 (or v4p, which is easier to compute) makes an acute angle with the normal vector of the plane determined by creases v2 and v3. 
(Below is a pic that might help visualize this.)

Graphics[{Line[{{0,0},Drop[v1,-1]}],Line[{{0,0},Drop[v2[0],-1]}],Line[{{0,0},Drop[v3[0.001,67°],-1]}],Line[{{0,0},Drop[v4p[0],-1]}],Line[{{0,0},Drop[v5,-1]}],
Text[Style["v1",Large],Drop[v1,-1]],Text[Style["v2",Large],Drop[v2[0],-1]],Text[Style["v3",Large],Drop[v3[0.001,67°],-1]],Text[Style["v4",Large],Drop[v4p[0],-1]],Text[Style["v5",Large],Drop[v5,-1]]}]

When the kinematic equations fail, tho, the ph1 and ph2 functions can return complex values.  So we need to filter for that and then check the dot product:
Cross[v2,v3].v4
I made an auxiliary function, GEQ, to help with this:

GEQ[a_]:=If[Element[a,Reals],
a≥0,False]

GEQ[Cross[v2[1],v3[1,65°]].v4p[1]]

Then here’s the configuration space when θ = 67°:

RegionPlot3D[GEQ[Cross[v2[b],v3[b,67°]].v4p[a]]&&GEQ[Cross[v2[a],v3[a,67°]].v4p[c]]&&GEQ[Cross[v2[c],v3[c,67°]].v4p[b]],{a,0.001,Pi},{b,0.001,Pi},{c,0.001,Pi},PlotPoints->50]

So it seems that no self-intersection happens along the α = β = γ line when θ = 67°.

But now let’s try the θ = 66.7163° critical value:

RegionPlot3D[GEQ[Cross[v2[b],v3[b,66.7163°]].v4p[a]]&&GEQ[Cross[v2[a],v3[a,66.7163°]].v4p[c]]&&GEQ[Cross[v2[c],v3[c,66.7163°]].v4p[b]],{a,0.001,Pi},{b,0.001,Pi},{c,0.001,Pi},PlotPoints->100]

Oh boy!  There might a self-intersection happening along the α = β = γ line here!  Just barely?

But notice, the bridges are gone.  I don’t think they can be folded without self-intersections.