/* bezier-circle.c
   by Mike Bertrand
   Mar 16, 2015

   Calculates how close an approximating Bezier curve is to a quarter circle.
   Consider the unit circle in the first quadrant and point a horizontal vector
   to the right from (0,1) and another vertical vector up from (1,0). The Bezier
   curve handles are a = 4/3*(sqrt(2) - 1) ~ 0.5523 units along those vectors.

   Sample at equally spaced intervals for 0 < t < 1. Calculate the corresponding
   point on the Bezier curve, connect it to the origin by a straight line, then
   see where that line intersects the unit circle. Then calculate the distance
   between the corresponing points on the Bezier curve and the circle.

   The maximum distance occurrs at about t = 0.21, angle apprimately 19 degrees
   (and 71 degrees, the situation is symmetric around the line at 45 degrees).

   The maximum distance is about 0.00027. We're dealing with a unit circle, so
   that's the maximum divergence, less than 3 parts in 10,000.
 */

#include "stdio.h"
#include "math.h"     // sqrt, hypot
#include "stdlib.h"   // fabs

// Number of equally spaced values in [0, 1] to sample.
#define N 50

double calcBezCoord(double p[], double t);

void main(void)
{
  int k = -1;

  // Parameter in [0, 1].
  double t = -1.0;

  // Points on the Bezier curve and circle respectively.
  double xBez, yBez;
  double xCircle, yCircle;

  // Length of Bezier curve handles when approximating quarter circle.
  double a = (double)4/3*(sqrt(2) - 1);

  // Bezier curve x and y coordinates respectively.
  double px[] = {0.0, a, 1.0, 1.0};
  double py[] = {1.0, 1.0, a, 0.0};

  // Slope of line to origin and sqrt(1 + slope^2).
  double m = -1.0;
  double n = -1;

  // Deltas between Bezier curve and circle.
  double dx = -1;
  double dy = -1;
  double distance = -1;

  printf("    t      xBez     yBez   xCircle  yCircle  distance\n");
  printf("-----------------------------------------------------\n");

  for (k = 1; k < N; k++)
  {
	// Space t at equal intervals of size 1/N, calculate corresponding point on Bezier.
    t = (double) k / N;
    xBez = calcBezCoord(px, t);
    yBez = calcBezCoord(py, t);

    // Connect point (xBez, yBez) to origin - m is slope of that line.
    m = yBez / xBez;
    n = sqrt(1 + m*m);

    // (xCircle, yCircle) is intersection of the line and unit circle.
    xCircle = 1 / n;
    yCircle = m / n;

    // dx and dy are deltas between Bezier point and circle point, distance = how far appart they are.
    dx = xBez - xCircle;
    dy = yBez - yCircle;
    distance = hypot(dx, dy);

    printf("%7.5f  %7.5f  %7.5f  %7.5f  %7.5f  %7.5f\n", t, xBez, yBez, xCircle, yCircle, distance);
  }
}


double calcBezCoord(double p[], double t)
{
  double s = 1.0 - t;
  return p[0]*s*s*s + 3*p[1]*s*s*t + 3*p[2]*s*t*t + p[3]*t*t*t;
}