(* www.jdawiseman.com/2016/20160420_Arc_Radius_Droog.nb *)
(* www.jdawiseman.com/2016/20160420_Arc_Radius_Droog.pdf *)
(* www.jdawiseman.com/author.html *)

Quiet[Remove["Global`*"],{Remove::rmnsm}];Print["Mathematica $Version = “",$Version,"”"];Print["Execution time = ",DateString[DateList[],{"Hour",":","Minute"," on ","DayNameShort"," ","Day"," ","MonthNameShort"," ","Year"}]];

KnotsFromCoeffs[{c0_,c1_,c2_,c3_}]={c0,c0+c1/3,c0+(2c1+c2)/3,c0+c1+c2+c3};
CoeffsFromKnots[{z0_,z1_,z2_,z3_}]={z0,-3 z0+3 z1,3 z0-6 z1+3 z2,-z0+3 z1-3 z2+z3};

(* Simple checks *)
And[
	{z0,z1,z2,z3}==KnotsFromCoeffs[CoeffsFromKnots[{z0,z1,z2,z3}]],
	{c0,c1,c2,c3}==CoeffsFromKnots[KnotsFromCoeffs[{c0,c1,c2,c3}]]
]//Expand

(* https://groups.google.com/forum/#!topic/comp.lang.postscript/YiuM6lj5ngY *)

xAdobe=Table[t^i,{i,0,3}].CoeffsFromKnots[{1,1,552/1000,0}];
yAdobe=Table[t^i,{i,0,3}].CoeffsFromKnots[{0,552/1000,1,1}];
rSquaredAdobe=(xAdobe^2+yAdobe^2)//FullSimplify

Sqrt[rSquaredAdobe]/.{t->1/2}

N[Sqrt[rSquaredAdobe/.{t->1/2}],25]

Plot[Sqrt[rSquaredAdobe]-1,{t,0,1}]

{xDroog,yDroog}={
	Table[t^i,{i,0,3}].CoeffsFromKnots[{x0,x1,x2,x3}],
	Table[t^i,{i,0,3}].CoeffsFromKnots[{y0,y1,y2,y3}]
}//.{
	x0->Cos[halfA],
	y0->Sin[halfA],
	x1->(4-Cos[halfA])/3,
	y1->((3-Cos[halfA])(1-Cos[halfA]))/Sin[halfA],
	x2->x1, y2->-y1,
	x3->x0 ,y3->-y0
};
rSquaredDroog=(xDroog^2+yDroog^2)/.{halfA->π/4}//Simplify

Plot[Sqrt[rSquaredDroog]-1,{t,0,1}]

(* The above chart shows that the latest incarnation of the Droog formula is wrong. Very wrong. So let’s help. *)

{xMiddle,yMiddle}=({
	Table[t^i,{i,0,3}].CoeffsFromKnots[{x0,x1,x2,x3}],
	Table[t^i,{i,0,3}].CoeffsFromKnots[{y0,y1,y2,y3}]
}//.{
	z->Tan[a/4]×4/3,
	x0->Cos[a/2],
	y0->Sin[a/2],
	x1->x0+z×Sin[a/2],
	y1->y0-z×Cos[a/2],  (* Droog erroneously has a ‘+’ sign here *)
	x2->x1, y2->-y1,
	x3->x0 ,y3->-y0
})//Simplify;
rSquaredMiddle=(xMiddle^2+yMiddle^2)/.{a->π/2}//FullSimplify

Plot[{Sqrt[rSquaredAdobe]-1,Sqrt[rSquaredMiddle]-1},{t,0,1}]

Solve[0==D[rSquaredMiddle/.{a->π/2},t],t]

(rSquaredMiddle/.{a->π/2}/.%[[4]])//FullSimplify

N[Sqrt[%]-1,25] (* Worst of middle-match errors *)