Cubic case

(Subscript[B, 0])[u_]=(1-u)^3;(Subscript[B, 1])[u_]=3(1-u)^2u;(Subscript[B, 2])[u_]=3(1-u)u^2;(Subscript[B, 3])[u_]=u^3;
(Subscript[b, 0])={Cos[φ],-Sin[φ]};(Subscript[b, 1])={ξ,-η};
(Subscript[b, 2])={ξ,η};(Subscript[b, 3])={Cos[φ],Sin[φ]};
b30[t_,φ_,ξ_,η_]=Sum[(Subscript[b, j])(Subscript[B, j])[(t+1)/2],{j,0,3}];
e30[t_,φ_,ξ_,η_]=b30[t,φ,ξ,η][[1]]^2+b30[t,φ,ξ,η][[2]]^2-1//Simplify;

ψ[t_]:=(1/16)((t^2)-1)((((3 η-s))^2)(t^4)+(16(s^2)-9 (((η+s))^2)+9(((ξ-c))^2))(t^2)+(16-(((3 ξ+c))^2)));

ψ[t] -e30[t,φ,ξ,η]/.{c->Cos[φ],s->Sin[φ]}//Simplify

Error function of the best interpolant

Solve[ChebyshevT[6,t]==0,t]

χ[t_]:=2(26-15 (Sqrt[3]))ChebyshevT[6,((Sqrt[2+(Sqrt[3])])/2)t]//Simplify

((t^2)-1) ((t^2)-((((Sqrt[3])-1))^2))((t^2)-(((2-(Sqrt[3])))^2))-χ[t]//Simplify

Function f

We eliminate (η^2)

η1=η/.Solve[(((1-(v^2)))^2)(v^2)ψ[u]-(((1-(u^2)))^2)(u^2)ψ[v]==0/.{u->(Sqrt[3])-1,v->2-(Sqrt[3])},η][[1,1]]//FullSimplify

Equation (2)

(1/s) (((2+(Sqrt[3]))/8)(((3ξ+c))^2)-ξ c-3-2(Sqrt[3]))-η1/.s->(Sqrt[1-(c^2)])//Simplify

f3:=(s^2)(ψ[u]-ψ[v])/.{u->(Sqrt[3])-1,v->2-(Sqrt[3]),η->(1/s) (((2+(Sqrt[3]))/8)(((3ξ+c))^2)-ξ c-3-2(Sqrt[3]))}/.s->(Sqrt[1-(c^2)])//Simplify

f[ξ_]:=243(ξ^3)-27c(11-16(Sqrt[3]))(ξ^2)-3(32(1+2(Sqrt[3]))-3(81-32(Sqrt[3]))(c^2))ξ-32(13+2(Sqrt[3]))c-(163-112(Sqrt[3]))(c^3)

(7 (Sqrt[3])-12)(9/512) (ξ-c)f[ξ] -f3//Simplify

Position of zeros of function f for c ∈ (0,1)

Plot[{-(1/3)(4+c),-(1/9)(8(Sqrt[3])-1)c,-(1/3)c,c,(1/3)(4-c)},{c,0,1}]

N[7-4 (Sqrt[3])]

f[-(1/3)(4+c)]-(-64 (7-4 (Sqrt[3])) (((1+c))^3))//Simplify

f[-(1/9)(8(Sqrt[3])-1)c]-(256/3) c (1-(c^2))//Simplify

f[-(1/3)c]-(-(64 c (6+(7-4 (Sqrt[3])) (c^2))))//Simplify

f[c]-(-256 (2+(Sqrt[3])) c (1-(c^2)))//Simplify

f[(1/3)(4-c)]-64 (7-4 (Sqrt[3])) (((1-c))^3)//Simplify

Zeros of function f for c=0

Solve[f[ξ]==0/.c->0,ξ]

{-(1/3)(4+c),-(4/9) (Sqrt[2 (1+2 (Sqrt[3]))]),-(1/9)(8(Sqrt[3])-1)c}/.c->0//N

{c,(4/9) (Sqrt[2 (1+2 (Sqrt[3]))]),(1/3)(4-c)}/.c->0//N

Necessary condition on p to be good approximant of c (d) is x(0) > 0 (x(0) < 0)

b30[0,φ,ξ,η][[1]]-(1/4)( 3ξ+c)/.c->Cos[φ]//Simplify

Coefficient[ψ[t],t,0]-(-(1/16) (16-(((3 ξ+c))^2)))//Simplify

Parameters for the best interpolant when c=1, i.e., 2π-φ=2π

Solve[(((1-(v^2)))^2)(v^2)ψ[u]-(((1-(u^2)))^2)(u^2)ψ[v]==0/.{u->(Sqrt[3])-1,v->2-(Sqrt[3]),c->1,s->0},ξ]//Simplify

Solve[ψ[u]==0/.{u->(Sqrt[3])-1,v->2-(Sqrt[3]),c->1,s->0,ξ->-(1/9) (8 (Sqrt[3])-1)},η]//Simplify

Table of best interpolants of c with corresponding errors

angles={(π/2),(π/3),(π/4),(π/6),(π/8),(π/12)};

sez3={{"φ","ξ","η","error"}};For[i=1,i≤Length[angles],i++,φ0=angles[[i]];ξ0=Select[ξ/.NSolve[f[ξ]==0/.c->Cos[φ0],ξ,Reals,WorkingPrecision->30],#>Cos[φ0]&&#<(1/3)(4+Cos[φ0])&][[1]];η0=(1/(Sin[φ0])) (((2+(Sqrt[3]))/8)(((3ξ0+Cos[φ0]))^2)-ξ0 Cos[φ0]-3-2(Sqrt[3]));AppendTo[sez3,{φ,NumberForm[ξ0,{6,5}],NumberForm[η0,{5,5}],ScientificForm[N[-e30[0,φ,ξ0,η0]],6]}/.φ->angles[[i]]]];

Grid[sez3,Frame->All]

You can draw the best interpolant of c, the corresponding error function and the curvature for arbitrary angle φ0 ∈ (0,(π/2)]

φ0=π/3;
ξ0=Select[ξ/.NSolve[f[ξ]==0/.c->Cos[φ0],ξ,Reals,WorkingPrecision->30],#>Cos[φ0]&&#<(1/3)(4+Cos[φ0])&][[1]];
η0=(1/(Sin[φ0])) (((2+(Sqrt[3]))/8)(((3ξ0+Cos[φ0]))^2)-ξ0 Cos[φ0]-3-2(Sqrt[3]));
GraphicsRow[{Show[ParametricPlot[{Cos[φ],Sin[φ]},{φ,-φ0,φ0},PlotStyle->{Blue,Dashed},Ticks->{{0,0.5,1,1.5},{-1,-0.5,0,0.5,1}}],ParametricPlot[{b30[t,φ0,ξ0,η0][[1]],b30[t,φ0,ξ0,η0][[2]]},{t,-1,1},PlotStyle->Red],ListPlot[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])},PlotStyle->{PointSize[0.02],Black}],Graphics[{Black,Line[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])}]}],AspectRatio->Automatic,PlotRange->All],Plot[e30[t,φ0,ξ0,η0],{t,-1,1}],Plot[Evaluate[(D[b30[t,φ0,ξ0,η0][[1]],t]D[b30[t,φ0,ξ0,η0][[2]],{t,2}]-D[b30[t,φ0,ξ0,η0][[2]],t]D[b30[t,φ0,ξ0,η0][[1]],{t,2}])/(((((D[b30[t,φ0,ξ0,η0][[1]],t])^2)+((D[b30[t,φ0,ξ0,η0][[2]],t])^2)))^(3/2))],{t,-1,1}]}]/.{φ->φ0,ξ->ξ0,η->η0}

Table of best interpolants of d with corresponding errors

angles={(π/2),(π/3),(π/4),(π/6),(π/8),(π/12)};

sez3={};For[i=1,i≤Length[angles],i++,φ0=angles[[i]];ξ0=Select[ξ/.NSolve[f[ξ]==0/.c->Cos[φ0],ξ,Reals,WorkingPrecision->30],#>-(1/3)(4+Cos[φ0])&&#<-(1/9)(8(Sqrt[3])-1)Cos[φ0]&][[1]];η0=(1/(Sin[φ0])) (((2+(Sqrt[3]))/8)(((3ξ0+Cos[φ0]))^2)-ξ0 Cos[φ0]-3-2(Sqrt[3]));PrependTo[sez3,{π-φ,NumberForm[ξ0,{6,5}],NumberForm[η0,{5,5}],ScientificForm[N[-e30[0,φ,ξ0,η0]],6]}/.φ->angles[[i]]]];
PrependTo[sez3,{"φ","ξ","η","error"}];

Grid[sez3,Frame->All]

You can draw the best interpolant of d, the corresponding error function and the curvature for arbitrary angle φ0∈ [0,(π/2)]

φ0=π/3;
ξ0=Select[ξ/.NSolve[f[ξ]==0/.c->Cos[φ0],ξ,Reals,WorkingPrecision->30],#>-(1/3)(4+Cos[φ0])&&#<-(1/9)(8(Sqrt[3])-1)Cos[φ0]&][[1]];
η0=(1/(Sin[φ0])) (((2+(Sqrt[3]))/8)(((3ξ0+Cos[φ0]))^2)-ξ0 Cos[φ0]-3-2(Sqrt[3]));
GraphicsRow[{Show[ParametricPlot[{Cos[φ],Sin[φ]},{φ,-φ0,φ0-2π},PlotStyle->{Blue,Dashed},Ticks->{{0,0.5,1,1.5},{-1,-0.5,0,0.5,1}},AxesOrigin->{0,0}],ParametricPlot[{b30[t,φ0,ξ0,η0][[1]],b30[t,φ0,ξ0,η0][[2]]},{t,-1,1},PlotStyle->Red],ListPlot[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])},PlotStyle->{PointSize[0.02],Black}],Graphics[{Black,Line[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])}]}],AspectRatio->Automatic,PlotRange->All],Plot[e30[t,φ0,ξ0,η0],{t,-1,1}],Plot[Evaluate[-(D[b30[t,φ0,ξ0,η0][[1]],t]D[b30[t,φ0,ξ0,η0][[2]],{t,2}]-D[b30[t,φ0,ξ0,η0][[2]],t]D[b30[t,φ0,ξ0,η0][[1]],{t,2}])/(((((D[b30[t,φ0,ξ0,η0][[1]],t])^2)+((D[b30[t,φ0,ξ0,η0][[2]],t])^2)))^(3/2))],{t,-1,1}]}]/.{φ->φ0,ξ->ξ0,η->η0}

The best interpolant, the corresponding error function and the curvature for whole unite circle

φ0=0;
ξ0=-(1/9) (8 (Sqrt[3])-1);
η0=(4/9) (Sqrt[38+22 (Sqrt[3])]);
GraphicsRow[{Show[ParametricPlot[{Cos[φ],Sin[φ]},{φ,-φ0,φ0-2π},PlotStyle->{Blue,Dashed},Ticks->{{0,0.5,1,1.5},{-1,-0.5,0,0.5,1}},AxesOrigin->{0,0}],ParametricPlot[{b30[t,φ0,ξ0,η0][[1]],b30[t,φ0,ξ0,η0][[2]]},{t,-1,1},PlotStyle->Red],ListPlot[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])},PlotStyle->{PointSize[0.02],Black}],Graphics[{Black,Line[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])}]}],AspectRatio->Automatic,PlotRange->All],Plot[e30[t,φ0,ξ0,η0],{t,-1,1}],Plot[Evaluate[-(D[b30[t,φ0,ξ0,η0][[1]],t]D[b30[t,φ0,ξ0,η0][[2]],{t,2}]-D[b30[t,φ0,ξ0,η0][[2]],t]D[b30[t,φ0,ξ0,η0][[1]],{t,2}])/(((((D[b30[t,φ0,ξ0,η0][[1]],t])^2)+((D[b30[t,φ0,ξ0,η0][[2]],t])^2)))^(3/2))],{t,-1,1}]}]/.{φ->φ0,ξ->ξ0,η->η0}

Quartic case

(Subscript[B, 0])[u_]=(1-u)^4;(Subscript[B, 1])[u_]=4(1-u)^3u;(Subscript[B, 2])[u_]=6(1-u)^2u^2;(Subscript[B, 3])[u_]=4(1-u)u^3;(Subscript[B, 4])[u_]=u^4;
(Subscript[b, 0])={Cos[φ],-Sin[φ]};(Subscript[b, 1])={α,-β};
(Subscript[b, 2])={γ,0};(Subscript[b, 4])={Cos[φ],Sin[φ]};(Subscript[b, 3])={α,β};
b40[t_,φ_,α_,β_,γ_]=Sum[(Subscript[b, j])(Subscript[B, j])[(t+1)/2],{j,0,4}];//Simplify;
e40[t_,φ_,α_,β_,γ_]=b40[t,φ,α,β,γ][[1]]^2+b40[t,φ,α,β,γ][[2]]^2-1//Simplify;

ψ[t_]:=-1+(1/64) (((4 (1-(t^4)) α+3 (((1-(t^2)))^2) γ+(1+6 (t^2)+(t^4)) c))^2)+(1/4) (t^2) (((2 (1-(t^2)) β+(1+(t^2)) s))^2);

ψ[t]-e40[t,φ,α,β,γ]/.{c->Cos[φ],s->Sin[φ]}//Simplify

Error function of the best interpolant

Solve[ChebyshevT[8,t]==0,t]

χ[t_]:=(2/(((2+(Sqrt[2+(Sqrt[2])])))^4))ChebyshevT[8,(1/2) (Sqrt[2+(Sqrt[2+(Sqrt[2])])])t];

(Subscript[u, 1])=(Sqrt[2(2+(Sqrt[2]))])-1-(Sqrt[2]);(Subscript[u, 2])=(Sqrt[2+(Sqrt[2])])-1;(Subscript[u, 3])=1+(Sqrt[2])-(Sqrt[2+(Sqrt[2])]);

χ[t]-((t^2)-1)((t^2)-((Subscript[u, 1])^2))((t^2)-((Subscript[u, 2])^2))((t^2)-((Subscript[u, 3])^2))//Simplify

(Subscript[σ, 1])=(1-((Subscript[u, 1])^2))+(1-((Subscript[u, 2])^2))+(1-((Subscript[u, 3])^2));
(Subscript[σ, 2])=(1-((Subscript[u, 1])^2))(1-((Subscript[u, 2])^2))+(1-((Subscript[u, 1])^2))(1-((Subscript[u, 3])^2))+(1-((Subscript[u, 2])^2))(1-((Subscript[u, 3])^2));
(Subscript[σ, 3])=(1-((Subscript[u, 1])^2))(1-((Subscript[u, 2])^2))(1-((Subscript[u, 3])^2));

Grid[{{"(Subscript[σ, 1])","(Subscript[σ, 2])","(Subscript[σ, 3])"},{FullSimplify[(Subscript[σ, 1])],ToRadicals[FullSimplify[(Subscript[σ, 2])]],ToRadicals[FullSimplify[(Subscript[σ, 3])]]},{N[(Subscript[σ, 1])],N[(Subscript[σ, 2])],N[(Subscript[σ, 3])]}},Frame->All]

Function f

We eliminate the parameter β and we get the equation (4)

term4=64((((Subscript[u, 2]^2)(Subscript[u, 3]^2) ψ[(Subscript[u, 1])])/(((Subscript[u, 2]^2)-(Subscript[u, 1]^2))((Subscript[u, 1]^2)-(Subscript[u, 3]^2))(1-(Subscript[u, 1]^2))))+(((Subscript[u, 1]^2)(Subscript[u, 3]^2) ψ[(Subscript[u, 2])])/(((Subscript[u, 3]^2)-(Subscript[u, 2]^2))((Subscript[u, 2]^2)-(Subscript[u, 1]^2))(1-(Subscript[u, 2]^2))))+(((Subscript[u, 1]^2)(Subscript[u, 2]^2) ψ[(Subscript[u, 3])])/(((Subscript[u, 1]^2)-(Subscript[u, 3]^2))((Subscript[u, 3]^2)-(Subscript[u, 2]^2))(1-(Subscript[u, 3]^2)))));term4simpl=64+(Subscript[u, 1]^2)(Subscript[u, 2]^2)(Subscript[u, 3]^2)(((4α-3γ-c))^2)-(((4α+3γ+c))^2);

term4-term4simpl/.s->(Sqrt[1-(c^2)])//Simplify

Solve[4α-3γ-c==x&&4α+3γ+c==y,{α,γ}]

Equation (5)

term5=(((Subscript[u, 3]^2)(((1-(Subscript[u, 3]^2)))^2)-(Subscript[u, 2]^2)(((1-(Subscript[u, 2]^2)))^2))ψ[(Subscript[u, 1])]+((Subscript[u, 1]^2)(((1-(Subscript[u, 1]^2)))^2)-(Subscript[u, 3]^2)(((1-(Subscript[u, 3]^2)))^2))ψ[(Subscript[u, 2])]+((Subscript[u, 2]^2)(((1-(Subscript[u, 2]^2)))^2)-(Subscript[u, 1]^2)(((1-(Subscript[u, 1]^2)))^2))ψ[(Subscript[u, 3])])-((Subscript[σ, 2])/64)((Subscript[u, 1]^2)-(Subscript[u, 2]^2))((Subscript[u, 2]^2)-(Subscript[u, 3]^2))((Subscript[u, 3]^2)-(Subscript[u, 1]^2))(64+(Subscript[u, 1]^2)(Subscript[u, 2]^2)(Subscript[u, 3]^2)(x^2)-(y^2))/.{α->((x+y)/8),γ->(1/6) (-2 c-x+y)};
term5simpl=-1+(1/128)(Subscript[σ, 3])(x^2)+(1/8)(x+y)c+β s;

term5-2((Subscript[u, 1]^2)-(Subscript[u, 2]^2))((Subscript[u, 2]^2)-(Subscript[u, 3]^2))((Subscript[u, 3]^2)-(Subscript[u, 1]^2))((Subscript[σ, 2])-(Subscript[σ, 1])+1)term5simpl/.s->(Sqrt[1-(c^2)])//Simplify

Equation before (6)

term6a=(((ψ[(Subscript[u, 1])])/(1-(Subscript[u, 1]^2)))+((ψ[(Subscript[u, 2])])/(1-(Subscript[u, 2]^2)))+((ψ[(Subscript[u, 3])])/(1-(Subscript[u, 3]^2))))+(1/64) (Subscript[σ, 1])(64+(Subscript[u, 1]^2)(Subscript[u, 2]^2)(Subscript[u, 3]^2)(x^2)-(y^2))-((Subscript[σ, 1]^2)-3(Subscript[σ, 1])-2(Subscript[σ, 2])+6)(-1+(1/128)(Subscript[σ, 3])(x^2)+(1/8)c(x+y)+β s)/.{α->((x+y)/8),γ->(1/6) (-2 c-x+y)};
term6asimpl= 96-4 x y-128 (β^2)+16 c(x-y) +32 (c^2)-(4-2(Subscript[σ, 1])+(Subscript[σ, 3]))(x^2);

term6a-(1/128) ((Subscript[σ, 1]^2)-2(Subscript[σ, 2])-(Subscript[σ, 1]))term6asimpl/.s->(Sqrt[1-(c^2)])//Simplify

Equation (6)

term6=(1/128)(s^2)term6asimpl-((((1-(1/128)(Subscript[σ, 3])(x^2)-(1/8)c(x+y)))^2)-(β^2)(s^2))-((c^2)/64)term4simpl/.{α->((x+y)/8),γ->(1/6) (-2 c-x+y)};
term6simpl=(-(x/32)+(c/8)+((c^3)/8)-(1/512) (Subscript[σ, 3])c (x^2))y-(1/4)(((((1/64)(Subscript[σ, 3])(x^2)+(1/4)c x-(s^2)))^2)-(1/16)(1-(Subscript[σ, 2])+2(Subscript[σ, 3]))(c^2)(x^2)+(1/16)(2-(Subscript[σ, 1]))(x^2)- c(x-8 c));

term6-term6simpl/.s->(Sqrt[1-(c^2)])//Simplify

f[x_]:=(1/16)(((((((1/64)(Subscript[σ, 3])(x^2)+(1/4)c x-1+(c^2)))^2)-(1/16)(1-(Subscript[σ, 2])+2(Subscript[σ, 3]))(c^2)(x^2)+(1/16)(2-(Subscript[σ, 1]))(x^2)- c(x-8 c) ))^2)-(64+(Subscript[u, 1]^2)(Subscript[u, 2]^2)(Subscript[u, 3]^2)(x^2))((((x/32)-(c/8)-((c^3)/8)+(1/512)  c (Subscript[σ, 3])(x^2)))^2)

Lemma 3

f[-(((1-c))^2)] < 0

PolynomialRemainder[f[-(((1-c))^2)],(((1-c))^4),c]

Plot[Evaluate[PolynomialQuotient[f[-(((1-c))^2)],(((1-c))^4),c]],{c,0,1}]

PolynomialQuotient[f[-(((1-c))^2)],(((1-c))^4),c]/.c->0//Simplify//N

f[0] > 0

f[0]//Simplify

f[4c(1+(c^2))] > 0

Plot[f[4c(1+(c^2))],{c,0,1},PlotRange->{0,1}]

f[(6-3(Sqrt[2])+(Sqrt[68-46(Sqrt[2])]))(((1+(c^4)))^2)] < 0

Plot[f[(6-3(Sqrt[2])+(Sqrt[68-46(Sqrt[2])]))(((1+(c^4)))^2)],{c,0,1}]

f[x]/.c->1/.x->4(6-3(Sqrt[2])+(Sqrt[68-46(Sqrt[2])]))//FullSimplify

Function f is increasing on (-(((1-c))^2),0)

g[t_]:=f[-t(((1-c))^2)];(Subscript[g, 1])[t_]:=(1/16)(((1+c))^2)(5+(c^2))c+(1/128)(12-56(c^2)+28(c^4)+2(Subscript[σ, 1])(1-(c^2))(1+5(c^2)+2(c^4))+2 (Subscript[σ, 2])(((1-(c^2)))^2)(c^2)+(Subscript[σ, 3])(((1-(c^2)))^3))t+(3/1024)(4(Subscript[σ, 1])(1-3(c^2))+4 (Subscript[σ, 2])(2-(c^2)+(c^4))+(1-(c^2))(-16+(-9+5(c^2))(Subscript[σ, 3])))(((1-c))^2)c (t^2);
(Subscript[g, 2])[t_]:=((((5/2)(Subscript[σ, 3])- (Subscript[σ, 2]) c))^2)(t^3);
(Subscript[g, 3])[t_]:=((Subscript[σ, 2])(32(1-(c^2))+(Subscript[σ, 3])(5-8c+8(c^3))c)-(1/4)(Subscript[σ, 3])(64(1-5(c^2)+2(c^4))+(Subscript[σ, 3])(37-24(c^2)+12(c^4)))+8(Subscript[σ, 1])(2 (Subscript[σ, 2]) (c^2)-(Subscript[σ, 3])(1+(c^2)-2(c^4)))-8(Subscript[σ, 1]^2)-((c^2)+8(c^4)) (Subscript[σ, 2]^2))(t^3)-(5/8)(Subscript[σ, 3])(8 (Subscript[σ, 2])+5(c^2)(Subscript[σ, 3]))(((1-c))^2)c (t^4)-(3/64)(Subscript[σ, 3]^2)(4+4(c^2)(Subscript[σ, 1])+(c^2)(Subscript[σ, 3]))(((1-c))^4)(t^5)-(1/4096)(Subscript[σ, 3]^4)(((1-c))^8)(t^7);
(Subscript[g, 4])[t_]:=(5/8)(Subscript[σ, 3])(4(Subscript[σ, 1])+4(c^2) (Subscript[σ, 2])+3(Subscript[σ, 3]))(((1-c))^2)c (t^4)+(3/64)(Subscript[σ, 3]^2)(2(Subscript[σ, 1])+2(c^2) (Subscript[σ, 2])+(Subscript[σ, 3]))(((1-c))^4)(t^5)+(7/512)(Subscript[σ, 3]^3)(((1-c))^6)c (t^6);

-((((1-c))^4)(Subscript[g, 1])[t]+(1/8192)(((1-c))^8)((Subscript[g, 2])[t]+(Subscript[g, 3])[t]+(Subscript[g, 4])[t]))-D[g[t],t]/.(Subscript[u, 1]^2) (Subscript[u, 2]^2) (Subscript[u, 3]^2)->(1 - (Subscript[σ, 1])+ (Subscript[σ, 2])- (Subscript[σ, 3]))/.s->(Sqrt[1-(c^2)])//Simplify

(Subscript[g, 2])[t]>0

(Subscript[g, 2])[t]//N

(Subscript[g, 4])[t]≥0

(Subscript[g, 4])[t]//N

(Subscript[g, 1])[t]≥0

CoefficientList[(Subscript[g, 1])[t],t][[1]]

Plot[Total[CoefficientList[(Subscript[g, 1])[t],t][[1;;2]]]//Simplify,{c,0,1}]

Plot[Total[CoefficientList[(Subscript[g, 1])[t],t][[1;;3]]]//Simplify,{c,0,1}]

(Subscript[g, 3])[t]≥0

Coefficient[(Subscript[g, 3])[t],(t^4)]//N
Coefficient[(Subscript[g, 3])[t],(t^5)]//N
Coefficient[(Subscript[g, 3])[t],(t^7)]//N

Plot[(Subscript[g, 3])[1],{c,0,1}]

Function f is increasing on (0,4(c^2))

(Subscript[g, 1])[x_]=-(1/(2^24))(21(Subscript[σ, 3]^4)(x^4)+84c (Subscript[σ, 3]^3)(10+c (Subscript[σ, 3]))(x^3)+48(Subscript[σ, 3]^2)(160+80(Subscript[σ, 1])-80(Subscript[σ, 2])+70(Subscript[σ, 3])+7(Subscript[σ, 3]^2))(x^2))(4(c^2)-x);
(Subscript[g, 2])[x_]=-((3(Subscript[σ, 3]^2))/(2^20))(80(1+c)(2(Subscript[σ, 1])-(Subscript[σ, 2]))+(Subscript[σ, 3])(10(11+11c+7(c^2))+7(1+c+(c^2)+(c^3))(Subscript[σ, 3])))(1-c)(x^3);
(Subscript[g, 3])[x_]=-((3(Subscript[σ, 3])(1-c))/(2^18))(80(Subscript[σ, 2])(4((c^2)+c-1)-(1+c)(Subscript[σ, 3]))+80(Subscript[σ, 1])(4+(1+c)(Subscript[σ, 3]))+(Subscript[σ, 3])(70(1+c)(Subscript[σ, 3])-80c(3+5c)+7(1+c)(Subscript[σ, 3]^2)))(x^2);
(Subscript[g, 0])[x_]=D[f[x],{x,3}]-(Subscript[g, 1])[x]-(Subscript[g, 2])[x]-(Subscript[g, 3])[x]//Simplify;

(Subscript[g, 1])[x] < 0

CoefficientList[21(Subscript[σ, 3]^4)(x^4)+84c (Subscript[σ, 3]^3)(10+c (Subscript[σ, 3]))(x^3)+48(Subscript[σ, 3]^2)(160+80(Subscript[σ, 1])-80(Subscript[σ, 2])+70(Subscript[σ, 3])+7(Subscript[σ, 3]^2))(x^2),x]//N

(Subscript[g, 2])[x] ≤ 0

CoefficientList[80(1+c)(2(Subscript[σ, 1])-(Subscript[σ, 2]))+(Subscript[σ, 3])(10(11+11c+7(c^2))+7(1+c+(c^2)+(c^3))(Subscript[σ, 3])),c]//N

(Subscript[g, 3])[x] ≤ 0

CoefficientList[80(Subscript[σ, 2])(4((c^2)+c-1)-(1+c)(Subscript[σ, 3]))+80(Subscript[σ, 1])(4+(1+c)(Subscript[σ, 3]))+(Subscript[σ, 3])(70(1+c)(Subscript[σ, 3])-80c(3+5c)+7(1+c)(Subscript[σ, 3]^2)),c]//N

(Subscript[g, 0])[x] ≤ 0

Length[ CoefficientList[(Subscript[g, 0])[x],x]]

CoefficientList[(Subscript[g, 0])[x],x][[2]]//N

(Subscript[g, 0])[0]//Simplify//N

(Subscript[g, 0])[4(c^2)]//Simplify//N

Plot[Evaluate[PolynomialQuotient[(Subscript[g, 0])[4(c^2)],c,c]],{c,0,1}]

Function f is concave on (4(c^2),8c)

(Subscript[g, 1])[x_]=-(1/(2^25))(7(Subscript[σ, 3]^4)(x^5)+56c (Subscript[σ, 3]^3)(6+(Subscript[σ, 3]))(x^4)+64(Subscript[σ, 3]^2)(60+30(Subscript[σ, 1])-30(Subscript[σ, 2])+42(Subscript[σ, 3])+7(Subscript[σ, 3]^2))(x^3)+512c (Subscript[σ, 3]) (7((((Subscript[σ, 3])+2))^3)+80 ((Subscript[σ, 1])-2 (Subscript[σ, 2]))+30((Subscript[σ, 1])-(Subscript[σ, 2]) )(Subscript[σ, 3])-56+36 (Subscript[σ, 3])+20(4 (Subscript[σ, 2])-5 (Subscript[σ, 3]))(c^2))(x^2))(8c-x);
(Subscript[g, 2])[x_]=-(1/(2^19))(Subscript[σ, 3]^2)(60 (Subscript[σ, 1])-30 (Subscript[σ, 2])+(Subscript[σ, 3]) (57+7 (Subscript[σ, 3])))(1-(c^2))(x^4);
(Subscript[g, 3])[x_]=-(1/(2^13))(7(Subscript[σ, 3])((((Subscript[σ, 3])+2))^3)+(30 (Subscript[σ, 1])-30 (Subscript[σ, 2])-73)(Subscript[σ, 3]^2)+24(4-2 (Subscript[σ, 1])+(Subscript[σ, 2]))(Subscript[σ, 2])+8(7 (Subscript[σ, 1])-10 (Subscript[σ, 2])-25)(Subscript[σ, 3]) + (24 (Subscript[σ, 2]^2)+96 (Subscript[σ, 3])-48 (Subscript[σ, 1]) (Subscript[σ, 3])+56 (Subscript[σ, 2]) (Subscript[σ, 3])-91 (Subscript[σ, 3]^2))(c^2))(1-(c^2))(x^2);
(Subscript[g, 0])[x_]=D[f[x],{x,2}]-(Subscript[g, 1])[x]-(Subscript[g, 2])[x]-(Subscript[g, 3])[x]//Simplify;

(Subscript[g, 1])[x] < 0

(Subscript[g, 1])[x]//N

(Subscript[g, 2])[x] < 0

(Subscript[g, 2])[x]//N

(Subscript[g, 3])[x] ≤ 0

(Subscript[g, 3])[x]//N

(Subscript[g, 0])[x] ≤ 0

Length[ CoefficientList[(Subscript[g, 0])[x],x]]

CoefficientList[(Subscript[g, 0])[x],x][[2]]//N

Plot[(Subscript[g, 0])[4(c^2)],{c,0,1}]

Plot[(Subscript[g, 0])[8c],{c,0,1}]

Function f is decreasing on (8c,4 (6-3 (Sqrt[2])+(Sqrt[68-46 (Sqrt[2])])))

g[t_]=f'[4 (6-3 (Sqrt[2])+(Sqrt[68-46 (Sqrt[2])])-2c)t+8c];

Length[CoefficientList[g[t],t]]

g[t]=(Subscript[a, 0])+(Subscript[a, 1])t+(Subscript[a, 2])(t^2)+(Subscript[a, 3])(t^3)+(Subscript[a, 4])(t^4)+(Subscript[a, 5])(t^5)+(Subscript[a, 6])(t^6)+(Subscript[a, 7])(t^7)<0 if (Subscript[a, 0])+...+(Subscript[a, k])<0 for all k=0,...,7

(Subscript[a, 0])<0

g[0]/.c->0//Simplify

Plot[Total[CoefficientList[g[t],t][[1;;1]]],{c,0,1},PlotRange->{-0.3,0.1}]

(Subscript[a, 0])+(Subscript[a, 1])<0

Plot[Total[CoefficientList[g[t],t][[1;;2]]],{c,0,1},PlotRange->{-1.9,0}]

(Subscript[a, 0])+(Subscript[a, 1])+(Subscript[a, 2])<0

Plot[Total[CoefficientList[g[t],t][[1;;3]]],{c,0,1},PlotRange->{-2.1,0}]

(Subscript[a, 0])+(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3])<0

Plot[Total[CoefficientList[g[t],t][[1;;4]]],{c,0,1},PlotRange->{-2.2,0}]

(Subscript[a, 0])+(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3])+(Subscript[a, 4])<0

Plot[Total[CoefficientList[g[t],t][[1;;5]]],{c,0,1},PlotRange->{-2,0}]

Total[CoefficientList[g[t],t][[1;;5]]]/.c->1//N

(Subscript[a, 0])+(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3])+(Subscript[a, 4])+(Subscript[a, 5])<0

Plot[Total[CoefficientList[g[t],t][[1;;6]]],{c,0,1},PlotRange->{-2,0}]

Total[CoefficientList[g[t],t][[1;;6]]]/.c->1//N

(Subscript[a, 0])+(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3])+(Subscript[a, 4])+(Subscript[a, 5])+(Subscript[a, 6])<0

Plot[Total[CoefficientList[g[t],t][[1;;7]]],{c,0,1},PlotRange->{-2,0}]

Total[CoefficientList[g[t],t][[1;;7]]]/.c->1//N

(Subscript[a, 0])+(Subscript[a, 1])+(Subscript[a, 2])+(Subscript[a, 3])+(Subscript[a, 4])+(Subscript[a, 5])+(Subscript[a, 6])+(Subscript[a, 7])<0

Plot[Total[CoefficientList[g[t],t][[1;;8]]],{c,0,1},PlotRange->{-2,0}]

g[t]/.c->1/.t->1//FullSimplify

Lemma 4

x  ∈ [-14,0)

(Subscript[q, 1])[x_]=-(x/32)+(c/8)+((c^3)/8)-(1/512) (Subscript[σ, 3])c (x^2);
(Subscript[q, 2])[x_]=(1/4)(((((1/64)(Subscript[σ, 3])(x^2)+(1/4)c x-1+(c^2)))^2));
(Subscript[q, 3])[x_]=(1/4)(-(1/16)(1-(Subscript[σ, 2])+2(Subscript[σ, 3]))(c^2)(x^2)+(1/16)(2-(Subscript[σ, 1]))(x^2)- c(x-8 c));

(Subscript[q, 1])>0

D[(Subscript[q, 1])[x],x]//N

D[(Subscript[q, 1])[x],x]/.x->-14//N

(Subscript[q, 1])[0]//Simplify

(Subscript[q, 3])>0

D[(Subscript[q, 3])[x],x]//Simplify//N

(Subscript[q, 3])[0]//Simplify

x  ∈ (4c(1+(c^2)),(6-3 (Sqrt[2])+(Sqrt[68-46 (Sqrt[2])]))(((1+(c^4)))^2)]

(Subscript[q, 1])<0

D[(Subscript[q, 1])[x],x]//N

(Subscript[q, 1])[8c]//Simplify//N

(Subscript[q, o])<0

(Subscript[q, 0])[x_]=(Subscript[q, 2])[x]+(Subscript[q, 3])[x]-(1/(2^14))(32c (Subscript[σ, 3])(x^3)+(Subscript[σ, 3]^2)(x^4))//Simplify;

Length[CoefficientList[(Subscript[q, 0])[x],x]]

Plot[Coefficient[(Subscript[q, 0])[x],x,2],{c,0.9,1}]

Plot[Discriminant[(Subscript[q, 0])[x],x],{c,0.9,1}]

CoefficientList[D[(Subscript[q, 0])[x],x],c]//Simplify//N

Plot[(Subscript[q, 0])[(6-3 (Sqrt[2])+(Sqrt[68-46 (Sqrt[2])]))(((1+(c^4)))^2)],{c,0,0.9},PlotRange->{0,0.25}]

Table of best interpolants of arc c with corresponding errors

angles={(π/2),(π/3),(π/4),(π/6),(π/8),(π/12)};

sez4={{"φ","α","β","γ","error"}};For[i=1,i≤Length[angles],i++,φ0=angles[[i]];
x0=Select[x/.NSolve[f[x]==0/.{c->Cos[φ0],s->Sin[φ0]},x,Reals,WorkingPrecision->30],#>-(((1-Cos[φ0]))^2)&&#<0&][[1]];
y0=y/.Solve[term6simpl==0/.{x->x0,c->Cos[φ0],s->Sin[φ0]},y][[1]];
β0=β/.Solve[term5simpl==0/.{x->x0,y->y0,c->Cos[φ0],s->Sin[φ0]},β][[1]];AppendTo[sez4,{φ0,NumberForm[(1/8)(x0+y0),{6,5}],NumberForm[β0,{6,5}],NumberForm[(1/6)(y0-x0-2Cos[φ0]),{6,5}],ScientificForm[N[e40[0,φ0,(1/8)(x0+y0),β0,(1/6)(y0-x0-2Cos[φ0])]],6]}/.φ->φ0]];

Grid[sez4,Frame->All]

You can draw the best interpolant of c, the corresponding error function and the curvature for arbitrary angle φ0 ∈ (0,(π/2)]

φ0=(π/3);
x0=Select[x/.NSolve[f[x]==0/.{c->Cos[φ0],s->Sin[φ0]},x,Reals,WorkingPrecision->30],#>-(((1-Cos[φ0]))^2)&&#<0&][[1]];
y0=y/.Solve[term6simpl==0/.{x->x0,c->Cos[φ0],s->Sin[φ0]},y][[1]];
β0=β/.Solve[term5simpl==0/.{x->x0,y->y0,c->Cos[φ0],s->Sin[φ0]},β][[1]];
α0=(1/8)(x0+y0);
γ0=(1/6)(y0-x0-2Cos[φ0]);
GraphicsRow[{Show[ParametricPlot[{Cos[φ],Sin[φ]},{φ,-φ0,φ0},PlotStyle->{Blue,Dashed},Ticks->{{0,0.5,1,1.5},{-1,-0.5,0,0.5,1}}],ParametricPlot[{b40[t,φ0,α0,β0,γ0][[1]],b40[t,φ0,α0,β0,γ0][[2]]},{t,-1,1},PlotStyle->Red],ListPlot[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4])},PlotStyle->{PointSize[0.02],Black}],Graphics[{Black,Line[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4])}]}],AspectRatio->Automatic,PlotRange->All],Plot[e40[t,φ0,α0,β0,γ0],{t,-1,1}],Plot[Evaluate[(D[b40[t,φ0,α0,β0,γ0][[1]],t]D[b40[t,φ0,α0,β0,γ0][[2]],{t,2}]-D[b40[t,φ0,α0,β0,γ0][[2]],t]D[b40[t,φ0,α0,β0,γ0][[1]],{t,2}])/(((((D[b40[t,φ0,α0,β0,γ0][[1]],t])^2)+((D[b40[t,φ0,α0,β0,γ0][[2]],t])^2)))^(3/2))],{t,-1,1}]}]/.{φ->φ0,α->α0,β->β0,γ->γ0}

If φ0 < 0.6772 there are two admissible solutions with alternating simplified radial error function

φ0=(π/5);
x0={x/.FindRoot[f[x]/.{c->Cos[φ0],s->Sin[φ0]},{x,-(1/2)(((1-Cos[φ0]))^2)}],x/.FindRoot[f[x]/.{c->Cos[φ0],s->Sin[φ0]},{x,-1}]};
ima={{},{}};
For[i=1,i≤2,i++,
y0=y/.Solve[term6simpl==0/.{x->x0[[i]],c->Cos[φ0],s->Sin[φ0]},y][[1]];
β0=β/.Solve[term5simpl==0/.{x->x0[[i]],y->y0,c->Cos[φ0],s->Sin[φ0]},β][[1]];
α0=(1/8)(x0[[i]]+y0);
γ0=(1/6)(y0-x0[[i]]-2Cos[φ0]);
ima[[i]]={Show[ParametricPlot[{Cos[φ],Sin[φ]},{φ,-φ0,φ0},PlotStyle->{Blue,Dashed},Ticks->{{0,0.5,1,1.5},{-1,-0.5,0,0.5,1}}],ParametricPlot[{b40[t,φ0,α0,β0,γ0][[1]],b40[t,φ0,α0,β0,γ0][[2]]},{t,-1,1},PlotStyle->Red],ListPlot[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4])},PlotStyle->{PointSize[0.02],Black}],Graphics[{Black,Line[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4])}]}],AspectRatio->Automatic,PlotRange->All],Plot[e40[t,φ0,α0,β0,γ0],{t,-1,1}],Plot[Evaluate[(D[b40[t,φ0,α0,β0,γ0][[1]],t]D[b40[t,φ0,α0,β0,γ0][[2]],{t,2}]-D[b40[t,φ0,α0,β0,γ0][[2]],t]D[b40[t,φ0,α0,β0,γ0][[1]],{t,2}])/(((((D[b40[t,φ0,α0,β0,γ0][[1]],t])^2)+((D[b40[t,φ0,α0,β0,γ0][[2]],t])^2)))^(3/2))],{t,-1,1}]}/.{φ->φ0,α->α0,β->β0,γ->γ0};]
GraphicsGrid[ima]

Table of best interpolants of arc d with corresponding errors

angles=Reverse[{(π/2),(π/3),(π/4),(π/6),(π/8),(π/12)}];

sez4={{"φ","α","β","γ","error"}};For[i=1,i≤Length[angles],i++,φ0=angles[[i]];
x0=Select[x/.NSolve[f[x]==0/.{c->Cos[φ0],s->Sin[φ0]},x,Reals,WorkingPrecision->30],#>0&&#<14&][[1]];
y0=y/.Solve[term6simpl==0/.{x->x0,c->Cos[φ0],s->Sin[φ0]},y][[1]];
β0=β/.Solve[term5simpl==0/.{x->x0,y->y0,c->Cos[φ0],s->Sin[φ0]},β][[1]];AppendTo[sez4,{2π-2φ0,NumberForm[(1/8)(x0+y0),{6,5}],NumberForm[β0,{6,5}],NumberForm[(1/6)(y0-x0-2Cos[φ0]),{6,5}],ScientificForm[N[e40[0,φ0,(1/8)(x0+y0),β0,(1/6)(y0-x0-2Cos[φ0])]],6]}/.φ->φ0]];

Grid[sez4,Frame->All]

You can draw the best interpolant of d,  the corresponding error function and the curvature for arbitrary angle φ0 ∈ [0,(π/2)]

φ0=(π/3);
x0=Select[x/.NSolve[f[x]==0/.{c->Cos[φ0],s->Sin[φ0]},x,Reals,WorkingPrecision->30],#>0&&#<14&][[1]];
y0=y/.Solve[term6simpl==0/.{x->x0,c->Cos[φ0],s->Sin[φ0]},y][[1]];
β0=β/.Solve[term5simpl==0/.{x->x0,y->y0,c->Cos[φ0],s->Sin[φ0]},β][[1]];
α0=(1/8)(x0+y0);
γ0=(1/6)(y0-x0-2Cos[φ0]);
GraphicsRow[{Show[ParametricPlot[{Cos[φ],Sin[φ]},{φ,-φ0,φ0-2π},PlotStyle->{Blue,Dashed},Ticks->{{0,0.5,1,1.5},{-1,-0.5,0,0.5,1}}],ParametricPlot[{b40[t,φ0,α0,β0,γ0][[1]],b40[t,φ0,α0,β0,γ0][[2]]},{t,-1,1},PlotStyle->Red],ListPlot[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4])},PlotStyle->{PointSize[0.02],Black}],Graphics[{Black,Line[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4])}]}],AspectRatio->Automatic,PlotRange->All],Plot[e40[t,φ0,α0,β0,γ0],{t,-1,1}],Plot[Evaluate[-(D[b40[t,φ0,α0,β0,γ0][[1]],t]D[b40[t,φ0,α0,β0,γ0][[2]],{t,2}]-D[b40[t,φ0,α0,β0,γ0][[2]],t]D[b40[t,φ0,α0,β0,γ0][[1]],{t,2}])/(((((D[b40[t,φ0,α0,β0,γ0][[1]],t])^2)+((D[b40[t,φ0,α0,β0,γ0][[2]],t])^2)))^(3/2))],{t,-1,1}]}]/.{φ->φ0,α->α0,β->β0,γ->γ0}

The best interpolant, the corresponding error function and the curvature for whole unite circle

φ0=0;
x0=4 (6-3 (Sqrt[2])+(Sqrt[68-46 (Sqrt[2])]));
y0=y/.NSolve[term6simpl==0/.{x->x0,c->Cos[φ0],s->Sin[φ0]},y][[1]];
β0=β/.NSolve[ψ[(Subscript[u, 1])]==0/.{α->((x0+y0)/8),γ->(1/6) (-2 Cos[φ0]-x0+y0),x->x0,y->y0,c->Cos[φ0],s->Sin[φ0]},β][[1]];
α0=(1/8)(x0+y0);
γ0=(1/6)(y0-x0-2Cos[φ0]);GraphicsRow[{Show[ParametricPlot[{Cos[φ],Sin[φ]},{φ,-φ0,φ0-2π},PlotStyle->{Blue,Dashed},Ticks->{{0,0.5,1,1.5},{-1,-0.5,0,0.5,1}}],ParametricPlot[{b40[t,φ0,α0,β0,γ0][[1]],b40[t,φ0,α0,β0,γ0][[2]]},{t,-1,1},PlotStyle->Red],ListPlot[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4])},PlotStyle->{PointSize[0.02],Black}],Graphics[{Black,Line[{(Subscript[b, 0]),(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4])}]}],AspectRatio->Automatic,PlotRange->All],Plot[e40[t,φ0,α0,β0,γ0],{t,-1,1}],Plot[Evaluate[(D[b40[t,φ0,α0,β0,γ0][[1]],t]D[b40[t,φ0,α0,β0,γ0][[2]],{t,2}]-D[b40[t,φ0,α0,β0,γ0][[2]],t]D[b40[t,φ0,α0,β0,γ0][[1]],{t,2}])/(((((D[b40[t,φ0,α0,β0,γ0][[1]],t])^2)+((D[b40[t,φ0,α0,β0,γ0][[2]],t])^2)))^(3/2))],{t,-1,1}]}]/.{φ->φ0,α->α0,β->β0,γ->γ0}