// I wrote one of these from scratch once starting with a 2d spline procedure,
// but I seem to have lost that source.  This version appears to have come from
//                   http://paulbourke.net/geometry/spline/
// - I didn't have the patience to re-implement it from first principles I guess.
// FOUND IT:
#ifdef NEVER
/* The procedure below is based on */
/* http://stackoverflow.com/questions/5443653/opengl-coordinates-from-bezier-curves */
/* - it's a 1D Bezier than can be called twice to create a 2D Bezier spline and */
/* 3 times to create a 3D Bezier spline etc.  Each axis is independent. */

double bezier(double A,  /* Start value */
              double B,  /* First control value */
              double C,  /* Second control value */
              double D,  /* Ending value */
              int t, int max_slots, /* eg 0..255/256 slots */
              double *delta)  /* Parameter 0 <= t <= 1 */
{
  double s = (max_slots-1) - t;
  double AB = (A*s + B*t)/max_slots;
  double BC = (B*s + C*t)/max_slots;
  double CD = (C*s + D*t)/max_slots;
  //  double ABC = (AB*s + CD*t)/max_slots; - this 'bug fix' was wrong...
  double ABC = (AB*s + BC*t)/max_slots; // original code was correct...
  double BCD = (BC*s + CD*t)/max_slots;
  if (delta) *delta = ABC-BCD;
  return (ABC*s + BCD*t)/max_slots;
}
#endif
//---------------------------------------------------- FLIGHTPATH GENERATOR
#include <stdio.h>
// once code is performing properly, replace this with a version using fixed-
// point arithmetic (which is currently broken :-( )

typedef struct XYZ {
  double x,y,z;
} XYZ;

double SplineBlend(int k,int t,int *u,double v);
/*
   This returns the point "output" on the spline curve.
   The parameter "v" indicates the position, it ranges from 0 to n-t+2
   
*/
void SplinePoint(int *u,int n,int t,double v,XYZ *control,XYZ *output)
{
   int k;
   double b;

   output->x = 0;
   output->y = 0;
   output->z = 0;

   for (k=0;k<=n;k++) {
      b = SplineBlend(k,t,u,v);
      output->x += control[k].x * b;
      output->y += control[k].y * b;
      output->z += control[k].z * b;
   }
}

/*
   Calculate the blending value, this is done recursively.
   
   If the numerator and denominator are 0 the expression is 0.
   If the deonimator is 0 the expression is 0
*/
double SplineBlend(int k,int t,int *u,double v)
{
   double value;

   if (t == 1) {
      if ((u[k] <= v) && (v < u[k+1]))
         value = 1;
      else
         value = 0;
   } else {
      if ((u[k+t-1] == u[k]) && (u[k+t] == u[k+1]))
         value = 0;
      else if (u[k+t-1] == u[k]) 
         value = (u[k+t] - v) / (u[k+t] - u[k+1]) * SplineBlend(k+1,t-1,u,v);
      else if (u[k+t] == u[k+1])
         value = (v - u[k]) / (u[k+t-1] - u[k]) * SplineBlend(k,t-1,u,v);
     else
         value = (v - u[k]) / (u[k+t-1] - u[k]) * SplineBlend(k,t-1,u,v) + 
                 (u[k+t] - v) / (u[k+t] - u[k+1]) * SplineBlend(k+1,t-1,u,v);
   }
   return(value);
}

/*
   The positions of the subintervals of v and breakpoints, the position
   on the curve are called knots. Breakpoints can be uniformly defined
   by setting u[j] = j, a more useful series of breakpoints are defined
   by the function below. This set of breakpoints localises changes to
   the vicinity of the control point being modified.
*/
void SplineKnots(int *u,int n,int t)
{
   int j;

   for (j=0;j<=n+t;j++) {
      if (j < t)
         u[j] = 0;
      else if (j <= n)
         u[j] = j - t + 1;
      else if (j > n)
         u[j] = n - t + 2;	
   }
}

/*-------------------------------------------------------------------------
   Create all the points along a spline curve
   Control points "inp", "n" of them.
   Knots "knots", degree "t".
   Output curve "outp", "res" of them.
*/
void SplineCurve(XYZ *inp,int n,int *knots,int t,XYZ *outp,int res)
{
   int i;
   double interval,increment;

   interval = 0;
   increment = (n - t + 2) / (double)(res - 1);
   for (i=0;i<res-1;i++) {
      SplinePoint(knots,n,t,interval,inp,&(outp[i]));
      interval += increment;
   }
   outp[res-1] = inp[n];
}

/*
   Example of how to call the spline functions
	Basically one needs to create the control points, then compute
   the knot positions, then calculate points along the curve.
*/

#define SCL 2048   // x/y fudge factor
#define ZSCL 512   // distance away
#define PSCL 1024  // perspective scale

#define MAXN 10
XYZ inp[] = {
//X        Y        Z
  450.0, 1000.0, 10.0,        // final coord is center of 1280x1024 screen
  450.0, 1000.0, 6.66,        // final coord is center of 1280x1024 screen
  450.0, 1000.0, 3.33,        // final coord is center of 1280x1024 screen
  450.0, 1000.0, 0.0,        // final coord is center of 1280x1024 screen
};

#define T 3
int knots[MAXN+T+1];
#define RESOLUTION 512
XYZ outp[RESOLUTION];

void newpaths(void) // yet another quick hack...
{
   int i, N;
   double x,y,z,x2,y2;

   for (i = 0; i < MAXN; i++) {
     if (inp[i].z < 0.0001) break;
   }
   N = i-1;

   SplineKnots(knots,N,T);
   SplineCurve(inp,N,knots,T,outp,RESOLUTION);
return;
   printf("static int xyflightpath[RES*2] = {\n");
   for (i=0;i<RESOLUTION;i++) {

      x=outp[i].x; y=outp[i].y; z=outp[i].z;

      x2 = ((PSCL*x)/(PSCL+x)); y2 = ((PSCL*y)/(PSCL+y));
      printf("  %d, %d,\n",(int)x2, (int)y2);
   }
   printf("};\n");

return;

   /* Display the curve */
   printf("#define RES %d\n",RESOLUTION);
#ifdef VERBOSE
   /* Control point polygon */
   printf("#define POINTS %d\nstatic double controls[POINTS*3] = {\n", N+1);
   for (i=0;i<=N;i++)
      printf("  %g, %g, %g,\n",inp[i].x,inp[i].y,inp[i].z);
   printf("}\n");



   /* Display the curve */
   printf("static double curve[RES*3] = {\n",RESOLUTION);
   for (i=0;i<RESOLUTION;i++)
      printf("  %g, %g, %g,\n",outp[i].x,outp[i].y,outp[i].z);
   printf("}\n");
#endif
   printf("static int xyflightpath[RES*2] = {\n");
   for (i=0;i<RESOLUTION;i++) {

      x=outp[i].x; y=outp[i].y; z=outp[i].z;

      x2 = ((PSCL*x)/(PSCL+x)); y2 = ((PSCL*y)/(PSCL+y));
      printf("  %d, %d,\n",(int)x2, (int)y2);
   }
   printf("};\n");
}
//------------------------------------------------ END FLIGHTPATH GENERATOR