(*
 *  
 * CSPlotter.m - A plotter for Contact structures in R^3
 *
 *    Matias Dahl  2002-2005
 *
 *
 * Mathematica package for plotting contact structures
 * in R^3. (Tested with Mathematica v5.0.1.)
 *
 * For more information on the use of this package, see 
 * the project page at  
 *     www.math.hut.fi/~fdahl/CSPlotter/
 *
 * Please feel free to do whatever you want with this
 * code. If you think it is appropriate, you may mention
 * this package, but that is not required. For example, 
 * in a poster or on slides, this would be completely 
 * irrelevant information. 
 * 
 *)

<< Graphics`Shapes`
<< Calculus`VectorAnalysis`


(*
 *  Rotvect[v]
 *     
 *  Suppose v is a non-zero vector in R^3. Then 
 *  Rotvect[v] is a 3x3 matrix representing the
 *  rotational matrix (an orthogonal matrix with det = 1)
 *  such that R.(0, 0, 1) points in the direction of v.
 * 
 *)

Rotvect[v_] :=
  Module[
    {ex, ey, ez},
    If[v[[1]] == 0 && v[[2]] == 0,
      Return[{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}];
      ]; 
    ez = v/Sqrt[v.v];
    ey = 1/Sqrt[v[[1]]^2 + v[[2]]^2] {-v[[2]], v[[1]], 0};
    ex = Cross[ey, ez];
    Transpose[{ex, ey, ez}]  
]

(*
 *  contactelement[disp, covect, scale]
 *   
 *  Returns a rectange centered at 'disp', normal to 'covect' 
 *  and of size 'scale' as a Mathematica Graphics3D object.
 *
 *)
contactelement[disp_, covect_, scale_] :=
  Module[
    {P, b1, b2, b3, b4},
    b1 = {-1/2, -1/2, 0}; b2 = {-1/2, 1/2, 0}; b3 = {1/2, 1/2, 0}; 
    b4 = {1/2, -1/2, 0};
    P = Rotvect[covect];
    gs = Graphics3D[Polygon[scale{ P.b1, P.b2, P.b3, P.b4}], Boxed -> False, 
        PlotRange -> All];
    gs = TranslateShape[gs, disp]
    ]


(*
 *  coords[xmax, ymax, zmax]
 *   
 *  returns (as a Mathematica Graphics3D object) coordinate
 *  axes of lengths 'xmax', 'ymax', and 'zmax', respectively. 
 *)

coords[xmax_, ymax_, zmax_] :=
  Graphics3D[
    {
      Line[{{0, 0, 0}, {0, 0, zmax}}],
      Line[{{0, 0, 0}, {0, ymax, 0}}],
      Line[{{0, 0, 0}, {xmax, 0, 0}}],
      Text["x", {xmax + 0.25, 0, 0}],
      Text["y", {0, ymax + 0.25, 0}],
      Text["z", {0, 0, zmax + 0.25}]
    },
    PlotRange -> All]