# $Source: /u/maple/research/lib/algcurves/src/RCS/plot_knot,v $
# $Notify: mvanhoei@daisy.uwaterloo.ca $

#########################
#   Singularity knots   #
#########################

macro(	plot_knot=`algcurves/plot_knot`,
	all_substitutions=`algcurves/pl_allsub`,
	curve=`algcurves/pl_curve`,

	`puiseux/technical_answer`=`algcurves/puiseux_te`,

	g_conversion1=`algcurves/g_conversion1`,
	g_conversion2=`algcurves/g_conversion2`,
	g_expand=`algcurves/g_expand`,
	g_ext=`algcurves/g_ext`
):

# All possible floating point substitutions for algebraic numbers.
# v: a list of RootOf's
all_substitutions:=proc(v)
	local w,i,j;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	if v=[] then RETURN({v}) fi;
	w:=evalf(allvalues(v));
	{seq([seq(v[i]=j[i],i=1..nops(v))],j=w)}
end:

# f is a polynomial in x en y
plot_knot:=proc(f,x,y)
	global color,colours,epsilon,g_conversion1,g_conversion2,radius,tubepoints;
	local ext_f,pui,ext_pui,i,pui2,curven,t,opt,eps,cols;
option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Authors: H. Derksen and M. van Hoeij`;

if has(evala(Content(f,y)),x) then
	# If f has a component which is independent of y, this component
	# will be ignored by the puiseux algorithm, but I want to
	# include it in the computation. This is a fix for this case
	RETURN(procname(subs(x=x+y,f),args[2..nargs]))
fi;

	opt:={op(4..nargs,[args])};
	# default options
if not member(numpoints,indets(opt)) then opt:=opt union {numpoints=150} fi;
if not member(radius,indets(opt)) then opt:=opt union {radius=0.05} fi;
if not member(tubepoints,indets(opt)) then opt:=opt union {tubepoints=5} fi;
if not member(scaling,indets(opt)) then opt:=opt union {scaling=CONSTRAINED} fi;
if not member(style,indets(opt)) then opt:=opt union {style=PATCHNOGRID} fi;
	eps:=1;
	if member(epsilon,indets(opt)) then
		eps:=subs(opt,epsilon);
		opt:=opt minus {epsilon=eps}
	fi;
	cols:=[];

	# A quick fix to bring plot_knot in line with the options for plot

	if member( colour, indets(opt) ) then
		cols:=subs(opt,colour);
		opt:=opt minus {colour=cols}
	end if;

	if member( color, indets(opt) ) then
		cols:=subs(opt,color);
		opt:=opt minus {color=cols}
	end if;

	if not type( cols, 'list' ) then
		cols := [cols];
	end if;

	ext_f:=g_ext(f);
	pui2:=`puiseux/technical_answer`(g_expand(subs(g_conversion1,f),ext_f)
	 ,x,y,0,ext_f);
	pui:={};
	for i in pui2 do if coeff(i[1],x,0)=0 and ldegree(i[1],x)>=0 then
		pui:=pui union {i}
	fi od;
	ext_f:=subs(g_conversion2,ext_f);
	ext_f:=all_substitutions(ext_f);
	ext_f:=ext_f[1];
	ext_pui:=subs(ext_f,subs(g_conversion2,g_ext(pui)));
	ext_pui:=all_substitutions(ext_pui);
	pui:=subs(ext_f,subs(g_conversion2,pui));
	# Now replace the RootOf's by all corresponding floating point numbers
	pui2:={};
	for i in ext_pui do
		pui2:=pui2 union subs(i,pui)
	od;
	curven:={};
	for i in pui2 do
		curven:=curven union {curve(subs(x=x*evalf(
		 lcoeff(i[3],x)^(-1/degree(i[3],x))),i[1])
		,degree(i[3],x),t,eps,x)}
	od;
	if printlevel>1 then print(`Number of branches:`,nops(curven)) fi;
	if cols=[] then
		plots['tubeplot'](curven,t=0..2*Pi,op(opt))
	else
		curven:=[op(curven)];
		plots['display']([seq(plots['tubeplot'](curven[i],t=0..2*Pi,op(opt),
		 color=cols[1+irem(i,nops(cols))]),i=1..nops(curven))])
	fi
end:

curve:=proc(functie_in_x,d,t,eps,x) local f,fre,fim,factor;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	f:=expand(subs(x=eps*(cos(t)+I*sin(t)),functie_in_x));
        fre:=Re(f)/eps;fim:=Im(f)/eps;
	factor:=1/(sqrt(1+fre^2+fim^2)-fre);
	[factor*cos(d*t),factor*sin(d*t),factor*fim];
end:

#savelib('`algcurves/pl_allsub`','`algcurves/pl_curve`','plot_knot'):