# Frenet Frame Display Package Version 0.8

with(linalg): with(plots):

# Generate the tangent to a curve as a list.
# Typical usage: tang([f(var),g(var),h(var)],var);

tang := proc(gamma,t);
	map(diff,gamma,t)
end:

# A helix paramterized by x.
#sigma := [cos(x),sin(x),x];
#map(diff,sigma,x);

#tang(sigma,x);

# Generate a symbolic unit tangent as a list.
# Typical usage: T([f(var),g(var),h(var)],var);
T := proc(gamma,t)
	local T_len,T_before;
	T_before := vector(3,tang(gamma,t));
	T_len := sqrt(norm(T_before,2));
	evalm((1/T_len) * T_before)
	end:

subs2 := proc(a,b) subs(b,a) end:

# Return a Frenet frame  and curvature at a parameter value t=t0.

# Typical usage: TNB([f(var),g(var),h(var)],var,value).

TNB := proc(gamma,t,t0)
	local gamma_dot,gamma_ddot,gamma_dot_val,gamma_ddot_val,curvature,N_before,T,N,N_len,B,T_len;
	gamma_dot:= vector(3,map(diff,gamma,t));
	gamma_ddot := map(diff,gamma_dot,t);
	gamma_dot_val := map(evalf @ subs2,gamma_dot,t=t0);
	gamma_ddot_val := map(evalf @ subs2,gamma_ddot,t=t0);
	T_len := norm(gamma_dot_val,2);
	N_before:= evalm((1/T_len^3)*(T_len*gamma_ddot_val- dotprod(gamma_dot_val,gamma_ddot_val)));
	curvature := norm(N_before,2);
	T:= evalm((1/T_len)*gamma_dot_val);
	N := evalm((1/curvature)*N_before);
	B := crossprod(T,N);
	map(eval,[T,N,B,curvature])
	end:

#TNB(sigma,x,Pi/2);

# Return a list of num_steps  [T,N,B,curvature] values for parameters in t_range.
# Typical usage: TNB2([f(var),g(var),h(var)],var= val1 .. val2,num_steps,).
TNB2 := proc(gamma,t_range,num_steps)
	local point,gamma_dot,gamma_ddot,gamma_dot_val,gamma_ddot_val,curvature,N_before,T,N,N_len,B,t,t0,t1,t_val,i,T_len,result;
	t := lhs(t_range);
	t0 := op(1,rhs(t_range));
	t1 := op(2,rhs(t_range));
	
	gamma_dot:= vector(3,map(diff,gamma,t));
	gamma_ddot := map(diff,gamma_dot,t);
	
	for i from 1 to num_steps do
		if num_steps > 1 then 
			t_val := t0 + ((t1-t0)*i)/(num_steps-1)
		else
			t_val := t0
		fi;
		point := map(evalf @ subs2,gamma,t=t_val);
		gamma_dot_val := map(evalf @ subs2,gamma_dot,t=t_val);
		gamma_ddot_val := map(evalf @ subs2,gamma_ddot,t=t_val);
		T_len := norm(gamma_dot_val,2);
		N_before:= evalm((1/T_len^3)*(T_len*gamma_ddot_val- dotprod(gamma_dot_val,gamma_ddot_val)));
		curvature := norm(N_before,2);
		T:= evalm((1/T_len)*gamma_dot_val);
		N := evalm((1/curvature)*N_before);
		B := crossprod(T,N);
		result[i] := map(eval,[point,T,N,B,curvature]);
	od;
	seq(result[i],i=1..num_steps)
	end:

#TNB2(sigma,x=Pi/2 .. Pi,2);

# Return a graph of  num_steps Frenet frames, equally spaced along the curve gamma.
# Each frame is drawn with T in red, N in blue, and B in green.
# The length of each vector is vec_len.
# Typical usage: TNB2([f(var),g(var),h(var)],var= val1 .. val2,num_steps,vec_len).
TNB3 := proc(gamma,t_range,num_steps,vec_len)
	local point,gamma_dot,gamma_ddot,gamma_dot_val,gamma_ddot_val,curvature,N_before,T,N,N_len,B,t,t0,t1,t_val,i,T_len,result,gr,T_end,N_end,B_end;
	t := lhs(t_range);
	t0 := op(1,rhs(t_range));
	t1 := op(2,rhs(t_range));
	
	gamma_dot:= vector(3,map(diff,gamma,t));
	gamma_ddot := map(diff,gamma_dot,t);
	
	for i from 1 to num_steps do
		if num_steps > 1 then 
			t_val := t0 + ((t1-t0)*(i-1))/(num_steps-1)
		else
			t_val := t0
		fi;
		point := map(evalf @ subs2,gamma,t=t_val);
		gamma_dot_val := map(evalf @ subs2,gamma_dot,t=t_val);
		gamma_ddot_val := map(evalf @ subs2,gamma_ddot,t=t_val);
		T_len := norm(gamma_dot_val,2);
		N_before:= evalm((1/T_len^3)*(T_len*gamma_ddot_val- dotprod(gamma_dot_val,gamma_ddot_val)));
		curvature := norm(N_before,2);
		T:= evalm((1/T_len)*gamma_dot_val);
		N := evalm((1/curvature)*N_before);
		B := crossprod(T,N);
		result[i] := map(eval,[point,T,N,B,curvature]);
		T_end := convert(evalm(point + vec_len *T),list);
		N_end := convert(evalm(point + vec_len* N),list);
		B_end := convert(evalm(point + vec_len*B),list);
		gr[i] := polygonplot3d([point,T_end],color=red,thickness=3),polygonplot3d([point,N_end],color=green,thickness=3),polygonplot3d([point,B_end],color=blue,thickness=3);
	od;
	{seq(gr[i],i=1..num_steps)}
	end:
#frames := TNB3(sigma,x=0 .. Pi,4,.3):
#curve_gr := spacecurve(sigma,x= 0 ..Pi,thickness =3,axes=framed):
# display3d(frames union {curve_gr});