with(linalg):
interface(verboseproc=2):

sub := proc(x,y) matadd(x,y,1,-1)
#  sub(x,y) = x - y  (x, y vectors or matrices)
   end:

len := proc(v)
# return the length of a real vector
local i,n,s:
n:=vectdim(v):
s:=0:
for i from 1 to n do
s := s + (v[i]^2):
od:
RETURN(simplify(sqrt(s)))
end:

triple := proc(x,y,z) dotprod(x,crossprod(y,z)):
#  triple(x,y,z) = triple scalar product of vectors in 3-space
   end:

angcos := proc(x,y) simplify((dotprod(x,y))/(len(x)*len(y))):
#  angcos(x,y) = cosine of the angle between x and y
   end:

sproj := proc(x,y) dotprod(x,y)/len(y):
#  sproj(x,y) = scalar projection of vector x on vector y
   end:

proj := proc(x,y) scalarmul(y, dotprod(x,y)/(len(y)^2) ):
#  proj(x,y) = vector projection of vector x on vector y
   end:

perp := proc(x,y) matadd(x, proj(x,y), 1, -1):
#  perp(x,y) = perpendicular projection of vector x on vector y
   end:

dist := proc(x,y) len(sub(y,x)):
#  dist(x,y) = distance from point x to point y
   end:

ldist := proc(z,x,y)
#  ldist(ext. pt., pt. 1 on line, pt. 2 on line)
   local v,w : w:=matadd(x,y,-1,1): v:=matadd(x,z,-1,1):
   len(perp(v,w)): end:

pdist := proc(x,y,v) 
#  pdist(ext. pt., pt. in plane, attitude)
   abs(sproj(sub(x,y),v)) end:

vecsub := proc(M,v)
#  vecsub(M,v) -- substitute value v for "t" in vector M
local i,n,N;
    n := vectdim(M);
    N := vector(n);
    for i to n do  N[i] := simplify(subs(t = v,M[i])) od;
    RETURN(op(N))
end:

div := proc(v,x) simplify(diverge(v,x)) end:
lapl := proc(f,x) simplify(laplacian(f,x)) end:
curl2 := proc(v,p)
  if vectdim(v) <> 2 then
    print("curl2 is for planar vectors only"):
    RETURN()
  elif vectdim(p) <> 2 then 
    print("curl2 is for planar vectors only"):
    RETURN()
  else RETURN(simplify(diff(v[2],p[1])-diff(v[1],p[2])))
  fi:
  end:

# "rsubs" is "subs" with parameters switched (use with "map")
rsubs := proc(x,y) subs(y,x) end:

wint := proc(F,r,a,b)
# wint(F,r,a,b) = path integral of the vector field  F  over the
#         path parameterized by  r(t)  for  t  from  a  to  b
#         Note: 't', 'x', 'y', and 'z' must have global meaning
#         in the session
local V, V1, W, W1, W2, d, j, sl, sp:
d := vectdim(F);
if vectdim(r) <> d then
  print("wint: first two args. must be vectors of the same length");
  RETURN()
fi;
V1 := map(diff,r,t);
# print(V1);
if d = 1 then
  W1 := map(rsubs, F, {x=r[1]})
fi;
if d = 2 then
  W1 := map(rsubs, F, {x=r[1], y=r[2]})
fi;
if d = 3 then
  W1 := map(rsubs, F, {x=r[1], y=r[2], z = r[3]})
fi;
# print(W1);
V  := map(simplify,V1);
W2 := map(simplify, W1);
W  := map(eval, W2);
# print(V, " ", W);
sp := eval(simplify(dotprod(W,V)));
# print(sp);
int(sp,t=a..b)
end: