with(plots):
with(LinearAlgebra):

EFG:= proc(X)
   local Xu, Xv, E, F, G;
   Xu:= <diff(X[1], u), diff(X[2], u), diff(X[3], u)>;   
   Xv:= <diff(X[1], v), diff(X[2], v), diff(X[3], v)>;
   E:= DotProduct(Xu, Xu, conjugate=false);
   F:= DotProduct(Xu, Xv, conjugate=false);
   G:= DotProduct(Xv, Xv, conjugate=false);
   simplify([E, F, G]);
end:

UN:= proc(X)
   local Xu, Xv, Z, s;
   Xu:= <diff(X[1], u), diff(X[2], u), diff(X[3], u)>;   
   Xv:= <diff(X[1], v), diff(X[2], v), diff(X[3], v)>;
   Z:= CrossProduct(Xu, Xv);
   s:= VectorNorm(Z, Euclidean, conjugate=false);
   simplify(<Z[1]/s|Z[2]/s|Z[3]/s>, sqrt, trig, symbolic);
end:

lmn:= proc(X)
   local Xu, Xv, Xuu, Xuv, Xvv, U, l, m, n;
   Xu:= <diff(X[1], u), diff(X[2], u), diff(X[3], u)>;   
   Xv:= <diff(X[1], v), diff(X[2], v), diff(X[3], v)>;
   Xuu:= <diff(Xu[1], u), diff(Xu[2], u), diff(Xu[3], u)>;   
   Xuv:= <diff(Xu[1], v), diff(Xu[2], v), diff(Xu[3], v)>;   
   Xvv:= <diff(Xv[1], v), diff(Xv[2], v), diff(Xv[3], v)>;   
   U:= UN(X);
   l:= DotProduct(U, Xuu, conjugate=false);
   m:= DotProduct(U, Xuv, conjugate=false);
   n:= DotProduct(U, Xvv, conjugate=false);
   simplify([l,m,n], sqrt, trig, symbolic);
end:

shape:= proc(X)
   local Y, Z, a, b, c, d;
   Y:= EFG(X);
   Z:= lmn(X);
   a:= simplify((Z[1]*Y[3] - Z[2]*Y[2])/(Y[1]*Y[3] - Y[2]^2));
   b:= simplify((Z[2]*Y[1] - Z[1]*Y[2])/(Y[1]*Y[3] - Y[2]^2));
   c:= simplify((Z[2]*Y[3] - Z[3]*Y[2])/(Y[1]*Y[3] - Y[2]^2));
   d:= simplify((Z[3]*Y[1] - Z[2]*Y[2])/(Y[1]*Y[3] - Y[2]^2));
   [S(x_u) = a*x_u + b*x_v, S(x_v) = c*x_u + d*x_v];
end: 

shape_matrix:= proc(X)
   local Y, Z, a, b, c, d;
   Y:= EFG(X);
   Z:= lmn(X);
   a:= simplify((Z[1]*Y[3] - Z[2]*Y[2])/(Y[1]*Y[3] - Y[2]^2));
   b:= simplify((Z[2]*Y[1] - Z[1]*Y[2])/(Y[1]*Y[3] - Y[2]^2));
   c:= simplify((Z[2]*Y[3] - Z[3]*Y[2])/(Y[1]*Y[3] - Y[2]^2));
   d:= simplify((Z[3]*Y[1] - Z[2]*Y[2])/(Y[1]*Y[3] - Y[2]^2));
   Matrix([[a,c], [b,d]]);
end:

normal_curvature:= proc(X, u0, v0, t)
   local S, evals, k;
   S:= subs({u=u0, v=v0}, shape_matrix(X));
   evals:= Eigenvalues(S);
   if evalf(evals[1]) < evalf(evals[2]) then
      k:= [evals[2], evals[1]];
   else
      k:= [evals[1], evals[2]];
   end if;
   k[1]*cos(t)^2 + k[2]*sin(t)^2;
end:

rule_surf:= proc(beta, delta)
   simplify(beta + v*delta, symbolic);
end: