(* ::Package:: *)

BeginPackage["quadrisecant`"]

(* version 1.0 *)

(* variables *)
KVS;  (* vertices of the knot *)
QSS;  (* quadrisecants of the knot *)
K14;	(* the pre-inputed knot *)
K14Generic;
OQSS;	(* the quadisecant points in the order along the knot *)

(* functions *)
QuadriInit;
InputKnot;
DrawKnot;
CheckGeneralPosition;
GetQuadrisecantsForEdges;
GetQuadrisecantsForKnot;
DrawQuadrisecantApproximation;

Begin["`Private`"]

CearSymbols[] := Module[
	{},
	Clear[
		KVS,
		K14,
		QSS,
		OQSS
	];
	Clear[
		QuadriInit,
		InputKnot,
		DrawKnot,
		CheckGeneralPosition,
		GetQuadrisecantsForEdges,
		GetQuadrisecantsForKnot,
		DrawQuadrisecantApproximation
	];
]

QuadriInit[] := Module[
	{},
	K14 = {{-(3/5),-(1/5),0},{-(19/25),11/50,-(99/25)},{228/25,-(66/25),-(112/25)},{143/20,-(121/40),-(109/20)},{-(13/10),11/20,-(61/10)},{-1,1/10,-(171/20)},{10,-1,-8},{10,-1,1},{28/5,6/5,-(1/10)},{81/10,1/20,-(21/20)},{59/10,-(21/20),1/20},{12,-(11/5),-(33/10)},{12,1/5,3/10},{6,2,0}};
]

(* check whether the knot is in general position *)
CheckGeneralPosition[] := Module[
	{i, j, k, A, B, C, D, E, F, vnum},
	vnum = Length[KVS];
	(*check whether any four vertices are coplanar*)
	For[i=1, i<vnum, i++,
		For[j=i+1, j<vnum, j++, 
			For[k=j+1, k<vnum, k++,
				For[l=k+1, l<=vnum, l++,
					If[Det[{KVS[[j]]-KVS[[i]], KVS[[k]]-KVS[[i]], KVS[[l]]-KVS[[i]]}]==0,
						Print["The four vertices ", i, ", ", j, ", ", k, " and ", l, " are coplanar"]
					]
				]
			]
		]
	];

	(*check whether any three edge vectors are coplanar*)
	For[i=1, i<vnum, i++,
		For[j=i+1, j<vnum, j++, 
			For[k=j+1, k<vnum, k++,
				A = KVS[[i]];
				B = KVS[[i+1]];
				C = KVS[[j]];
				D = KVS[[j+1]];
				E = KVS[[k]];
				F = If[k==Length[KVS], KVS[[1]], KVS[[k+1]]];
				If[Det[{B-A, D-C, F-E}]==0, 
					Print["The edge vectors ", i, ", ", j, " and ", k, " are coplanar"]
				] 
			]
		]
	];

	Print["The knot is in general position"];
]

InputKnot[VS_] := Module[
	{},
	KVS = VS;
	QSS = {};
	OQSS = {};
]

GetQuadrisecantsForEdges[i_, j_, k_, l_] := Module[
	{A, B, C, D, E, F, G, H},
	A = KVS[[i]];
	B = KVS[[i+1]];
	C = KVS[[j]];
	D = KVS[[j+1]];
	E = KVS[[k]];
	F = KVS[[k+1]];
	G = KVS[[l]];
	H = If[l==Length[KVS], KVS[[1]], KVS[[l+1]]];
	(* Print[A, B, C, D, E, F, G, H]; *)
	GetQuadrisecantsForEdges[{i,j,k,l}, A, B, C, D, E, F, G, H];
]

GetQuadrisecantsForKnot[] := Module[
	{vnum, A, B, C, D, E, F, G, H, i, j, k, l, qs},
	vnum = Length[KVS];
	(*Print[VS, vnum];*)
	QSS = {};

	For[i=1, i<vnum, i++,
		Print["i=", i];
		For[j=i+1, j<vnum, j++, 
			For[k=j+1, k<vnum, k++,
				For[l=k+1, l<=vnum, l++,
					A = KVS[[i]];
					B = KVS[[i+1]];
					C = KVS[[j]];
					D = KVS[[j+1]];
					E = KVS[[k]];
					F = KVS[[k+1]];
					G = KVS[[l]];
					H = If[l==vnum, KVS[[1]], KVS[[l+1]]];
					(* Print[A, B, C, D, E, F, G, H]; *)
					GetQuadrisecantsForEdges[{i,j,k,l}, A, B, C, D, E, F, G, H];
				]
			]
		]
	];
	Print["There are ", Length[QSS], " quadrisecants"];
]

GetQuadrisecantsForEdges[IDX_, A_, B_, C_, D_, E_, F_, G_, H_] := Module[
	{a, b, c, d, e, f, g, h, aa, bb, cc, xx, p, q, r, s, x, y, qs, P, Q, R, S, SI},
	SI = False;
	If[SI, Print[IDX, A, B, C, D, E, F, G, H]];
	(* p = a/(1-x) + b = c/(1-y) + d*)
	If[Det[{B-A, D-C, F-E}]==0,
		Print[IDX, A, B, C, D, E, F, " the three vectors AB, CD, EF are coplanar"];
		Return[]
	];
	a = Det[{E-C, D-C, F-E}] / Det[{B-A, D-C, F-E}];
	b = Det[{C-A, D-C, F-E}] / Det[{B-A, D-C, F-E}];
	If[Det[{B-A, D-C, H-G}]==0,
		Print[IDX, A, B, C, D, G, H, " the three vectors AB, CD, GH are coplanar"];
		Return[]
	];
	c = Det[{G-C, D-C, H-G}] / Det[{B-A, D-C, H-G}];
	d = Det[{C-A, D-C, H-G}] / Det[{B-A, D-C, H-G}];
	(* q = e * (x-1)/x + f = g * (y-1)/y + h *)
	e = Det[{B-A, A-E, F-E}] / Det[{B-A, D-C, F-E}];
	f = Det[{B-A, E-C, F-E}] / Det[{B-A, D-C, F-E}];
	g = Det[{B-A, A-G, H-G}] / Det[{B-A, D-C, H-G}];
	h = Det[{B-A, G-C, H-G}] / Det[{B-A, D-C, H-G}];

	If[SI, Print[a, b, c, d, e, f, g, h]];

	(* x satisfies aa x^2 + bb x + cc \[Equal] 0*)
	aa = b*e-c*e-d*e+b*f-c*f-d*f+c*g-b*h+c*h+d*h;
	(* If[aa==0, Print[aa]; Return[]]; *)
	bb = -a*e-2*b*e+2*c*e+2*d*e-a*f-b*f+c*f+d*f-c*g+a*h+b*h-c*h-d*h;
	cc = a*e+b*e-c*e-d*e;
	If[SI, Print["equation for x: ", aa, "x^2+", bb, "x+", cc]];
	If[ bb^2 - 4*aa*cc < 0, Return[{}]];
	xx = Solve[ aa * x^2 + bb * x + cc == 0, x][[All, 1, 2]];
	(*Print["xx=", xx];*)

	qs = {};
	If[SI, Print[xx, xx[[1]]]];
	For[i=1, i<=Length[xx], i++, 
		x = xx[[i]];
		If[x==1 || x==0, Continue[]];
		If[a+b-d-b*x+d*x == 0, Continue[]];
		y = (a+b-c-d-b*x+c*x+d*x)/(a+b-d-b*x+d*x);
		If[y==1, Continue[]];
		p = a / (1-x) + b;
		q = e * (x-1)/x + f;
		r = (Det[{B-A, D-C, C-A}] * x + Det[{B-A, D-C, A-E}])/(Det[{B-A, D-C, F-E}]);
		s = (Det[{B-A, D-C, C-A}] * y + Det[{B-A, D-C, A-G}])/(Det[{B-A, D-C, H-G}]);
		(*Print["p=", p, ", q=", q, ", r=", r, ", s=", s];*)
		If[p<=0 || p>=1 || q<=0 || q>=1 || r<=0 || r>=1 || s<=0 || s>=1, Continue[]];
		P = A + p*(B-A);
		Q = C + q*(D-C);
		R = E + r*(F-E);
		S = G + s*(H-G);
		AppendTo[QSS, {IDX, x, y, p, q, r, s, P, Q, R, S}];
		If[SI, Print[QSS]];
	]
]

DrawKnot[] := Module[
	{i, vnum, ps, v, vec},
	vnum = Length[KVS];
	ps = {};
	ps = Append[KVS, KVS[[1]]];
	(*
	For[i=1, i<=vnum, i++,
		v = KVS[[i]];
		vec = If[i==vnum, KVS[[1]]-v, KVS[[i+1]]-v];
		ps = Join[ps, Table[v + (t/200)*vec, {t, 0, 199}]];
	];*)
	Graphics3D[{CapForm[None], 
		Lighter@Orange, Tube[ps, 0.05], 
		FaceForm[None, Glow[White]], Tube[ps, 0.15]},  Boxed -> False
	]
]

(* Draw the quadrisecant approximation *)
DrawQuadrisecantApproximation[] := Module[
	{pos, i, idx, vs, from, to, vnum, k, p},
	pos = {};
	For[i=1, i<=Length[QSS], i++,
		AppendTo[pos, {QSS[[i, 1, 1]], N[QSS[[i, 4]]], QSS[[i, 8]]}];
		AppendTo[pos, {QSS[[i, 1, 2]], N[QSS[[i, 5]]], QSS[[i, 9]]}];
		AppendTo[pos, {QSS[[i, 1, 3]], N[QSS[[i, 6]]], QSS[[i, 10]]}];
		AppendTo[pos, {QSS[[i, 1, 4]], N[QSS[[i, 7]]], QSS[[i, 11]]}];
	];
	pos = Sort[pos];
	(* Print[Simplify[pos]]; *)
	OQSS = {};
	vnum = Length[KVS];
	For[i=1, i<=Length[pos], i++,
		AppendTo[OQSS, pos[[i,3]]];
	];
	
	AppendTo[OQSS, OQSS[[1]]];
	(* Print[Simplify[vs]]; *)
	Graphics3D[{CapForm[None], 
		Lighter@Orange, Tube[OQSS, 0.05], 
		FaceForm[None, Glow[White]], Tube[OQSS, 0.15]},  Boxed -> False
	]
]

End[];

EndPackage[]