Crinkled Rods by Will Webber

Use this function view crinkled rods

CylPlot[3,1,1.]        
(* The first input is the Symmetry.  
 The second input is the number of times to wrap around.
 The third input is the opacity.  
 Click and drag to change the point of view. *)

Use this Function for Stereo images

StereoCylPlot[3,1,11.]
(* The inputs are the same as above.  *)

In the hidden cells bellow are all the calculations

(* They were created for My Thesis at UW  *)

(*

Variables:
    r      radius of a point from the z-axis
    n      The symmetry of the figure. 
    s,t    Each torus is made of copies of the same scalene triangle
           having base 1 unit.  t is where the foot of the altitude is
           and s is the length of the altitude.
    qlist  This is a list of plus and minus ones which determines the 
           fold sequence.  In the routine Z3b we are finding the 
           solution of a quadratic equation. Thus we have a choice 
           of which solution.
    q      One element of qlist.
    xs     This stands for cross section.
    v      View point.

It should be noted that these procedures produce tori with antiprismatic
symmetry groups.*)

A[r1_,r2_,z2_,n_] := N[G[r1,r2,n]^2+z2^2]

B[z2_,t_] := N[-2*t*z2]

K[r1_,r2_,n_,t_] := N[t^2-G[r1,r2,n]^2]

G[r1_,r2_,n_] := N[r1 - r2*Cos[Pi/n]]

P[r1_,r2_,n_,q_] := N[If[(Re[r2]*Cos[Pi/n] > Re[r1]), -1*q, q]]

Z2[r1_,r2_,n_,s_,t_] := N[Sqrt[s^2 + t^2 - r1^2 - r2^2 + 2*r1*r2*Cos[Pi/n]]]

Z3b[r1_,r2_,z2_,n_,t_,p_] := 
         N[(-1*B[z2,t] + p*Sqrt[B[z2,t]^2 - 4*A[r1,r2,z2,n]*K[r1,r2,n,t]]) /
         (2*A[r1,r2,z2,n])]

R3b[r1_,r2_,z2_,z3_,n_,t_] := N[r1 + (t-z2*z3)/(r2*Cos[Pi/n] - r1)]

(* Fold2 returns a list of coordinates of points in a 
   double cross section around the bottom of the torus. *)

Fold2[n_,s_,t_,r_,qlist_] := 
        N[Module[{z2},
                 z2 = Z2[r,r,n,s,t] ;
                 Prepend[Prepend[Foldloop3[n,r,r,z2,0,t,qlist],
                                 {r,0}],
                         {r,z2}]
         ]      ]

(* Fold3 returns a list of coordinates of points in a 
   double cross section around the top of the torus. *)

Fold3[n_,s_,t_,r_,qlist_] := 
                        N[Module[{z2},
                                 z2 = Z2[r,r,n,s,t];
                                 Prepend[Prepend[Foldloop3[n,r,r,0,z2,t,qlist],
                                                 {r,z2}],
                                         {r,0}]
                                ]]

Folder3[n_,s_,t_,r_,qlist_] := 
            N[{Fold3[n,s,t,r,qlist],Fold2[n,s,t,r,-1*qlist]}]

(* This is the routine that actually calculates the cross section. *)
(* Notice that t is toggled each time.  *)

Foldloop3[n_,r1_,r2_,z1_,z2_,t_,qlist_] := 
         N[If[qlist == {},{},
                Module[{r3,z3},
                       z3 = z1+Z3b[r1,r2,(z2-z1),n,t,P[r1,r2,n,First[qlist]]];
                       r3 = R3b[r1,r2,(z2-z1),(z3-z1),n,t];
                       Prepend[Foldloop3[n,r2,r3,z2,z3,1-t,Rest[qlist]],{r3,z3}]
                      ]
             ] ] 

(* The following stuff is to find a cylindrical polyhedron.*)


Cylrad[n_,s_,q_] := 
     N[Sqrt[(16*s^2+q*Sqrt[256*s^4-64*(4*s^2-1)*Sin[Pi/n]^2])/32/Sin[Pi/n]^2]]

Cyl[n_,s_,q_] := N[Fold3[n,s,.5,Cylrad[n,s,q],{-1,1,1,-1}]]

Cylplot[n_,j_,s_,m_] := 
     TorusPlot3[n,j,Xsection[Folder3[n,s,.5,Cylrad[n,s,-1],{-1,1,1,-1,-1,1,1,-1,-1,1,1}]],{100,10,50},m]

CylPlot[n_,j_,m_] :=Cylplot[n,j,Sqrt[3]/2,m]

Cylplot4[n_,j_,s_,m_] := 
     TorusPlot4[n,j,Xsection[Folder3[n,s,.5,Cylrad[n,s,-1],{-1,1,1,-1,-1,1,1,-1}]],{100,10,50},m]

Xsection[list_] := Module[{tlist},
                          tlist = Join[Reverse[list[[1]]],Drop[list[[2]],2]] ;
                          Table[Table[tlist[[k]],{k,2+l,Length[tlist]+l,2}],
                                {l,-1,0}]
                         ]

StereoCylPlot[n_,j_,m_] := StereoCyl[n,j,Sqrt[3]/2,m]

StereoCyl[n_,j_,s_,m_] :=
    Module[{temp},
           temp = Cylplot[n,j,s,m] ;
           Show[Graphics[{Rectangle[{0,0},{1,1},temp],
                          Rectangle[{1,0},{2,1},Show[temp,ViewPoint -> 
                                                     RotateW[{100,10,50},5]]]}
                        ],
                AspectRatio -> Automatic]
          ]

Coord[l_,k_,n_,xs_]:=N[{xs[[(-1)^l/2+1.5]][[k]][[1]]*Cos[Pi*(l-1)/n],xs[[(-1)^l/2+1.5]][[k]][[1]]*Sin[Pi*(l-1)/n],xs[[(-1)^l/2+1.5]][[k]][[2]]}]

TorusPlot3[n_,j_,xs_,v_,m_]:=Show[Graphics3D[{Table[{GrayLevel[.95],Opacity[m],Polygon[{Coord[l,k,n,xs],Coord[l,k+1,n,xs],Coord[l+1,k,n,xs]}]},{l,1,2*n*j-1,2},{k,1,Length[xs[[1]]]-1}],Table[{GrayLevel[.70],Opacity[m],Polygon[{Coord[l,k,n,xs],Coord[l,k+1,n,xs],Coord[l+1,k+1,n,xs]}]},{l,2,2n*j,2},{k,1,Length[xs[[1]]]-1}],Table[{GrayLevel[.55],Opacity[m],Polygon[{Coord[l,k,n,xs],Coord[l,k+1,n,xs],Coord[l-1,k+1,n,xs]}]},{l,2,2*n*j,2},{k,1,Length[xs[[1]]]-1}],Table[{GrayLevel[.20],Opacity[m],Polygon[{Coord[l,k,n,xs],Coord[l,k+1,n,xs],Coord[l-1,k,n,xs]}]},{l,3,2*n*j+1,2},{k,1,Length[xs[[1]]]-1}]}],Boxed->False,AspectRatio->Automatic,(*Lighting->False,*) DisplayFunction->Identity,ViewPoint->v,SphericalRegion->True]

TorusPlot4[n_,j_,xs_,v_,m_]:=Show[Graphics3D[{Table[{GrayLevel[.95],Opacity[m],Polygon[{Coord[l,k,n,xs]+{0.5Cos[Pi*(l-1)/n], -0.5Sin[Pi*l],0},                                                                                                                            Coord[l,k+1,n,xs]+{0.5Cos[Pi*l],-0.5Sin[Pi*l],0},                                                                         Coord[l+1,k,n,xs]+{0.5Cos[Pi*l],-0.5Sin[Pi*l],0}}]},{l,1,2*n*j-1,2},{k,1,Length[xs[[1]]]-1}],Table[{GrayLevel[.90],Opacity[m],Polygon[{Coord[l,k,n,xs]-{0.5Cos[Pi*l], 0.5Sin[Pi*l],0},Coord[l,k+1,n,xs]+{-0.5Cos[Pi*l], 0.5Sin[Pi*l],0},Coord[l+1,k+1,n,xs]+{-0.5Cos[Pi*l], 0.5Sin[Pi*l],0}}]},{l,2,2n*j,2},{k,1,Length[xs[[1]]]-1}],Table[{GrayLevel[.85],Opacity[m],Polygon[{Coord[l,k,n,xs],Coord[l,k+1,n,xs],Coord[l-1,k+1,n,xs]}]},{l,2,2*n*j,2},{k,1,Length[xs[[1]]]-1}],Table[{GrayLevel[.80],Opacity[m],Polygon[{Coord[l,k,n,xs],Coord[l,k+1,n,xs],Coord[l-1,k,n,xs]}]},{l,3,2*n*j+1,2},{k,1,Length[xs[[1]]]-1}]}],Boxed->False,AspectRatio->Automatic,(*Lighting->False,*) DisplayFunction->Identity,ViewPoint->v,SphericalRegion->True]

(* Shader is a procedure that determines the GrayLevel of a 
   plane that contains points p1, p2 and p3 if the sun is at
   sun.  *)

Shader[p1_,p2_,p3_,sun_] :=
      N[Module[{nor},
               nor = CrossP[p2-p1,p3-p1] ;
               .6+nor.sun/3/Len[nor]/Len[sun]
       ]      ]

CrossP[u_,v_] := N[{u[[2]]*v[[3]]-u[[3]]*v[[2]],
                   u[[3]]*v[[1]]-u[[1]]*v[[3]],
                   u[[1]]*v[[2]]-u[[2]]*v[[1]]}]

Len[v_] := N[Sqrt[v[[1]]^2+v[[2]]^2+v[[3]]^2]]

(* Rotate rotates the vector v around the z-axis through an angle s.
   Notice that s is in degrees.   *)

RotateW[v_,s_] :=
               N[{{Cos[s*Pi/180],-Sin[s*Pi/180],0},
                  {Sin[s*Pi/180], Cos[s*Pi/180],0},
                  {     0       ,      0       ,1}}.v]

Sh[glist_] :=
     N[Show[Graphics[Table[Rectangle[{0,i},{2,i+1},glist[[i]]],
                           {i,1,Length[glist]}]],
            AspectRatio -> Automatic
      ]    ]

PSPdown[x_] := Display["!psfix | lpr -h -Plw340", x]

MyArcTan[x_,y_] := 
      N[If[x<0 && y<0,
           2*Pi+ArcTan[x,y],
           ArcTan[x,y]
       ]  ]