ek-curvature.mac

/* -*- mode: maxima -*- */

/*
 Input:
   (x0,y0) and (x2,y2) - the endpoints
   (px,py)             - the point to interpolate at the maximal curvature
   a                   - alpha, the ratio for the conversion between quadratic & cubic
 Output:
   <a 9th-degree polynomial in t, having one real solution in [0,1]>
 */

/* (cx,cy) is the curve, (dx,dy) and (ddx,ddy) are the 1st and 2nd derivatives */
cx: (1-t)^3*x0+3*(1-t)^2*t*((1-a)*x0+a*x1)+3*(1-t)*t^2*((1-a)*x2+a*x1)+t^3*x2$
cy: (1-t)^3*y0+3*(1-t)^2*t*((1-a)*y0+a*y1)+3*(1-t)*t^2*((1-a)*y2+a*y1)+t^3*y2$
dx: diff(cx,t)$
dy: diff(cy,t)$
ddx: diff(dx,t)$
ddy: diff(dy,t)$

/* n and d are the numerator and denominator of the curvature, respectively */
n: dx*ddy-ddx*dy$
d: (dx^2+dy^2)^(3/2)$

/* The numerator of the curvature's derivative; we need to solve dk = 0 */
dk: diff(n,t)*d-n*diff(d,t)$
/* Looking at factor(dk), we can see that there is some room for simplification */
dk1: factor(dk/(162*a*(x1*y2-x0*y2-x2*y1+x0*y1+x2*y0-x1*y0)*sqrt(dx^2+dy^2)))$
solution: rhs(solve(dk1,t)[1])$

/* (x1,y1) is set s.t. the curve interpolates (px,py) */
x1: (px-((1-t)^3+3*(1-t)^2*t*(1-a))*x0-(3*(1-t)*t^2*(1-a)+t^3)*x2)/(3*(1-t)*t*a)$
y1: (py-((1-t)^3+3*(1-t)^2*t*(1-a))*y0-(3*(1-t)*t^2*(1-a)+t^3)*y2)/(3*(1-t)*t*a)$

/* Generate a string representation that can be inserted in a program */
display2d: false$ /* programming-friendly output */
collectterms(expand(num(xthru(ev(solution,x1=x1,y1=y1)))),t);