/****************************************************************************
*
*  This module implements the code for Bezier bicubic patch shapes
*
*  This file was written origionally by Alexander Enzmann.  He wrote the code 
*  for bezier bicubic patches and generously provided povray with these 
*  enhancements.
*
*  This file was converted and edited for use in a povray include file by
*  Loial Raven.
*
*****************************************************************************/
/* bicubic_line.inc
 this is an exact conversion from the .c file, except it makes a line instead
 of a surface. 

the use is:
  bicubic_line (points, rad)
where points is an array of four points and the rad is the radius of the 
spheresweep

the bezier curve created by this touches either end point and are waited
directly towards the two middle points.

there is also a quality modifier, all you have to do is #declare it before you 
use the macro. the number is the inverse of the number of segments in the
sweep.
you can define it like this:
  #declare bicubic_line_quality = 1/20;
for a 20 segment line
*/


#declare C = array[4] { 1.0, 3.0, 3.0, 1.0 }

#macro bezier_value(points, u0)
  #local P = <0,0,0>;
	#local i = 0;
  #local c = 0.0;
	#local _t = 0.0; #local ut = 0.0;
  #local _u = array[4]  { 0.0, 0.0, 0.0, 0.0 }
	#local uu  = array[4] { 0.0, 0.0, 0.0, 0.0 }
  #local du = array[4]  { 0.0, 0.0, 0.0, 0.0 }
	#local duu = array[4] { 0.0, 0.0, 0.0, 0.0 }
  #local U1 = <0,0,0>;
  
  /* Calculate binomial coefficients times coordinate positions. */
  
  #local _u[0] = 1.0; 
	#local uu[0] = 1.0; 
	#local du[0] = 0.0; 
	#local duu[0] = 0.0;

  #local i = 1;
	#while (i < 4)
    #local _u[i] = _u[i - 1] * u0;  
		#local uu[i] = uu[i - 1] * (1.0 - u0);

    #local du[i] = i * _u[i - 1];  
		#local duu[i] = -i * uu[i - 1];
		#local i = i + 1;
	#end

	/* Now evaluate position and tangents based on control points. */
	
	#local P  = <0,0,0>;
	#local U1 = <0,0,0>;
	
	#local i = 0;
	#while (i < 4)
		#local c = C[i];
		
		#local ut = _u[i] * uu[3 - i];
		
		#local _t = c * ut;
		
		#local P = P + _t*points[i];
		
		#local _t = c * (du[i] * uu[3 - i] + _u[i] * duu[3 - i]);
		
		#local U1 = U1 + _t*points[i];
		
		#local i = i + 1;
	#end
	P
#end

#macro bicubic_line (points, rad)
	#ifndef (bicubic_line_quality) #declare bicubic_line_quality = 0.05; #end
	sphere_sweep {
		linear_sphere_sweep
		#local step = bicubic_line_quality;
		1/step+1
		#local i = 0.0;
		#while (i <= 1)
			bezier_value(points,i), rad
			#local i = i + step;
		#end
		points[3], rad
		sphere_sweep_depth_tolerance 1.0e-3
	}	
#end