/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 0.8.4 ] */

/* [wxMaxima: title   start ]
Approximating a Spiral with a Cubic Bezier Curve
   [wxMaxima: title   end   ] */

/* [wxMaxima: comment start ]
Spirals generally do not exist as primitives in maistream graphics APIs. Many systems do have
support for Bezier curves, however. This worksheet shows how to determine the control points of a cubic Bezier curve
that approximates a segment of an Archimedean spiral. The control points must be recomputed for each segment in order
to draw a full spiral. It is not correct to simply reuse a scaled and rotated version of a single segment.
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Definitions
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
We start by defining the equation for a general Bezier curve from Berstein polynomials.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
bernstein(i, n, t) := binomial(n, i)*t^i*(1-t)^(n-i);
bezier(n) := sum(P[i]*bernstein(i,n,t), i, 0, n);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Specifically, we are interested in a cubic Bezier curve.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
C(t):=bezier(3)$
C(t);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Geometrical Configuration
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
The cubic Bezier curve is defined by four points, P[0..3]. P[0] and P[3] are endpoints, while
P[1] and P[2] are their respective control points.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
r(%theta) := 1 + %theta/(2*%pi);
spiral[0](w) := r(w)*cos(w);
spiral[1](w) := r(w)*sin(w);
wxplot2d ([parametric, spiral[0](w), spiral[1](w), [w, -%pi, 2*%pi*2], [nticks, 100]])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Endpoints
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Consider a spiral arc A of sweep < 180°, rising from the x-axis.
Force the endpoints of C to the endpoints of A.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
numsegments : 10;
%theta : 2*%pi/numsegments;
curve_start : 0;
curve_end : curve_start + %theta;
curve_midpoint : (curve_start + curve_end) / 2;

eq0 : x[0] = spiral[0](curve_start);
eq1 : y[0] = spiral[1](curve_start);
eq2 : x[3] = spiral[0](curve_end);
eq3 : y[3] = spiral[1](curve_end);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Midpoint
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
Force the midtime of the Bezier curve to pass through the midpoint of the spiral arc segment.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
s(w) := 1/2 * ( (w * sqrt(1 + w^2)) + log(w + sqrt(1 + w^2)) );
q : (s(curve_end) - s(curve_midpoint))/(s(curve_end) - s(curve_start))$
q : float(q);

eq4 : at(subst(x,P,C(t)),t=q) = spiral[0](curve_midpoint);
eq5 : at(subst(y,P,C(t)),t=q) = spiral[1](curve_midpoint);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Slope
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
We ensure that the spiral's slope at each endpoint is equal to the the slope of the segment connecting each endpoint
to its control point.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
dy : diff(spiral[1](w), w);
dx : diff(spiral[0](w), w);
slope(w) := dy/dx;
cp0_slope : (y[0]-y[1])/(x[0]-x[1]);
cp1_slope : (y[3]-y[2])/(x[3]-x[2]);
eq6 : at(slope(w), w=curve_start) = cp0_slope;
eq7 : at(slope(w), w=curve_end) = cp1_slope;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: subsect start ]
Solution
   [wxMaxima: subsect end   ] */

/* [wxMaxima: comment start ]
We now solve a system of eight variables (X[0..3] and Y[0..3]) and eight equations.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
ans : float(linsolve(
    [eq0, eq1, eq2, eq3, eq4, eq5, eq6, eq7],
    [x[0], y[0], x[1], y[1], x[2], y[2], x[3], y[3]]
));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
dist(dx, dy) := sqrt(dx*dx + dy*dy);
cp0_length : at( dist(x[1]-x[0], y[1]-y[0]), ans );
cp1_length : at( dist(x[2]-x[3], y[2]-y[3]), ans );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
/home/kostmo/hacks/snippets/python/bezier/circular_arc/stink.png
   [wxMaxima: comment end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$