(************************************************************************)
(*                                                                      *)
(*  Braids.m                                                            *)
(*  Draws braids                                                        *)
(*                                                                      *)
(*  version 2.0                                                         *)
(*  Joe Christy               10/25/91                                  *)
(*                                                                      *)
(************************************************************************)

(*
 * REVISION HISTORY
 *
 * version 1.0                4/15/90
 * version 1.1                4/17/90
 *    added Ribbons, made open braids and plats run from top to bottom.
 * version 1.2                5/3/90
 *    now handles framings of closed braids.
 * version 1.3                5/6/90
 *    through the magic of GrayLevel, crossings now visible in
 *    Curve Style.
 * version 1.4                2/5/91
 *    make default rendering Quick and dirty.
 * version 1.5                2/27/91
 *    make radius of closed braid depend on the number of crossings,
 *    so as to make crossings better looking.
 * version 2.0               10/25/91
 *	various Mathematica 2.0 updates
 *	JPC
 * version 2.1               11/4/91
 *  My Cross[], using Det[] broke with M--2.2, call standard
 *  LinearAlgebra package CrossProduct instead
 *  JPC
 *
 *)

(*
 *   TO DO
 *
 * Faster
 *
 *)


BeginPackage["Braids`",
	"SymmetricGroup`",
	"BraidGroup`"]

(* Messages *)

OpenBraid::usage =
	"OpenBraid[b, (options)] constructs a Graphic3D of b as an open braid.
	It takes a single option, Style (Automatic) which tells it whether to draw
	quick and dirty curves, nice curves, ribbons, or tubes."

Options[OpenBraid] = {Style -> Automatic}

ClosedBraid::usage =
	"ClosedBraid[b, (options)] constructs a Graphic3D of b as a closed braid.
	It takes the following options:\n
	Style (Automatic) which tells it whether to draw quick and dirty curves,
	nice curves, ribbons, or tubes.\n
	Canon (1) which tells it what the convention for the canonical framing is.\n
	Framing (Automatic) which twists the framing by some number of full twists relative
	to the canonical framing"
	
Options[ClosedBraid] = {Style -> Automatic, Canon->1, Framing->Automatic}

Plat::usage =
	"Plat[b, (options)] constructs a Graphic3D of b as a plat.
	It takes a single option, Style (Automatic) which tells it whether to draw
	quick and dirty curves, nice curves, ribbons, or tubes."

Options[Plat] = {Style -> Automatic}

Components::usage =
	"Components[b] returns the number of components of the closed braid b."
	
AlexPoly::usage =
	"AlexPoly[b] returns the Multivariable Alexander Polynomial in 
	variables T[1], T[2], ... , T[n] which represent the different components of
	the closed braid b."
	
T::usage = "The variables of the Alexander Polynomial."

Style::usage =
	"Style is an option to the braid drawing functions.\n
	The default is to draw quick and dirty curves, which is fast, but ugly.\n
	Style -> Curve draws smooth looking curves, which are slow, but cute.\n
	The two curve styles hide the framing.\n
	Style -> Ribbon draws ribbons, which are slower, but pretty.\n
	Style -> Tube draws tubes, which are very slow, but beautiful."

Curve::usage =
	"Style -> Curve draws smooth looking curves, which are slow, but cute."

Ribbon::usage =
	"Style -> Ribbon draws ribbons, which are pretty slow, but also pretty."

Tube::usage =
	"Style -> Tube draws tubes, which are very slow, but beautiful."

Canon::usage =
	"Canon (1) is an option to ClosedBraid which specifies the canonical framing.\n
	Canon -> 0 specifies the physicist's framing, determined by the normal to
	the projection plane.\n
	Canon -> 1 specifies the topologist's framing, determined by the boundary
	of the Seifert surface."
	
Framing::usage =
	"Framing (Automatic) is an option to ClosedBraid which specifies a framing
	by giving a list of the numbers of full twists around each component
	relative to the canonical framing."
	
Style::bogus =
	"`` is not a valid Style for ``!"

SetOptions[Graphics3D, Boxed -> False]

Begin["`Private`"]    (* begin the private context *)

(******************************************************************)
(* The functions below, based on Henry Cejtin's Tube package,     *)
(* take care of drawing tubular neighborhoods.                    *)
(******************************************************************)

(*
 *
 * Given a (symbolic) 3-space curve, a radius and a range over which to 
 * vary an independent variable, return a 3-D graphics object which 
 * represents a tubular neighborhood of the curve of the given radius.
 * Note, the radius may be depend on the independent variable.
 * the curve.
 * Note, the curve must be differentiable (by Mathematica) and the 
 * derivative vector must not vanish.
 * parity, either 1 or -1, is a kludge to adjust the framing to the 
 * physicists vertical framing.
 * Several options are understood by tblrNbhd:
 *	CurvePoints	the number of points to sample along the curve.
 *	TubePoints	the number of points to sample around the disks
 *				forming the tube.
 *	Twist		an angle by which to twist the mesh on the surface
 *				(defaulting to 0, or no twisting).
 *  Torsion     an angle by which to twist the mesh along the curve
 *				(defaulting to 0, or no twisting).
 *	Closed		if True, (the default) then the curve is assumed to
 *				be closed (i.e., end where it starts).
 *
 *)

Options[tblrNbhd] = {
		CurvePoints -> 50,
		TubePoints -> 8,
		Twist -> 0,
		Torsion -> 0,
		Closed -> True
	}

tblrNbhd[curve_, radius_, {var_, low_, high_}, parity_, opts___] :=
	Block[{
		CurvePoints,
		TubePoints,
		Twist,
		Torsion,
		Closed,
		angles,
		cpoints,
		vpoints,
		tpoints,
		tube,
		mesh,
		theta,
		delta,
		i,
		j,
		res,
		rot,
		rotm
		},
		{CurvePoints, TubePoints, Twist, Torsion, Closed} =
			{CurvePoints, TubePoints, Twist, Torsion, Closed}
				/. Options[{opts}]
				/. Options[tblrNbhd];
		cpoints = Table[	(* points along curve *)
				curve,
				{var, low, high, (high-low)/CurvePoints}
			] //N;
		If[Closed,	(* match ends of curve *)
			cpoints[[-1]] = cpoints[[1]]
		];
		tpoints = Table[	(* tangent vectors to curve *)
				Release[
					D[curve, var]
				],
				{var, low, high, (high-low)/CurvePoints}
			] //N;
		vpoints = tpoints = Map[	(* unit tangent vectors to curve *)
				Unitize,
				tpoints
			] // N;
		rot = Angle[	(* between initial and final tangents *)
					vpoints[[1]],
					vpoints[[-1]]
				] //N;
		rotm = {{Cos[rot],Sin[rot],0},
			{-Sin[rot],Cos[rot],0},
			{0,0,1}};
			(* rotation in x-z plane by angle between initial and final tangents *)
		tpoints = CTrans[	(* translates of initial normal *)
				tpoints,
				Cross[Perp[First[tpoints]],First[tpoints]]
			] //N;
		tpoints *= Table[	(* scaled normal field *)
				N[radius],
				{var, low, high, (high-low)/CurvePoints}
			];
		angles = Table[	(* constant twist field *)
				N[Twist],
				{var, low, high, (high-low)/CurvePoints}
			];
		If[Closed,
			delta = Angle[	(* angle between initial and final normals *)
				tpoints[[1, 1]],
				tpoints[[-1, 1]]
			] //N,
			delta = parity * Angle[	(* angle between initial and final normals after *)
				tpoints[[1, 1]],	(* turning of tangent is taken into account *)
				rotm . tpoints[[-1, 1]]
			] //N
		];
		If[delta != 0,	(* compensate for torsion of curve *)
			angles -= Table[
				theta,
				{theta, 0, delta, delta/CurvePoints}
			]
		];
		angles += Table[	(* add external torsion *)
				N[Torsion*dt],
				{dt, 0, 1, 1/CurvePoints}
			];
		tpoints = Map[	(* twist normal field *)
				Twist2,
				Transpose[{
					angles,
					tpoints
				}]
			];
		If[Closed,	(* force normals to match at ends *)
			tpoints[[-1]] = tpoints[[1]]
		];
		tube = Table[
				{Cos[var], Sin[var]},
				{var, 0, 2 Pi, 2 Pi / TubePoints}
			] //N;
		mesh = Re[Table[	(* Somehow these vectors can get to be complex! *)
				cpoints[[i]] +
					tube[[j]] . tpoints[[i]],
				{i, Length[cpoints]},
				{j, Length[tube]}
			]];
		res = ToGraphics[mesh];
		If[! Closed,
			res = Graphics3D[{
				Polygon[First[mesh]],
				res[[1]],
				Polygon[Last[mesh]]
			}]
		];
		res
	]

(*
 * Given a list of unit 3-dimensional vectors, and another vector,
 * return a list of pairs of vectors such that corresponding elements
 * of the original list and the result list form an orthonormal,
 * right-handed coordinate system.
 * This is done by assuming that the list represents tangents to a curve,
 * and parallely translating the orthonormal component of the given vector
 * along the curve.
 *)
CTrans[tlist_, vec_] :=
	Module[{
		v = vec
		},
		Map[
			{
				v = Unitize[v - (v . #) #],
				Cross[#, v]
			}&,
			tlist
		]
	]

(*
 * Given a 2-dimensional mesh of points on a surface, produce a
 * Graphics3D object which contains the polygons forming the
 * surface.
 *)
ToGraphics[mesh_] :=
	Module[{
		i,
		j
		},
		Graphics3D[
			Table[
				Polygon[{
					mesh[[i, j]],
					mesh[[i, j+1]],
					mesh[[i+1, j+1]],
					mesh[[i+1, j]]
				}],
				{i, Length[mesh] - 1},
				{j, Length[mesh[[1]]] - 1}
			]
		]
	]

(*
 * Given a non-zero vector, return any vector not parallel to it.
 * This is done by returning the coordinate vector for the coordinate
 * with minimal absolute value in the original vector.
 *)
Perp[u_] :=
	Module[{
		coord = 1,
		value = Abs[First[u]],
		i
		},
		Do[
			If[Abs[u[[i]]] <= value,
				value = Abs[u[[i]]];
				coord = i
			],
			{i, 2, Length[u]}
		];
		Table[
			If[i == coord,
				1,
				0
			],
			{i, Length[u]}
		]
	]

(*
 * Given an angle and two orthogonal vectors, return the result of rotating
 * the vectors through the angle.
 *)
Twist2[{angle_, vecs_}] :=
	Module[{
		cos = Cos[angle]//N,
		sin = Sin[angle]//N
		},
		{{cos,	sin},
		 {-sin,	cos}} . vecs
	];

(*
 * Given a (non-zero) vector, return a vector that points in the same
 * direction but is of unit length.
 *)
Unitize[vec_] := vec / Sqrt[vec . vec]


(*
 * Given two (non-zero) vectors, return the angle between them.
 *)
Angle[u_, v_] := ArcCos[Unitize[u] . Unitize[v]]

(**************************)
(* Some constants we need *)
(**************************)

twopi = N[2*Pi];

(***********************************************)
(* These functions take care of drawing curves *)
(* with the GrayLevel proportional to height.  *)
(***********************************************)

gly[{{x1_,y1_,z1_},{x2_,y2_,z2_}}]:=
	GrayLevel[2 (y1+y2)/3 + 0.4]

glz[{{x1_,y1_,z1_},{x2_,y2_,z2_}}]:=
	GrayLevel[0.4 - 2 (z1+z2)/3]

(*********************************************************************)
(* The functions dc* handle the construction of a single crossing.   *)
(* They are called repeatedly by the top level functions.            *)
(*********************************************************************)

(***************************)
(* open braided Q&D curves *)
(***************************)
(*
 * dcOQ[m, sn, cn]
 * draws the sn-th strand passing through the elementary
 * braid m at level cn. Open Curve.
 *)

dcoq[m_,sn_,cn_]:=
	Table[{m + 1/2 - Cos[Pi Sin[t]^2]/2,
			0.3 Sin[Pi Sin[t]^2],
			1 - 2 t/Pi -cn},
		{t, 0, Pi/2, Pi/4}] /; sn == m

dcoq[m_,sn_,cn_]:=
	Table[{m + 1/2 - Cos[Pi Sin[t]^2 + Pi]/2,
			0.3 Sin[Pi Sin[t]^2 + Pi],
			1 - 2 t/Pi -cn},
		{t, 0, Pi/2, Pi/4}] /; sn == m + 1

dcoq[m_,sn_,cn_]:=
	Table[{-m + 1/2 - Cos[Pi Sin[t]^2]/2,
			-0.3 Sin[Pi Sin[t]^2],
			1 - 2 t/Pi -cn},
		{t, 0, Pi/2, Pi/4}] /; sn == -m

dcoq[m_,sn_,cn_]:=
	Table[{-m + 1/2 - Cos[Pi Sin[t]^2 + Pi]/2,
			-0.3 Sin[Pi Sin[t]^2 + Pi],
			1 - 2 t/Pi -cn},
		{t, 0, Pi/2, Pi/4}] /; sn == -m + 1

dcoq[m_,sn_,cn_]:=
	Table[{sn, 0, 1 - 2 t/Pi -cn}, {t, 0, Pi/2, Pi/2}]

dcOQ[m_,sn_,cn_]:=
	Module[{pq3=
			Transpose[{
				Drop[dcoq[m,sn,cn],-1],
				Drop[RotateLeft[dcoq[m,sn,cn]],-1]
			}]
		},
		Graphics3D[Flatten[Transpose[{Map[gly,pq3],Map[Line,pq3]}],1]]
	]

(***********************)
(* open braided curves *)
(***********************)
(*
 * dcOC[m, sn, cn]
 * draws the sn-th strand passing through the elementary
 * braid m at level cn. Open Curve.
 *)

dco[m_,sn_,cn_]:=
	Table[{m + 1/2 - Cos[Pi Sin[t]^2]/2,
			0.3 Sin[Pi Sin[t]^2],
			1 - 2 t/Pi -cn},
		{t, 0, Pi/2, Pi/40}] /; sn == m

dco[m_,sn_,cn_]:=
	Table[{m + 1/2 - Cos[Pi Sin[t]^2 + Pi]/2,
			0.3 Sin[Pi Sin[t]^2 + Pi],
			1 - 2 t/Pi -cn},
		{t, 0, Pi/2, Pi/40}] /; sn == m + 1

dco[m_,sn_,cn_]:=
	Table[{-m + 1/2 - Cos[Pi Sin[t]^2]/2,
			-0.3 Sin[Pi Sin[t]^2],
			1 - 2 t/Pi -cn},
		{t, 0, Pi/2, Pi/40}] /; sn == -m

dco[m_,sn_,cn_]:=
	Table[{-m + 1/2 - Cos[Pi Sin[t]^2 + Pi]/2,
			-0.3 Sin[Pi Sin[t]^2 + Pi],
			1 - 2 t/Pi -cn},
		{t, 0, Pi/2, Pi/40}] /; sn == -m + 1

dco[m_,sn_,cn_]:=
	Table[{sn, 0, 1 - 2 t/Pi -cn}, {t, 0, Pi/2, Pi/2}]

dcOC[m_,sn_,cn_]:=
	Module[{pq3=
			Transpose[{
				Drop[dco[m,sn,cn],-1],
				Drop[RotateLeft[dco[m,sn,cn]],-1]
			}]
		},
		Graphics3D[Flatten[Transpose[{Map[gly,pq3],Map[Line,pq3]}],1]]
	]

(*****************************)
(* closed braided Q&D curves *)
(*****************************)
(*
 * dcCQ[m, sn, cn, c]
 * draws the sn-th strand passing through the 
 * elementary braid m at level cn out of c total crossings. Closed Curve.
 *)

dccq[m_,sn_,cn_,c_]:=
	Table[{(c/6 + m + 1/2 - Cos[Pi Sin[t]^2]/2)	Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + m + 1/2 - Cos[Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			0.3 Sin[Pi Sin[t]^2]},
		{t, 0, Pi/2, Pi/8}] /; sn == m

dccq[m_,sn_,cn_,c_]:=
	Table[{(c/6 + m + 1/2 - Cos[Pi + Pi Sin[t]^2]/2) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + m + 1/2 - Cos[Pi + Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			0.3 Sin[Pi + Pi Sin[t]^2]},
		{t, 0, Pi/2, Pi/8}] /; sn == m + 1

dccq[m_,sn_,cn_,c_]:=
	Table[{(c/6 + -m + 1/2 -Cos[Pi Sin[t]^2]/2)	Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + -m + 1/2 -Cos[Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			-0.3 Sin[Pi Sin[t]^2]},
		{t, 0, Pi/2, Pi/8}] /; sn == -m

dccq[m_,sn_,cn_,c_]:=
	Table[{(c/6 + -m + 1/2 -Cos[Pi + Pi Sin[t]^2]/2) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + -m + 1/2 -Cos[Pi + Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			-0.3 Sin[Pi + Pi Sin[t]^2]},
		{t, 0, Pi/2, Pi/8}] /; sn == -m + 1

dccq[m_,sn_,cn_,c_]:=
	Table[{(c/6 + sn) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + sn) Sin[2(cn - 1)Pi/c + 4 t/c],
			0},
		{t, 0, Pi/2, Pi/4}]

dcCQ[m_,sn_,cn_,c_]:=
	Module[{pq3=
			Transpose[{
				Drop[dccq[m,sn,cn,c],-1],
				Drop[RotateLeft[dccq[m,sn,cn,c]],-1]
			}]
		},
		Graphics3D[Flatten[Transpose[{Map[glz,pq3],Map[Line,pq3]}],1]]
	]

(*************************)
(* closed braided curves *)
(*************************)
(*
 * dcCC[m, sn, cn, c]
 * draws the sn-th strand passing through the 
 * elementary braid m at level cn out of c total crossings. Closed Curve.
 *)

dcc[m_,sn_,cn_,c_]:=
	Table[{(c/6 + m + 1/2 - Cos[Pi Sin[t]^2]/2)	Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + m + 1/2 - Cos[Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			0.3 Sin[Pi Sin[t]^2]},
		{t, 0, Pi/2, Pi/40}] /; sn == m

dcc[m_,sn_,cn_,c_]:=
	Table[{(c/6 + m + 1/2 - Cos[Pi + Pi Sin[t]^2]/2) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + m + 1/2 - Cos[Pi + Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			0.3 Sin[Pi + Pi Sin[t]^2]},
		{t, 0, Pi/2, Pi/40}] /; sn == m + 1

dcc[m_,sn_,cn_,c_]:=
	Table[{(c/6 + -m + 1/2 -Cos[Pi Sin[t]^2]/2) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + -m + 1/2 -Cos[Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			-0.3 Sin[Pi Sin[t]^2]},
		{t, 0, Pi/2, Pi/40}] /; sn == -m

dcc[m_,sn_,cn_,c_]:=
	Table[{(c/6 + -m + 1/2 -Cos[Pi + Pi Sin[t]^2]/2) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + -m + 1/2 -Cos[Pi + Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			-0.3 Sin[Pi + Pi Sin[t]^2]},
		{t, 0, Pi/2, Pi/40}] /; sn == -m + 1

dcc[m_,sn_,cn_,c_]:=
	Table[{(c/6 + sn) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + sn) Sin[2(cn - 1)Pi/c + 4 t/c],
			0},
		{t, 0, Pi/2, Pi/20}]

dcCC[m_,sn_,cn_,c_]:=
	Module[{pq3=
			Transpose[{
				Drop[dcc[m,sn,cn,c],-1],
				Drop[RotateLeft[dcc[m,sn,cn,c]],-1]
			}]
		},
		Graphics3D[Flatten[Transpose[{Map[glz,pq3],Map[Line,pq3]}],1]]
	]

(**********************)
(* open braided tubes *)
(**********************)
(*
 * dcOTN[m, sn, cn, tpn]
 * draws the sn-th strand passing through the 
 * elementary braid m at level cn.
 * Open Tubular Neighborhood with tpn points meridionally.
 *)

dcOTN[m_,sn_,cn_,tpn_]:=
	tblrNbhd[{m + 1/2 - Cos[Pi Sin[t]^2]/2,
			0.3 Sin[Pi Sin[t]^2],
			1 - 2 t/Pi -cn},
		0.2,
		{t, 0, Pi/2},
		1,
		CurvePoints->15, TubePoints->tpn, Closed->False] /; sn == m

dcOTN[m_,sn_,cn_,tpn_]:=
	tblrNbhd[{m + 1/2 - Cos[Pi Sin[t]^2 + Pi]/2,
			0.3 Sin[Pi Sin[t]^2 + Pi],
			1 - 2 t/Pi -cn},
		0.2,
		{t, 0, Pi/2},
		1,
		CurvePoints->15, TubePoints->tpn, Closed->False] /; sn == m + 1

dcOTN[m_,sn_,cn_,tpn_]:=
	tblrNbhd[{-m + 1/2 - Cos[Pi Sin[t]^2]/2,
			-0.3 Sin[Pi Sin[t]^2],
			1 - 2 t/Pi -cn},
		0.2,
		{t, 0, Pi/2},
		-1,
		CurvePoints->15, TubePoints->tpn, Closed->False] /; sn == -m

dcOTN[m_,sn_,cn_,tpn_]:=
	tblrNbhd[{-m + 1/2 - Cos[Pi Sin[t]^2 + Pi]/2,
			-0.3 Sin[Pi Sin[t]^2 + Pi],
			1 - 2 t/Pi -cn},
		0.2,
		{t, 0, Pi/2},
		-1,
		CurvePoints->15, TubePoints->tpn, Closed->False] /; sn == -m + 1

dcOTN[m_,sn_,cn_,tpn_]:=
	tblrNbhd[{sn, 0, 1 - 2 t/Pi -cn}, 0.2, {t, 0, Pi/2}, 1,
		CurvePoints->2, TubePoints->tpn, Closed->False]

(************************)
(* closed braided tubes *)
(************************)
(*
 * dcCTN[m, sn, cn, c, tpn, tw, skew, screw]
 * draws the sn-th strand passing through the 
 * elementary braid m at level cn out of c total crossings.
 * Closed Tubular Neighborhood with tpn points meridionally.
 * tw is an amount by which to twist the framing in a crossing.
 * skew is the amount to rotate the framing,
 * and screw is the amount by which to twist it.
 *)

dcCTN[m_,sn_,cn_,c_,tpn_,tw_,skew_,screw_]:=
	tblrNbhd[{(c/6 + m + 1/2 -Cos[Pi Sin[t]^2]/2) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + m + 1/2 -Cos[Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			0.3 Sin[Pi Sin[t]^2]},
		0.2,
		{t, 0, Pi/2},
		1,
		CurvePoints->15, TubePoints->tpn, Closed->False,
		Twist->N[twopi*skew], Torsion->N[twopi*(screw-tw)]] /; sn == m

dcCTN[m_,sn_,cn_,c_,tpn_,tw_,skew_,screw_]:=
	tblrNbhd[{(c/6 + m + 1/2 -Cos[Pi + Pi Sin[t]^2]/2) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + m + 1/2 -Cos[Pi + Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			0.3 Sin[Pi + Pi Sin[t]^2]},
		0.2,
		{t, 0, Pi/2},
		1,
		CurvePoints->15, TubePoints->tpn, Closed->False,
		Twist->N[twopi*skew], Torsion->N[twopi*(screw-tw)]] /; sn == m + 1

dcCTN[m_,sn_,cn_,c_,tpn_,tw_,skew_,screw_]:=
	tblrNbhd[{(c/6 + -m + 1/2 -Cos[Pi Sin[t]^2]/2) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + -m + 1/2 -Cos[Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			-0.3 Sin[Pi Sin[t]^2]},
		0.2,
		{t, 0, Pi/2},
		-1,
		CurvePoints->15, TubePoints->tpn, Closed->False,
		Twist->N[twopi*skew], Torsion->N[twopi*(screw+tw)]] /; sn == -m

dcCTN[m_,sn_,cn_,c_,tpn_,tw_,skew_,screw_]:=
	tblrNbhd[{(c/6 + -m + 1/2 -Cos[Pi + Pi Sin[t]^2]/2) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + -m + 1/2 -Cos[Pi + Pi Sin[t]^2]/2) Sin[2(cn - 1)Pi/c + 4 t/c],
			-0.3 Sin[Pi + Pi Sin[t]^2]},
		0.2,
		{t, 0, Pi/2},
		-1,
		CurvePoints->15, TubePoints->tpn, Closed->False,
		Twist->N[twopi*skew], Torsion->N[twopi*(screw+tw)]] /; sn == -m + 1

dcCTN[m_,sn_,cn_,c_,tpn_,tw_,skew_,screw_]:=
	tblrNbhd[{(c/6 + sn) Cos[2(cn - 1)Pi/c + 4 t/c],
			(c/6 + sn) Sin[2(cn - 1)Pi/c + 4 t/c],
			0},
		0.2,
		{t, 0, Pi/2},
		1,
		CurvePoints->10, TubePoints->tpn, Closed->False,
		Twist->N[twopi*skew], Torsion->N[twopi*screw]]

(**********************************************************************)
(* The *Cap* functions handle the construction of the caps of a plat. *)
(* They are called by Plat.                                           *)
(**********************************************************************)

topCapQ[st_] :=
	Table[
		Graphics3D[{GrayLevel[0.4],
			Line[Table[
				{Cos[theta]/2 + 2 m - 1/2,
				0,
				Sin[theta]/2},
			{theta,0,Pi,Pi/3}
			]]
		}],
		{m,1,st/2}
	]
				
				
bottomCapQ[st_,xn_] :=
	Table[
		Graphics3D[{GrayLevel[0.4],
			Line[Table[
					{Cos[theta]/2 + 2 m - 1/2,
					0,
					-xn - Sin[theta]/2},
				{theta,0,Pi,Pi/3}
			]]
		}],
		{m,1,st/2}
	]

topCapC[st_] :=
	Table[
		Graphics3D[{GrayLevel[0.4],
			Line[Table[
				{Cos[theta]/2 + 2 m - 1/2,
				0,
				Sin[theta]/2},
			{theta,0,Pi,Pi/15}
			]]
		}],
		{m,1,st/2}
	]
				
				
bottomCapC[st_,xn_] :=
	Table[
		Graphics3D[{GrayLevel[0.4],
			Line[Table[
					{Cos[theta]/2 + 2 m - 1/2,
					0,
					-xn - Sin[theta]/2},
				{theta,0,Pi,Pi/15}
			]]
		}],
		{m,1,st/2}
	]

topCapTN[st_,tpn_] :=
	Table[
		tblrNbhd[
			{Cos[theta]/2 + 2 m - 1/2,
			0,
			Sin[theta]/2},
			0.2,
			{theta, 0, Pi},
			0,
			CurvePoints->10, TubePoints->tpn, Closed->False
		],
		{m,1,st/2}
	]

bottomCapTN[st_,xn_,tpn_] :=
	Table[
		tblrNbhd[
			{Cos[theta]/2 + 2 m - 1/2,
			0,
			-xn - Sin[theta]/2},
			0.2,
			{theta, 0, Pi},
			0,
			CurvePoints->10, TubePoints->tpn, Closed->False
		],
		{m,1,st/2}
	]

(**********************************************)
(*          The top level functions.          *)
(**********************************************)

Components[b_] :=
	CycleLength[BraidToPermutation[b]] /; BraidQ[b] || Message[Braid::bogus, b]

OpenBraid[b_, options___] :=
	Module[{form = Style /. {options} /. Options[OpenBraid]},
		Switch[form,
			Automatic,
				Table[ dcOQ[b[[1]][[cn]], sn, cn],
					{sn, 1, b[[2]]}, {cn, 1, CrossingNumber[b]}
				],
			Curve,
				Table[ dcOC[b[[1]][[cn]], sn, cn],
					{sn, 1, b[[2]]}, {cn, 1, CrossingNumber[b]}
				],
			Ribbon,
				Table[ dcOTN[b[[1]][[cn]], sn, cn, 2],
					{sn, 1, b[[2]]}, {cn, 1, CrossingNumber[b]}
				],
			Tube,
				Table[ dcOTN[b[[1]][[cn]], sn, cn, 8],
					{sn, 1, b[[2]]}, {cn, 1, CrossingNumber[b]}
				],
			_,
				Message[Style::bogus, form, "OpenBraid"]
		]
	] /; BraidQ[b] || Message[Braid::bogus, b]

ClosedBraid[b_, options___] :=
	Module[{form, cf, ff, ps, nx, nc, cl, pb},
		{form, cf, ff}={Style, Canon, Framing} /. {options} /. Options[ClosedBraid];
		nc = Components[b];
		cl = CycleLengths[BraidToPermutation[b]];
		cb = ToCycles[BraidToPermutation[b]];
		If[ ff == Automatic, ff = Table[0, {nc}]];
		If[ Length[ff] < nc,
			ff = Join[ff, Table[0, {nc - Length[ff]}]]];
		nx = CrossingNumber[b];
		ps = Table[ BraidToPermutation[{Take[ b[[1]], n], b[[2]]}],
			{n, 0, nx - 1}];
		Switch[form,
			Automatic,
				Table[
					dcCQ[b[[1]][[cn]], sn, cn, CrossingNumber[b]],
					{sn, 1, b[[2]]}, {cn, 1, CrossingNumber[b]}
				],
			Curve,
				Table[
					dcCC[b[[1]][[cn]], sn, cn, CrossingNumber[b]],
					{sn, 1, b[[2]]}, {cn, 1, CrossingNumber[b]}
				],
			Ribbon,
				Table[
					dcCTN[ b[[1, k]], ps[[ k, cb[[i, j]] ]], k, nx, 2, cf/2,
						(j-1)*ff[[i]]/cl[[i]] + (k-1)*ff[[i]]/(nx*cl[[i]]),
						ff[[i]]/(nx*cl[[i]])
					],
					{i,1,nc},{j,1,cl[[i]]},{k,1,nx}
				],
			Tube,
				Table[
					dcCTN[ b[[1, k]], ps[[ k, cb[[i, j]] ]], k, nx, 8, cf/2,
						(j-1)*ff[[i]]/cl[[i]] + (k-1)*ff[[i]]/(nx*cl[[i]]),
						ff[[i]]/(nx*cl[[i]])
					],
					{i,1,nc},{j,1,cl[[i]]},{k,1,nx}
				],
			_,
				Message[Style::bogus, form, "ClosedBraid"]
		]
	] /; BraidQ[b] || Message[Braid::bogus, b]

Plat[b_, options___] :=
	Module[{form = Style /. {options} /. Options[Plat],
		s = 2 Ceiling[ b[[2]]/2 ]},
		Switch[form,
			Automatic,
				{
					topCapQ[s],
					Table[ dcOQ[b[[1]][[cn]], sn, cn],
						{sn, 1, s}, {cn, 1, CrossingNumber[b]}
					],
					bottomCapQ[s, CrossingNumber[b]]
				},
			Curve,
				{
					topCapC[s],
					Table[ dcOC[b[[1]][[cn]], sn, cn],
						{sn, 1, s}, {cn, 1, CrossingNumber[b]}
					],
					bottomCapC[s, CrossingNumber[b]]
				},
			Ribbon,
				{
					topCapTN[s, 2],
					Table[ dcOTN[b[[1]][[cn]], sn, cn, 2],
						{sn, 1, s}, {cn, 1, CrossingNumber[b]}
					],
					bottomCapTN[s, CrossingNumber[b], 2]
				},
			Tube,
				{
					topCapTN[s, 8],
					Table[ dcOTN[b[[1]][[cn]], sn, cn, 8],
						{sn, 1, s}, {cn, 1, CrossingNumber[b]}
					],
					bottomCapTN[s, CrossingNumber[b], 8]
				},
			_,
				Message[Style::bogus, form, "Plat"]
		]
	] /; BraidQ[b] || Message[Braid::bogus, b]


(***********************************************************
Mathematica procedure for calculating the multivariable Alexander polynomial of a link, presented as a braid on n strings. Original version written in Maple by H.R.Morton and Julian Hodgson, Liverpool University,  January 1996.  Translated into Mathematica from Maple by Scott Baldridge, May 1999.
************************************************************)


AlexPoly[b_] := Module[

{m,r,funct,stringpermutation,burausub,buraumatrix,buraucurrent,burauproduct,forsubtraction,i,j,k,i1,x,y,count,plusminus,ithstring,offset,onestringlinks,lonenumber,trash,trashcount,traspointer,orbitbook,longlink,strings,flag,flag1,flag2,temp,longlinkpos,div,poly,result,resultfactors,resultcompare,resultfinal, tempM, polyfactor},

(***********************************************************
To use this procedure for calculating the Alexander polynomial of
a closed braid, start Mathematica and then read in this file. You can calculate the polynomial for the closure of a braid on n strings with braid word sigma_1 sigma_2^{-1} sigma_3, for example, by typing braid(1,-2,3);
**************************************************************)

(*************************************************************
Note that the n-array T is global so that when the final expression  is  returned to  the user it  can be manipulated. When the procedure has finished, besides the information on screen, the array T[1]..T[n] of variables in the eventual output can be altered to allow for substitutions of your choice.

Not implemented: The characteristic polynomial of the Burau matrix is also global, and is stored as charpoly, in terms of the variables in T, and a variable X.
**************************************************************)

Clear[T];

braidword = b[[1]];
narg = Length[braidword];

(* narg is the number of strings, need to check that braid group word does not extend past the number of strings *)

n=b[[2]];


(*
For[i=1, i <= narg, i++, If[ Abs[braidword[[i]]]> (n-1), Return["Error: String number out of range!"]]];

(* Need to check that the braid element is not trivial *)

If[narg < 1, Return["Error: No braid element given!"]];
*)


(************************************************************
Works out order of strings at the bottom of the braid. We note that the ordering is taken to be 1,...,n at the top of the braid (opposite to the convention in Preprint, H.R. Morton, `The Multivariable Alexander Polynomial for a Closed Braid`, 1998). The ordering at the bottom of the braid is then stored in string_permutation:
**************************************************************)

For[i=1,i<= narg,i++, stringpermutation[i]=i];

For[count=1, count<=narg, count++,
	i=braidword[[count]];
	ithstring=Abs[i];
	i1 = stringpermutation[ithstring];
	stringpermutation[ithstring] = stringpermutation[ithstring+1];
	stringpermutation[ithstring+1] = i1;
	];

(* Print[Table[stringpermutation[i], {i,1,narg}]]; *)

(* Create meridians T[1], T[2], ...., T[n] *)

lonenumber=0;
trashcount=0;
longlink=0;
For[i=1, i<= n, i++, 
	For[count =1, count<=(n+1), count++, 
		orbitbook[i,count]=0;
	]
];


(***********************************************
This rather long and complicated collection loops is responsible for ensuring that the correct T[i]'s are used in the final polynomial.
************************************************)

strings=0;

While[!(strings == n),
	strings++;
	flag=0;
	flag1=0;
	trashpointer=0;
	
	
	If[!(trashcount == 0),
	{	
		While[flag==0,
			trashpointer++;
			
			If[trashpointer == trashcount, flag=1];
			If[trash[trashpointer]==strings, flag=1;flag1=1];
			
		]; (* ends while statement *)
	}]; (*ends if statement*)

	If[flag1==0,
	{

		
		If[stringpermutation[strings]==strings,
		{
			lonenumber++;
			onestringlinks[lonenumber]=strings;
			
		},
		(*else*)
		{	
			flag=0;
			temp=strings;
			longlinkpos=1;
			longlink++;
			orbitbook[longlink,longlinkpos]=strings;
		
			
			While[flag ==0,
			{
				
				temp=stringpermutation[temp];
								
				If[!(temp == strings),
				{
					longlinkpos++;
					orbitbook[longlink,longlinkpos]=temp;
					trashcount++;
					trash[trashcount]=temp;
					
				},
				(*else*)
				{
					flag=1;
				}]; (*ends If statment*)
					
			}]; (*ends While Statement*)
		}
		];(*Ends If Statement*)
	}
	]; (*ends If statement*)

]; (*ends the while statement*)


(* Print[Table[orbitbook[i,j],{i,1,n},{j,1,n+1}]]; *)

For[i=1, i<=longlink, i++,
	{
	longlinkpos=1;
	flag=0;
	While[flag==0,
	{
		longlinkpos++;
		If[orbitbook[i,longlinkpos] == 0,
		{	
			flag=1;
		}, (*else*)
		{
			T[orbitbook[i,longlinkpos]]=T[orbitbook[i,1]];
		}]; (*end of if*)
	}]; (*end of While*)
}]; (*end of FOR statement*)

(****** Defines 3 by 3 submatrix used forming burau matrices ****)

tempM = IdentityMatrix[3];
tempM[[2,1]] = r; tempM[[2,2]] = -r; tempM[[2,3]] =1;
burausub = {Inverse[tempM], tempM};

buraumatrix = Table[IdentityMatrix[n], {i,1,2}, {j,1,(n-1)}]; 

(***** Create burausub matrices ****)

For[j=1, j<=2, j++,
	For[k=1,k<=(n-1), k++,
	{		
		If[k==1, offset=2, offset=1];
		
		For[x=offset, x<=3, x++,
			For[y=offset, y<=3, y++,
				{buraumatrix[[j,k,k+x-2,k+y-2]] = burausub[[j,x,y]]}
			]
		];
			
	}]; (*End of For statment*)
]; (*End of For statement*)

(* Print[MatrixForm[buraumatrix]]; *)

(****************************************************************
Here we work out each individual burau matrix with the corresponding meridian as its variable and keep a running total of the product.
****************************************************************)

burauproduct = IdentityMatrix[(n-1)];
forsubtraction = IdentityMatrix[(n-1)];
For[i=1,i<= narg,i++, stringpermutation[i]=i];


For[count=1, count<=narg, count++,
{
	i=braidword[[count]];
	ithstring = Abs[i];
	plusminus = (3+ithstring/i)/2; 
 
 	If[ Sign[i]==1,
	{
		r = T[stringpermutation[ithstring+1]];
	}, (*else *)
	{
		r = T[stringpermutation[ithstring]];
	}]; (* End of IF Statement *)

	buraucurrent = buraumatrix[[plusminus,ithstring]];
	buraucurrent = buraucurrent[[ Range[1, n-1], Range[1,n-1] ]];
  	(* Print[i, " ", ithstring, " ", plusminus]; *)

(**********************************************************
At this stage A is a matrix whose values are integers or functions of r. 
***********************************************************)

	buraucurrent = Map[	Evaluate, buraucurrent];

	
(**********************************************************
Now "map" has forced A to take on baord the actual values of r:=T[i].
************************************************************)
	
	If[Sign[i] == -1, 
	{
		buraucurrent = r*buraucurrent;
		forsubtraction = r*forsubtraction;
	}];
	

(***********************************************************
Here we multiply the matrix by the meridian whose value it takes only if it is the inverse form. This keeps fractions out of the  calculation rendering  it quicker. We also  keep track of all the factors we have used for later.
************************************************************)

	burauproduct = buraucurrent.burauproduct;
	i1 = stringpermutation[ithstring];
	stringpermutation[ithstring] = stringpermutation[ithstring +1];
	stringpermutation[ithstring +1] = i1;
	
}]; (*End of For *)

div = 1;

For[i=1, i<=n, i++,	div = div*T[i] ];

(* Print[Table[T[i], {i,1,n}]]; *)

div = 1-div;

poly = burauproduct - forsubtraction;
poly = Det[poly];

(* Print[Factor[poly]]; *)

div = Factor[div];
(* Print[div]; *)

result = poly / div;
result = Factor[result];


If[lonenumber==0, result = Simplify[result * (1-T[1])]];


(*********************************************************
This routine removes any factors of T[i] and should be removed if extra time is needed.  It takes a significant amount of time to run this part.
************************************************************)

polyfactor = FactorList[result];

For[i=1, i<=Length[polyfactor], i++,
{
	count=0;
	flag1=0;
	While[flag1==0,
	{
		count++;
		If[polyfactor[[i,1]] == T[count],
		{
			result = result / (T[count]^polyfactor[[i,2]]);
			flag1=1;
		}];
		If[count == n, flag1=1];
	}];
}];

result = Factor[result];
Return[result];

]  /; BraidQ[b] || Message[Braid::bogus, b] 


End[]         (* end the private context *)
EndPackage[]  (* end the package context *)