(* :Title: ParametrizeCurve *)

(* :Author: Steven Wilkinson, Northern Kentucky University *)

(* :Summary:
    ParametrizeCurve gives a parametrization of an implicitly defined curve. *)

(* :Context: `ParametrizeCurve` *)

(* :Copyright: © 1997 by Steven Wilkinson *)

(* :Package Version: 1.2.2 *)

(* :Mathematica Version: 2.- 5. *)

(* :Keywords: curve, implicit equations, parametrizations *)


BeginPackage["ParametrizeCurve`"];

Needs["Utilities`FilterOptions`"];

ParametrizeCurve::usage =
	"ParametrizeCurve[equations,variables,{t,tmin,tmax}] generates a parametrization for the curve defined implicitly by equations in the coordinates variables.";

VectorNorm::usage =
	"VectorNorm is an option for ParametrizeCurve that gives the norm used to calculate lengths.";

ParametrizationSpeed::usage =
	"ParametrizationSpeed is an option for ParametrizeCurve that gives the desired magnitude of the derivative of the parametrization.";

MinNorm::usage =
	"MinNorm is an option for ParametrizeCurve. Points on the curve with norm of the cross product of the gradients of the implicit equations under MinNorm are automatically assumed to be near a singular point.";

MaxNorm::usage =
	"MaxNorm is an option for ParametrizeCurve. Points on the curve with norm of the cross product of the gradients of the implicit equations over MaxNorm are automatically assumed to be far from singular points.";

InitialDirection::usage =
	"InitialDirection is an option for ParametrizeCurve. Its value is a vector that points in the same general direction that initial tangent vector points.";

CheckSingularPoints::usage =
	"CheckSingularPoints is an option for ParametrizeCurve. If set to True ParametrizeCurve will try to accurately parametrize through singular points. If set to False ParametrizeCurve will not check for singular points.";

UnitSpeed::usage =
	"UnitSpeed can be used in the option ParametrizationSpeed of ParametrizeCurve to specify the parametrization is to have constant speed 1.";


GradientSpeed::usage =
	"GradientSpeed can be used in the option ParametrizationSpeed of ParametrizeCurve to specify the parametrization's speed is equal to the magnitude of the cross product of the gradient vectors of the implicit equations for the curve.";

AdjustInitialPoint::usage =
	"AdjustInitialPoint is an option for ParametrizeCurve. If set to True, whenever the given initial point is not on the curve, the nearest point on the curve to this point will be used as the initial point.";

ParametrizeCurve::arg =
	"ParametrizeCurve must have at least 3 arguments where the first argument is a list of the implicit equations defining the curve and coordinates for the initial point, the second argument is a list of the coordinate variables, and the third argument is the parameter and its range, {t,tmin,tmax}.";

ParametrizeCurve::iniParameter =
	"The initial parameter value `1` must be a real number.";

ParametrizeCurve::eqnNo =
	"There must be `1` equations to implicitly define a curve in `2` coordinates.";

ParametrizeCurve::iniPt1 =
	"The list of defining equations, `1`, must include a valid initial point.";

ParametrizeCurve::iniPt2 =
	"The initial point `1` does not lie on the curve given by `2`.";

ParametrizeCurve::iniPt3 =
	"Using `1` as the initial point.";

ParametrizeCurve::iniPt4 =
	"The curve is not differentiable at the initial point `1`.";

ParametrizeCurve::singPt =
	"The point `1` appears to be near a degenerate singular point of the curve.";

ParametrizeCurve::isolatedPt =
	"The initial point `1` is an isolated singular point of the curve.";

Options[ParametrizeCurve] =
	{VectorNorm->2,
	ParametrizationSpeed->UnitSpeed,
	MinNorm->0,
	MaxNorm->Infinity,
	InitialDirection->Automatic,
	CheckSingularPoints->True,
	AdjustInitialPoint->True};

Begin["`Private`"];

(* NumericQ[] and Cross[] are functions built-in to version 3. but not version 2. *)
If[$VersionNumber<3.,
	NumericQ[expr_]:=NumberQ[N[expr]];
	Cross[m__]:=
		Module[{mt},
			mt=Transpose[{m}];
			Map[(-1)^(#+1)Det[Drop[mt,{#}]]&,Range[Length[mt]]]
		]
];

(* Norm[v,p] is a built-in function in version 5, but not before. *)
If[$VersionNumber<5.,
	Norm[v_,Infinity]:=Max[Abs[v]];
	Norm[v_,p_]:=Apply[Plus,Abs[v^p]]^(1/p)];

ParametrizeCurve[equations_List,variables_List,{t_,tmin_?NumericQ,tmax_?NumericQ},opts___] :=
	Module[{adjustIniPt,adjustInitialPointQ,checkSingularPointsQ,ComplexQ,
      directionIndex,dX,epsilon,eqns,equationNumber,firstDirection,
      firstIndex=1,firstOrderDirections,functions,Gradient,gradients,hessians,
      initialPoint,initialPtEqn,initialPtEqnStandard,initialValue,J,
      lastDirection,lastIndex=1,lastNorm=0,lastSwitchQ=False,lastT,length,
      maxDiv=0,maxEpsilon,minEpsilon,NInitialPoint,norm,NReal,perpFunction,
      rearrangeFunction,secondOrderDirections,speed,switchQ,t0,thisEpsilon,
      variableIndex,variableNumber,x},

		x = variables /. {a_[t]->a};
		variableNumber = Length[x];
		variableIndex = Range[variableNumber];
		equationNumber = variableNumber-1;

(* Options *)
        norm[v_]:= Norm[v,VectorNorm/.{opts}/.Options[ParametrizeCurve]];
		speed[s_]=
			If[
				(ParametrizationSpeed /.{opts} /.
					Options[ParametrizeCurve]) === GradientSpeed,
					Identity,
					ParametrizationSpeed /. {opts} /. Options[ParametrizeCurve] /.
						{t->s} /. {a_?(!NumericQ)->a[s]} /. {UnitSpeed->1} /.
						{a_->((a*#/norm[#])&)}
			];
		thisEpsilon = MinNorm /. {opts} /. Options[ParametrizeCurve];
		maxEpsilon = MaxNorm /. {opts} /. Options[ParametrizeCurve];
		firstDirection = InitialDirection /. {opts} /.
			Options[ParametrizeCurve] /. {Automatic->Table[0,{i,variableNumber}]};
		checkSingularPointsQ = CheckSingularPoints /. {opts} /.
			Options[ParametrizeCurve];
		rearrangeFunction = RearrangeFunction /. {opts} /.
			Options[ParametrizeCurve] /. {RearrangeFunction->Identity};

        adjustInitialPointQ=AdjustInitialPoint/.{opts}/.Options[ParametrizeCurve];
        
(* General command definitions. *)
		Gradient[func_,u_] := Map[D[func,#]&,u];
		J = {{0,-1},{1,0}};
		ComplexQ[z_] := (Head[z]==Complex);
		NReal[0] = 0.;
		NReal[z_?NumericQ] := N[z];
		NReal[z_List] := Map[NReal,z];

(* Initial point information. *)
		initialPtEqn =
			Select[equations,(#	/.	{(a_[b_?NumericQ]==c_?NumericQ)	:>
				True,(c_?NumericQ==a_[b_?NumericQ])	:>	True})&
			];

		initialPtEqnStandard = initialPtEqn /. {Equal[a_,b_[c_]]->Equal[b[c],a]};
		initialPoint = x /. (initialPtEqnStandard/.{Equal[b_[c_],a_]->Rule[b,a]});
		NInitialPoint=NReal[initialPoint];
		t0=Flatten[
			Union[Map[Level[initialPtEqnStandard[[#]],{2}]&,variableIndex]]];
		t0=First[NReal[t0]];
		lastT=t0-1;

(* Information on the generators of the curve. *)
		functions =
			Evaluate[
				Map[
					Function,
					Flatten[{(Complement[equations,initialPtEqn] /.
						MapThread[Rule,{x,Map[Slot,variableIndex]}]) /.
						{Equal->Subtract,a_[t]->a}}
					]
				]
			];
		gradients = Map[(Gradient[#,x]&),Map[#@@x&,functions]];
		hessians =
			Function[Evaluate[Map[Gradient[#,x]&,gradients,{2}] /.
				MapThread[Rule,{x,Map[Slot,variableIndex]}]]
			];
		perpFunction =
			Function[Evaluate[rearrangeFunction[Apply[Cross,gradients]] /.
				MapThread[Rule,{x,Map[Slot,variableIndex]}]]
			];
		length =
			Function[Evaluate[rearrangeFunction[norm[perpFunction@@x]] /.
				MapThread[Rule,{x,Map[Slot,variableIndex]}]]
			];
(* Error messages *)

(* Return error message if the number of equations is not one less than the
    number of coordinates. *)
		If[
			Length[functions]!=equationNumber,
			Return[
				Message[ParametrizeCurve::eqnNo,equationNumber,variableNumber]
			]
		];

(* Return error message if an initial point is not given. *)
		If[
			Length[NInitialPoint]!=variableNumber,
			Return[Message[ParametrizeCurve::iniPt1,equations]]
		];

(* Return error message if the initial point does not lie on the curve. Adjust
    the initial point to lie on the curve if the option allows it. *)
		initialValue = Map[#@@NInitialPoint&,functions];
		If[
			Chop[initialValue.initialValue]!=0,
			Message[ParametrizeCurve::iniPt2,
				initialPoint,Complement[equations,initialPtEqn]
			];
			If[
				adjustInitialPointQ,
				adjustIniPt[r_] :=
					(r-Map[#@@r&,functions].
						(Inverse[Outer[Dot,gradients,gradients,1]].
						gradients) /. MapThread[Rule,{x,r}]);
				NInitialPoint = FixedPoint[adjustIniPt,NInitialPoint,10];
				initialPtEqnStandard =
					MapThread[ReplacePart[#1,#2,-1]&,
						{initialPtEqnStandard,NInitialPoint}
					];
				Message[ParametrizeCurve::iniPt3,NInitialPoint];
			]
		];

(* Return error message if the initial parameter value is not real. *)
		If[Head[t0]!=Real,
			Return[Message[ParametrizeCurve::iniParameter,t0]]
		];

(* Return error message if the curve is not differentiable at the initial
    point. *)
		If[
			!MatrixQ[(gradients/.MapThread[Rule,{x,NInitialPoint}]),NumericQ],
			Return[Message[ParametrizeCurve::iniPt4,NInitialPoint]]
		];

(* The possible tangent directions if the point is assumed to be
    nonsingular. *)
		firstOrderDirections[s_] =
			Function[
				Evaluate[(
					Map[#*rearrangeFunction[speed[s][perpFunction@@x]]&,
						{-1,1}]) /. MapThread[Rule,{x,Map[Slot,variableIndex]}]
				]
			];

(* The possible tangent directions if the point is assumed to be singular. *)
		secondOrderDirections[s_][u__] :=
			Module[{U,sigma,V,K,VLast2,VLast2T,hess,eigvectors,fillOut},
				fillOut[m_,v_] := m /; Length[m]==equationNumber;
				fillOut[m_,v_] :=
					Module[{projection},projection=v-Plus@@Map[(#.v)#&,m];
						fillOut[Join[m,{projection/Sqrt[projection.projection]}],
							RotateLeft[v]
						] /; Chop[Plus@@Map[(#.v)#&,m]]!=Chop[v]
					];
				fillOut[m_,v_] := fillOut[m,RotateLeft[v]];
				{U,sigma,V} = SingularValues[gradients/.MapThread[Rule,{x,{u}}]];
				V = fillOut[V,Join[{1},Table[0,{equationNumber}]]];
				VLast2 = {Part[V,-1],Cross@@V};
				VLast2T = Transpose[VLast2];
				hess = Flatten[Take[U,-1].hessians[u],1];
				If[Re[Det[hess]]<$MachineEpsilon&&Length[sigma]<(equationNumber-1),
					Message[ParametrizeCurve::singPt,{u}]
				];
				eigvectors =
					Re[Map[speed[s],
						Eigenvectors[VLast2T.J.VLast2.hess.VLast2T.VLast2]]
					];
				Return[Join[eigvectors,-eigvectors]]
			];

(* Chooses the position of the desired tangent vector in the list of possible
    ones. *)
		directionIndex[point_,vectors_] :=
			Module[{newIndex},
				If[point==NInitialPoint,lastDirection=firstDirection];
				newIndex =
					First[
						Sort[Range[Length[vectors]],
							Less[Chop[(vectors[[#1]]-lastDirection).
							Conjugate[vectors[[#1]]-lastDirection]],
							Chop[(vectors[[#2]]-lastDirection).
							Conjugate[vectors[[#2]]-lastDirection]]]&]];
				lastIndex=newIndex;
				lastDirection=vectors[[newIndex]];
				Return[newIndex]
			];

(* The tangent vector used. *)
		dX[u_,s_?NumericQ,i_] :=
			Module[{possibleDirections,index},
				possibleDirections =
					If[
						switchQ[u,s],
						secondOrderDirections[s]@@u,
						firstOrderDirections[s]@@u
					];
				index = directionIndex[u,possibleDirections];
				Return[possibleDirections[[index,i]]]
			];

(* Determines whether the point is sufficiently near a singular point to
    switch between the nonsingular and singular point algorithms. *)
		switchQ[u_,s_] :=
			Module[{possibleEpsilon,thisNorm,thisSwitchQ=False},
				thisNorm = length@@Re[u];
				If[thisNorm==0,Return[True]];
			(* Definitely a singular point; switch *)
				If[s!=lastT,
					lastT=s;
					If[s!=t0,
						possibleEpsilon=(-thisNorm+lastNorm),
						possibleEpsilon=-thisNorm
					];
					lastNorm=thisNorm;
					If[possibleEpsilon>thisEpsilon &&
							possibleEpsilon<=maxEpsilon,
						thisEpsilon=possibleEpsilon];
					If[thisNorm<(thisEpsilon+possibleEpsilon)/2,
						thisSwitchQ=True,
						thisSwitchQ=False;
					];
					lastSwitchQ = thisSwitchQ,
					
					thisSwitchQ=lastSwitchQ
				];
				Return[thisSwitchQ]
			];

(* If you don't want to check for singular points,this does the integration. *)

		If[!checkSingularPointsQ,
			Module[{index,possibleDirections},
				possibleDirections = firstOrderDirections[t]@@Map[#[t]&,x];
				index = directionIndex[NInitialPoint,possibleDirections];
				Return[
					NDSolve[Flatten[
						Join[Map[(x[[#]]'[t]==possibleDirections[[index,#]])&,
							variableIndex],
						initialPtEqnStandard]],variables,{t,tmin,tmax},
						Utilities`FilterOptions`FilterOptions[NDSolve,opts]
					]//First
				]
			]
		];

(* If you do want to check for singular points, this does the integration. *)

(* Checks to see if theinitial point is an isolated point of the curve. *)
		If[NReal[Chop[length@@NInitialPoint]]==0. &&
				Or@@Map[VectorQ[#,ComplexQ]&,
			Chop[secondOrderDirections@@NInitialPoint]],
			Message[ParametrizeCurve::isolatedPt,initialPoint];
			If[variables==x,
				Return[{MapThread[Rule,{variables,Map[Function,NInitialPoint]}]}],
				Return[{MapThread[Rule,{variables,NInitialPoint}]}]
			]
		];
		
		NDSolve[
			Flatten[Join[
				Map[(x[[#]]'[t]==dX[Map[#[t]&,x],t,#])&,variableIndex],
				initialPtEqnStandard
			]],
			variables,{t,tmin,tmax},
			Utilities`FilterOptions`FilterOptions[NDSolve,opts]]//First
		];

ParametrizeCurve[exp___] :=
	Message[ParametrizeCurve::arg];

End[]; (* Private *)

EndPackage[];  (* ParametrizeCurve *)
