# bezcurve.m  - commands to compute bezier curves given matrix of
# control pts
printf("\nwith(plots): with(plottools): with(linalg):");
with(plots): with(plottools): with(linalg):

printf("\nCommands in bezcurve.txt");
printf("\n  BezCurve(matrix pts) - computes Bezier curve from matrix of control points");
printf("\n  PolyLine(matrix pts,color) - draws control polyline in the color from matrix of control points");
printf("\n  BezPlot(curve,color) - plots Bezier curve in the color given without axes");
printf("\n  BezSpline0(matrix pts,color) - plots piecewise Bezier curve from control pts");
printf("\n  BezSpline1(matrix pts,color) - plots C1 cubic Bezier spline from control pts"); 
printf("\n  BezSpline2(matrix pts,color) - plots C2 cubic Bezier spline from control pts");
printf("\n  BezContPolySpline1(matric pts) - plots control polyline in blue for C1 Bezier Spline");
printf("\n  BezContPolySpline2(matrix pts) - plots control polyline in red for C2 Bezier Spline");
printf("\n  BezierCurvature(matrix pts) - creates Curvature plot for Bezier Curve");
printf("\n  Spline1Curvature(matrix pts) - creates Curvature plot for C1 Bezier Spline");
printf("\n  Spline2Curvature(matrix pts) - creates Curvature plot for C2 Bezier Spline");
printf("\n  CurvDist1(curve1,curve2,n) - approximates the closeness of two curves based on distance");
printf("\n  CurvDist2(curve1,curve2,n) - approximates the closeness of two curves based on distance and tangents");

# bezcurv - computes bezier curve given matrix of control pts
BezCurve:=proc(pts)
  local d, i, m, xp, yp, x, y, ans;
   m:=rowdim(pts):
   d:=m-1:
   for i from 0 to d do
     xp[i]:=pts[i+1,1]:
     yp[i]:=pts[i+1,2]:
   od:
   unassign('i'):
   x:=evalm(sum(binomial(d,i)*s^i*(1-s)^(d-i)*xp[i],i=0..d)):
   y:=evalm(sum(binomial(d,i)*s^i*(1-s)^(d-i)*yp[i],i=0..d)):
   ans:=vector([x,y]):
   RETURN (evalm(ans));
  end:
  
# plots control polygon of bezier curve given matrix of control pts
PolyLine:=proc(pts,col)
  local a, b, i, m, xp, yp, L, ans;
  m:=rowdim(pts):
  for i from 1 to m do
    xp[i]:=pts[i,1]:
    yp[i]:=pts[i,2]:
  od:
  unassign('i'):
  for i from 1 to m-1 do
   
   L[i]:=line([xp[i],yp[i]],[xp[i+1],yp[i+1]],color=col,linestyle=2):
   od:
   unassign('i'):
  
  a:=display([seq(L[i],i=1..m-1)],insequence=false,scaling=constrained):
  b:=pointplot(pts,color=black):
  ans:=display({a,b},axes=none):
  RETURN (ans); 
  end:
  
BezPlot:=proc(curve,col)
 plot([curve[1],curve[2],s=0..1],scaling=constrained,axes=none,color=col);
 end: 

# procedure to draw a piecewise Bezier curve from control points
BezSpline0:=proc(pts,col)
  local i,m,n,px,py,px0,py0,px1,py1,px2,py2,px3,py3,c,P:
  m:=rowdim(pts):
  n:=(m-1)/3:
  for i from 0 to m-1 do
    px[i]:=pts[i+1,1]:
    py[i]:=pts[i+1,2]:
  od:
  unassign('i'):

  for i from 1 to n do
    px0:=px[3*i-3]: py0:=py[3*i-3]:
    px1:=px[3*i-2]: py1:=py[3*i-2]:
    px2:=px[3*i-1]: py2:=py[3*i-1]:
    px3:=px[3*i]:    py3:=py[3*i]:
    c:=BezCurve(matrix([[px0,py0],[px1,py1],[px2,py2],[px3,py3]])):
    P[i]:=BezPlot(c,col):
  od:
  RETURN (display([seq(P[i],i=1..n)],insequence=false,scaling=constrained,axes = none));
  end:
   
# procedure to draw a Bezier composite curve from spline control points   
BezSpline1:=proc(pts,col)
  local i,m,n,px,py,dx,dy,px0,py0,px1,py1,px2,py2,px3,py3,P,c:
  m:=rowdim(pts):
  n:=(m-2)/2:
  # loop to set control points for 
  for i from -1 to 2*n do
    dx[i]:=pts[i+2,1]:
    dy[i]:=pts[i+2,2]:
  od:
  unassign('i'):
   
  px[0]:=dx[-1]: py[0]:=dy[-1]:
  px[1]:=dx[0]: py[1]:=dy[0]:
  px[2]:=dx[1]: py[2]:=dy[1]:
  
  for i from 1 to n-1 do
    px[3*i+1]:=dx[2*i]:
    py[3*i+1]:=dy[2*i]:
    px[3*i+2]:=dx[2*i+1]:
    py[3*i+2]:=dy[2*i+1]:
  od:
  unassign('i'):
   
  px[3*n]:=dx[2*n]: py[3*n]:=dy[2*n]:
  
  for i from 1 to n-1 do
    px[3*i]:= (1/2)*px[3*i-1]+(1/2)*px[3*i+1]: 
    py[3*i]:= (1/2)*py[3*i-1]+(1/2)*py[3*i+1]: 
  od:
  unassign('i'):

  for i from 1 to n do   
    px0:=px[3*i-3]: py0:=py[3*i-3]:
    px1:=px[3*i-2]: py1:=py[3*i-2]:
    px2:=px[3*i-1]: py2:=py[3*i-1]:
    px3:=px[3*i]:   py3:=py[3*i]:
  
    c:=BezCurve(matrix([[px0,py0],[px1,py1],[px2,py2],[px3,py3]])):
    P[i]:=BezPlot(c,col):
  od:
   
  RETURN (display([seq(P[i],i=1..n)],insequence=false,scaling=constrained,axes = none));
  end:


# procedure to draw a Bezier composite curve from spline control points   
BezSpline2:=proc(pts,col)
  local i,m,L,px,py,dx,dy,px0,py0,px1,py1,px2,py2,px3,py3,P,c:
  m:=rowdim(pts):
  L:=m-3:
  # loop to set control points for 
  for i from -1 to L+1 do
    dx[i]:=pts[i+2,1]:
    dy[i]:=pts[i+2,2]:
  od:
  unassign('i'):
   
  px[0]:=dx[-1]: py[0]:=dy[-1]:
  px[1]:=dx[0]: py[1]:=dy[0]:
  px[2]:=(1/2)*dx[0]+(1/2)*dx[1]: py[2]:=(1/2)*dy[0]+(1/2)*dy[1]:
  
  for i from 1 to L-2 do
    px[3*i+1]:=(2/3)*dx[i]+(1/3)*dx[i+1]:
    py[3*i+1]:=(2/3)*dy[i]+(1/3)*dy[i+1]:
    px[3*i+2]:=(1/3)*dx[i]+(2/3)*dx[i+1]:
    py[3*i+2]:=(1/3)*dy[i]+(2/3)*dy[i+1]:
  od:
  unassign('i'):
   
  px[3*L-2]:=(1/2)*dx[L-1]+(1/2)*dx[L]: py[3*L-2]:=(1/2)*dy[L-1]+(1/2)*dy[L]:
  px[3*L-1]:=dx[L]: py[3*L-1]:=dy[L]:
  px[3*L]:=dx[L+1]: py[3*L]:=dy[L+1]:

  for i from 1 to L-1 do
    px[3*i]:= (1/2)*px[3*i-1]+(1/2)*px[3*i+1]: 
    py[3*i]:= (1/2)*py[3*i-1]+(1/2)*py[3*i+1]: 
  od:
  unassign('i'):

  for i from 1 to L do   
    px0:=px[3*i-3]: py0:=py[3*i-3]:
    px1:=px[3*i-2]: py1:=py[3*i-2]:
    px2:=px[3*i-1]: py2:=py[3*i-1]:
    px3:=px[3*i]:   py3:=py[3*i]:
  
    c:=BezCurve(matrix([[px0,py0],[px1,py1],[px2,py2],[px3,py3]])):
    P[i]:=BezPlot(c,col):
  od:
   
  RETURN (display([seq(P[i],i=1..L)],insequence=false,scaling=constrained,axes = none));
  end:


# procedure to compute and plot control points for composite curve from spline control points   
BezContPolySpline1:=proc(pts)
  local i,m,n,px,py,dx,dy,px0,py0,px1,py1,px2,py2,px3,py3,P,c:
 
  m:=rowdim(pts):
  n:=(m-2)/2:
  # loop to set control points for 
  for i from -1 to 2*n do
    dx[i]:=pts[i+2,1]:
    dy[i]:=pts[i+2,2]:
  od:
  unassign('i'):
   
  px[0]:=dx[-1]: py[0]:=dy[-1]:
  px[1]:=dx[0]: py[1]:=dy[0]:
  px[2]:=dx[1]: py[2]:=dy[1]:
  
  for i from 1 to n-1 do
    px[3*i+1]:=dx[2*i]:
    py[3*i+1]:=dy[2*i]:
    px[3*i+2]:=dx[2*i+1]:
    py[3*i+2]:=dy[2*i+1]:
  od:
  unassign('i'):
   
  px[3*n]:=dx[2*n]: py[3*n]:=dy[2*n]:
  
  for i from 1 to n-1 do
    px[3*i]:= (1/2)*px[3*i-1]+(1/2)*px[3*i+1]: 
    py[3*i]:= (1/2)*py[3*i-1]+(1/2)*py[3*i+1]: 
  od:
  unassign('i'):

  for i from 1 to n do   
    px0:=px[3*i-3]: py0:=py[3*i-3]:
    px1:=px[3*i-2]: py1:=py[3*i-2]:
    px2:=px[3*i-1]: py2:=py[3*i-1]:
    px3:=px[3*i]:   py3:=py[3*i]:
  
    P[i]:=PolyLine(matrix([[px0,py0],[px1,py1],[px2,py2],[px3,py3]]),blue):
  od:
   
  RETURN (display([seq(P[i],i=1..n)],insequence=false,scaling=constrained,axes = none));
  end:

# procedure to compute and plot control points for composite curve from spline control points   
BezContPolySpline2:=proc(pts)
  local i,m,L,px,py,dx,dy,px0,py0,px1,py1,px2,py2,px3,py3,P,c:
  m:=rowdim(pts):
  L:=m-3:
  # loop to set control points for 
  for i from -1 to L+1 do
    dx[i]:=pts[i+2,1]:
    dy[i]:=pts[i+2,2]:
  od:
  unassign('i'):
   
  px[0]:=dx[-1]: py[0]:=dy[-1]:
  px[1]:=dx[0]: py[1]:=dy[0]:
  px[2]:=(1/2)*dx[0]+(1/2)*dx[1]: py[2]:=(1/2)*dy[0]+(1/2)*dy[1]:
  
  for i from 1 to L-2 do
    px[3*i+1]:=(2/3)*dx[i]+(1/3)*dx[i+1]:
    py[3*i+1]:=(2/3)*dy[i]+(1/3)*dy[i+1]:
    px[3*i+2]:=(1/3)*dx[i]+(2/3)*dx[i+1]:
    py[3*i+2]:=(1/3)*dy[i]+(2/3)*dy[i+1]:
  od:
  unassign('i'):
 
  px[3*L-2]:=(1/2)*dx[L-1]+(1/2)*dx[L]:
  py[3*L-2]:=(1/2)*dy[L-1]+(1/2)*dy[L]:
  px[3*L-1]:=dx[L]: py[3*L-1]:=dy[L]:
  px[3*L]:=dx[L+1]: py[3*L]:=dy[L+1]:
 
  for i from 1 to L-1 do
    px[3*i]:=(1/2)*px[3*i-1]+(1/2)*px[3*i+1]:
    py[3*i]:=(1/2)*py[3*i-1]+(1/2)*py[3*i+1]:
  od:
  unassign('i'): 
 
  for i from 1 to L do   
    px0:=px[3*i-3]: py0:=py[3*i-3]:
    px1:=px[3*i-2]: py1:=py[3*i-2]:
    px2:=px[3*i-1]: py2:=py[3*i-1]:
    px3:=px[3*i]:   py3:=py[3*i]:
  
    P[i]:=PolyLine(matrix([[px0,py0],[px1,py1],[px2,py2],[px3,py3]]),red):
  od:
   
  RETURN (display([seq(P[i],i=1..L)],insequence=false,scaling=constrained,axes = none));
  end:

# procedure to create a curvature plot for a Bezier curve
BezierCurvature:=proc(pts)
  local i,c,v,speed,a,vxa,curv,t,arc,K:
  c:=BezCurve(pts):
  v:=vector([diff(c[1],s),diff(c[2],s),0]):
  speed:=sqrt(v[1]^2+v[2]^2):
  a:=vector([diff(v[1],s),diff(v[2],s),0]):
  vxa:=crossprod(v,a):
  curv:=vxa[3]/(sqrt(v[1]^2+v[2]^2))^3:
  for i from 0 to 100 do
    t:=i*0.01:
    arc[i]:=evalf(Int(speed,s=0..t)):
    K[i]:=evalf(subs(s=t,curv)):
  od:
  RETURN (display(curve([seq([arc[i],K[i]],i=0..100)],color=red,scaling=constrained)));
  end:

# procedure to create a curvature plot for a Bezier Spline
Spline1Curvature:=proc(pts)
  local i,j,m,n,px,py,dx,dy,px0,py0,px1,py1,px2,py2,px3,py3,P,c,v,speed,a,vxa,curv,arc,K,lngth:

  m:=rowdim(pts):
  n:=(m-2)/2:
  # loop to set control points for 
  for i from -1 to 2*n do
    dx[i]:=pts[i+2,1]:
    dy[i]:=pts[i+2,2]:
  od:
  unassign('i'):
   
  px[0]:=dx[-1]: py[0]:=dy[-1]:
  px[1]:=dx[0]: py[1]:=dy[0]:
  px[2]:=dx[1]: py[2]:=dy[1]:
  
  for i from 1 to n-1 do
    px[3*i+1]:=dx[2*i]:
    py[3*i+1]:=dy[2*i]:
    px[3*i+2]:=dx[2*i+1]:
    py[3*i+2]:=dy[2*i+1]:
  od:
  unassign('i'):
   
  px[3*n]:=dx[2*n]: py[3*n]:=dy[2*n]:
  
  for i from 1 to n-1 do
    px[3*i]:= (1/2)*px[3*i-1]+(1/2)*px[3*i+1]: 
    py[3*i]:= (1/2)*py[3*i-1]+(1/2)*py[3*i+1]: 
  od:
  unassign('i'):

    lngth:=0.0:
  for i from 1 to n do   
    px0:=px[3*i-3]: py0:=py[3*i-3]:
    px1:=px[3*i-2]: py1:=py[3*i-2]:
    px2:=px[3*i-1]: py2:=py[3*i-1]:
    px3:=px[3*i]:   py3:=py[3*i]:
  
    c:=BezCurve(matrix([[px0,py0],[px1,py1],[px2,py2],[px3,py3]])):
    v:=vector([diff(c[1],s),diff(c[2],s),0]):
    speed:=sqrt(v[1]^2+v[2]^2):
    a:=vector([diff(v[1],s),diff(v[2],s),0]):
    vxa:=crossprod(v,a):
    curv:=vxa[3]/(speed)^3:
    for j from 0 to 100 do
      arc[j]:=lngth+evalf(Int(speed,s=0..(j*0.01))):
      K[j]:=evalf(subs(s=j*0.01,curv)):
    od:
    unassign('j'):
    lngth:=lngth+evalf(Int(speed,s=0..1)):
    P[i]:=display(curve([seq([arc[j],K[j]],j=0..100)],color=red)):
  od:

  RETURN (display([seq(P[i],i=1..n)],insequence=false));
  end:

# procedure to create a curvature plot for a Bezier Spline
Spline2Curvature:=proc(pts)
  local i,j,m,L,px,py,dx,dy,px0,py0,px1,py1,px2,py2,px3,py3,P,c,v,speed,a,vxa,curv,arc,K,lngth:
  m:=rowdim(pts):
  L:=m-3:
  # loop to set control points for 
  for i from -1 to L+1 do
    dx[i]:=pts[i+2,1]:
    dy[i]:=pts[i+2,2]:
  od:
  unassign('i'):
   
  px[0]:=dx[-1]: py[0]:=dy[-1]:
  px[1]:=dx[0]: py[1]:=dy[0]:
  px[2]:=(1/2)*dx[0]+(1/2)*dx[1]: py[2]:=(1/2)*dy[0]+(1/2)*dy[1]:
  
  for i from 1 to L-2 do
    px[3*i+1]:=(2/3)*dx[i]+(1/3)*dx[i+1]:
    py[3*i+1]:=(2/3)*dy[i]+(1/3)*dy[i+1]:
    px[3*i+2]:=(1/3)*dx[i]+(2/3)*dx[i+1]:
    py[3*i+2]:=(1/3)*dy[i]+(2/3)*dy[i+1]:
  od:
  unassign('i'):
   
  px[3*L-2]:=(1/2)*dx[L-1]+(1/2)*dx[L]: py[3*L-2]:=(1/2)*dy[L-1]+(1/2)*dy[L]:
  px[3*L-1]:=dx[L]: py[3*L-1]:=dy[L]:
  px[3*L]:=dx[L+1]: py[3*L]:=dy[L+1]:

  for i from 1 to L-1 do
    px[3*i]:= (1/2)*px[3*i-1]+(1/2)*px[3*i+1]: 
    py[3*i]:= (1/2)*py[3*i-1]+(1/2)*py[3*i+1]: 
  od:
  unassign('i'):

    lngth:=0.0:
  for i from 1 to L do   
    px0:=px[3*i-3]: py0:=py[3*i-3]:
    px1:=px[3*i-2]: py1:=py[3*i-2]:
    px2:=px[3*i-1]: py2:=py[3*i-1]:
    px3:=px[3*i]:   py3:=py[3*i]:
  
    c:=BezCurve(matrix([[px0,py0],[px1,py1],[px2,py2],[px3,py3]])):
    v:=vector([diff(c[1],s),diff(c[2],s),0]):
    speed:=sqrt(v[1]^2+v[2]^2):
    a:=vector([diff(v[1],s),diff(v[2],s),0]):
    vxa:=crossprod(v,a):
    curv:=vxa[3]/(speed)^3:
    for j from 0 to 100 do
      arc[j]:=lngth+evalf(Int(speed,s=0..(j*0.01))):
      K[j]:=evalf(subs(s=j*0.01,curv)):
    od:
    unassign('j'):
    lngth:=lngth+evalf(Int(speed,s=0..1)):
    P[i]:=display(curve([seq([arc[j],K[j]],j=0..100)],color=red)):
  od:

  RETURN (display([seq(P[i],i=1..L)],insequence=false));
  end:


GSpline1:=proc(pts,kn,col)
  local i,t,m,n,px,py,dx,dy,xb,yb,x,y,P,c:
  m:=rowdim(pts):
  n:=(m-2)/2:
  # loop to set control points for 
  for i from -1 to 2*n do
    dx[i]:=pts[i+2,1]:
    dy[i]:=pts[i+2,2]:
  od:
  unassign('i'):
   
  px[0]:=dx[-1]: py[0]:=dy[-1]:
  px[1]:=dx[0]: py[1]:=dy[0]:
  px[2]:=dx[1]: py[2]:=dy[1]:
  
  for i from 1 to n-1 do
    px[3*i+1]:=dx[2*i]:
    py[3*i+1]:=dy[2*i]:
    px[3*i+2]:=dx[2*i+1]:
    py[3*i+2]:=dy[2*i+1]:
  od:
  unassign('i'):
   
  px[3*n]:=dx[2*n]: py[3*n]:=dy[2*n]:
  
  for i from 1 to n-1 do
    px[3*i]:= (kn[i+1]-kn[i])/(kn[i+2]-kn[i])*px[3*i-1]+(kn[i+2]-kn[i+1])/(kn[i+2]-kn[i])*px[3*i+1]: 
    py[3*i]:= (kn[i+1]-kn[i])/(kn[i+2]-kn[i])*py[3*i-1]+(kn[i+2]-kn[i+1])/(kn[i+2]-kn[i])*py[3*i+1]:  
  od:
  unassign('i'):

  for i from 1 to n do   
    xb[0]:=px[3*i-3]: yb[0]:=py[3*i-3]:
    xb[1]:=px[3*i-2]: yb[1]:=py[3*i-2]:
    xb[2]:=px[3*i-1]: yb[2]:=py[3*i-1]:
    xb[3]:=px[3*i]:   yb[3]:=py[3*i]:

    t:=(s-knots[i])/(knots[i+1]-knots[i]):
    x:=evalm(sum(binomial(3,j)*t^j*(1-t)^(3-j)*xb[j],j=0..3)):
    y:=evalm(sum(binomial(3,j)*t^j*(1-t)^(3-j)*yb[j],j=0..3)):
  
    c:=[x,y]:
    P[i]:=plot([c[1],c[2],s=knots[i]..knots[i+1]],color=col,scaling=constrained,axes=none):
  od:
   
  RETURN (display([seq(P[i],i=1..n)],insequence=false,scaling=constrained,axes = none));
  end:


GSpline1Poly:=proc(pts,kn,col)
  local i,t,m,n,px,py,dx,dy,xb,yb,x,y,P,c:
  m:=rowdim(pts):
  n:=(m-2)/2:
  # loop to set control points for 
  for i from -1 to 2*n do
    dx[i]:=pts[i+2,1]:
    dy[i]:=pts[i+2,2]:
  od:
  unassign('i'):
   
  px[0]:=dx[-1]: py[0]:=dy[-1]:
  px[1]:=dx[0]: py[1]:=dy[0]:
  px[2]:=dx[1]: py[2]:=dy[1]:
  
  for i from 1 to n-1 do
    px[3*i+1]:=dx[2*i]:
    py[3*i+1]:=dy[2*i]:
    px[3*i+2]:=dx[2*i+1]:
    py[3*i+2]:=dy[2*i+1]:
  od:
  unassign('i'):
   
  px[3*n]:=dx[2*n]: py[3*n]:=dy[2*n]:
  
  for i from 1 to n-1 do
    px[3*i]:= (kn[i+1]-kn[i])/(kn[i+2]-kn[i])*px[3*i-1]+(kn[i+2]-kn[i+1])/(kn[i+2]-kn[i])*px[3*i+1]: 
    py[3*i]:= (kn[i+1]-kn[i])/(kn[i+2]-kn[i])*py[3*i-1]+(kn[i+2]-kn[i+1])/(kn[i+2]-kn[i])*py[3*i+1]:  
  od:
  unassign('i'):

  for i from 1 to n do   
    xb[0]:=px[3*i-3]: yb[0]:=py[3*i-3]:
    xb[1]:=px[3*i-2]: yb[1]:=py[3*i-2]:
    xb[2]:=px[3*i-1]: yb[2]:=py[3*i-1]:
    xb[3]:=px[3*i]:   yb[3]:=py[3*i]:


    P[i]:=PolyLine(matrix([[xb[0],yb[0]],[xb[1],yb[1]],[xb[2],yb[2]],[xb[3],yb[3]]]),blue):
  od:
   
  RETURN (display([seq(P[i],i=1..n)],insequence=false,scaling=constrained,axes = none));
  end:

# determines the closeness of two curves with respect to distance n points
CurvDist1:=proc(curve1,curve2,n)
  local i, j, c1, c2, t, d, dst, ans;
  for i from 0 to n do
    t:=i/n:
    c1[i]:=evalf(subs(s=t,[curve1[1],curve1[2]])):
    c2[i]:=evalf(subs(s=t,[curve2[1],curve2[2]])):
  od:
  for i from 0 to n do
    for j from 0 to n do
      d[i,j]:=sqrt((c1[i][1]-c2[j][1])^2+(c1[i][2]-c2[j][2])^2):
    od:
  od:
  for i from 0 to n do
    ans:=d[i,i]:
    for j from 0 to n do
      ans:=min(d[i,j],ans):
    od:
       dst[i]:=ans:
  od:
   ans:=dst[0]:
   for i from 1 to n do
      ans:=max(ans,dst[i]):
   od:
   RETURN(ans)
  end:

# determines the closeness of two curves with respect to distance and tangents
CurvDist2:=proc(curve1,curve2,n)
  local i, j, c1, c2, dc1, dc2, T1, T2, t, d, dst, ans;
  for i from 0 to n do
    t:=i/n:
    c1[i]:=evalf(subs(s=t,[curve1[1],curve1[2]])):
    c2[i]:=evalf(subs(s=t,[curve2[1],curve2[2]])):
    dc1:=evalf(subs(s=t,[diff(curve1[1],s),diff(curve1[2],s)])):
    T1[i]:=evalf([dc1[1]/sqrt(dc1[1]^2+dc1[2]^2),dc1[2]/sqrt(dc1[1]^2+dc1[2]^2)]):
    dc2:=evalf(subs(s=t,[diff(curve2[1],s),diff(curve2[2],s)])):
    T2[i]:=evalf([dc2[1]/sqrt(dc2[1]^2+dc2[2]^2),dc2[2]/sqrt(dc2[1]^2+dc2[2]^2)]):
  od:
  for i from 0 to n do
    for j from 0 to n do
      d[i,j]:=sqrt((c1[i][1]-c2[j][1])^2+(c1[i][2]-c2[j][2])^2 + (T1[i][1]-T2[j][1])^2+(T1[i][2]-T2[j][2])^2):
    od:
  od:
  for i from 0 to n do
    ans:=d[i,i]:
    for j from 0 to n do
      ans:=min(d[i,j],ans):
    od:
       dst[i]:=ans:
  od:
  ans:=dst[0]:
  for i from 1 to n do
    ans:=max(ans,dst[i]):
  od:
  RETURN(ans)
  end: