Compute quadratic bezier approximation to a circle

bezier-circle.sage

var('a t')

# number of curves we want to approximate the circle with
for n_curves in range(3,11):

	# startpoint of quadratic bezier curve
	x0=cos(0)
	y0=sin(0)

	# control point of quadratic bezier curve
	x1=cos(pi/n_curves)
	y1=sin(pi/n_curves)

	# endpoint of quadratic bezier curve
	x2=cos(2*pi/n_curves)
	y2=sin(2*pi/n_curves)

	def f(a):
		#print("a:" + str(a))
		# quadratic bezier curve
		px(t)=(1-t)^2*x0+2*(1-t)*t*x1*a+t^2*x2
		py(t)=(1-t)^2*y0+2*(1-t)*t*y1*a+t^2*y2

		#val=numerical_integral(abs(px(t)^2+py(t)^2-1),0,1)
		val=numerical_integral(abs(sqrt(px(t)^2+py(t)^2)-1),0,1)
		#val=find_local_minimum(-abs((px(t)^2+py(t)^2)-1),0,1)
		#return -val[0]
		#print(val)
		return val[0]

	# fit derivative
	px(t,a)=(1-t)^2*x0+2*(1-t)*t*x1*a+t^2*x2
	b=diff(px,t)(t=0).roots(a)[0][0]

	# find bezier curve optimizing globally
	a=find_local_minimum(f,1,3)[1]

	print("vec2("+b.n(23).str()+","+str(a)+"),")

	#first quadratic bezier curve
	px(t)=(1-t)^2*x0+2*(1-t)*t*x1*a+t^2*x2
	py(t)=(1-t)^2*y0+2*(1-t)*t*y1*a+t^2*y2

	#p=parametric_plot((cos(t),sin(t)),(t,0,2*pi))
	#p+=line([(x0,y0), (a*x1,a*y1)])
	#p+=line([(a*x1,a*y1), (x2,y2)])
	#p+=parametric_plot((px(t),py(t)),(t,0,1))

	#p.save('/tmp/bezier_circle.svg')
	#os.system('xdg-open /tmp/bezier_circle.svg &')

	p=bezier_path([[(x0,y0),(a*x1,a*y1),(x2,y2)]],aspect_ratio=1)
	#p=parametric_plot((px(t),py(t)),(t,0,1))
	#p+=parametric_plot((cos(t),sin(t)),(t,0,2*pi))

	for i in range(n_curves-1):
		x0=x2
		y0=y2

		x1=cos((3+2*i)*pi/n_curves)
		y1=sin((3+2*i)*pi/n_curves)

		x2=cos((4+2*i)*pi/n_curves)
		y2=sin((4+2*i)*pi/n_curves)

		px(t)=(1-t)^2*x0+2*(1-t)*t*x1*a+t^2*x2
		py(t)=(1-t)^2*y0+2*(1-t)*t*y1*a+t^2*y2
		#p+=parametric_plot((px(t),py(t)),(t,0,1))
		#p+=bezier_path([[(x0,y0),(a*x1,a*y1),(x2,y2)]])

	#p.save('/tmp/bezier_circle2.svg')
	#os.system('xdg-open /tmp/bezier_circle2.svg &')