Geometric approximation of a sphere by parametric polynomial splines

3. Approximation of equilateral spherical triangles

(Subscript[v, 0])={c,0,(Sqrt[1-(c^2)])};
(Subscript[v, 1])={-(1/2)c,((Sqrt[3])/2)c,(Sqrt[1-(c^2)])};
(Subscript[v, 2])={-(1/2)c,-((Sqrt[3])/2)c,(Sqrt[1-(c^2)])};
For[n=2,n≤4,n++,For[i=0,i≤n,i++,For[j=0,j≤n-i,j++,(Subscript[B, i,j,n-i-j])=((n!)/(i! j!(n-i-j)!))(u^i)(v^j)(w^(n-i-j))]]]
cT=((2(Sqrt[2]))/3);
cO=((Sqrt[6])/3);
cI=(Sqrt[((2(5-(Sqrt[5])))/15)]);
R={{Cos[((2π)/3)],-Sin[((2π)/3)],0},{Sin[((2π)/3)],Cos[((2π)/3)],0},{0,0,1}};
RTet={{(5/6),-(1/(2 (Sqrt[3]))),((Sqrt[2])/3)},{(1/(2 (Sqrt[3]))),-(1/2),-(Sqrt[(2/3)])},{((Sqrt[2])/3),(Sqrt[(2/3)]),-(1/3)}};
ROct={{(2/3),-(1/(Sqrt[3])),((Sqrt[2])/3)},{(1/(Sqrt[3])),0,-(Sqrt[(2/3)])},{((Sqrt[2])/3),(Sqrt[(2/3)]),(1/3)}};
RIco={{(1/12) (9-(Sqrt[5])),-((3+(Sqrt[5]))/(4 (Sqrt[3]))),(1/3)},{((3+(Sqrt[5]))/(4 (Sqrt[3]))),(1/4) (-1+(Sqrt[5])),-(1/(Sqrt[3]))},{(1/3),(1/(Sqrt[3])),((Sqrt[5])/3)}};
Z={{-1,0,0},{0,-1,0},{0,0,-1}};
SymmetriesT={IdentityMatrix[3],RTet,RTet.RTet,R.RTet.RTet};
SymmetriesO={IdentityMatrix[3],ROct,ROct.ROct,ROct.ROct.ROct,Z,Z.ROct,Z.ROct.ROct,Z.ROct.ROct.ROct};
SymmetriesI={IdentityMatrix[3],RIco,R.RIco,R.R.RIco,RIco.RIco,R.RIco.RIco,R.R.RIco.RIco,RIco.RIco.RIco,R.RIco.RIco.RIco,R.R.RIco.RIco.RIco,Z,Z.RIco,Z.R.RIco,Z.R.R.RIco,Z.RIco.RIco,Z.R.RIco.RIco,Z.R.R.RIco.RIco,Z.RIco.RIco.RIco,Z.R.RIco.RIco.RIco,Z.R.R.RIco.RIco.RIco};
Off[LessEqual::meprec];
GaussianCurvature[f_,{u_,v_}]:=Simplify[(Det[{D[f,{u,2}],D[f,u],D[f,v]}] Det[{D[f,{v,2}],D[f,u],D[f,v]}]-Det[{D[f,u,v],D[f,u],D[f,v]}]^2)/(D[f,u].D[f,u] D[f,v].D[f,v]-(D[f,u].D[f,v])^2)^2];
GaussianCurvaturePoint[f_,{u_,v_},{a_,b_}]:=Simplify[(Det[{D[f,{u,2}],D[f,u],D[f,v]}] Det[{D[f,{v,2}],D[f,u],D[f,v]}]-Det[{D[f,u,v],D[f,u],D[f,v]}]^2)/(D[f,u].D[f,u] D[f,v].D[f,v]-(D[f,u].D[f,v])^2)^2/.{u->a,v->b}];

GraphicsRow[{Graphics[{Dashed,Line[{{0,0},{(1/2),(1/2)}}],Line[{{0,(1/2)},{1,0}}],Line[{{(1/2),0},{0,1}}],FaceForm[None],EdgeForm[Black],Triangle[{{0,0},{1,0},{0,1}}],FaceForm[Gray]}],
Show[{ParametricPlot[{u v+(u+v)(1-u-v),u v(1-u-v)},{u,0,1},{v,0,1-u},PlotStyle->Gray,PlotRange->{{0,0.4},{0,0.2}},Frame->False,Axes->False],
ParametricPlot[{u v+(u+v)(1-u-v),u v(1-u-v)}/.v->0,{u,0,1},PlotStyle->Black],ParametricPlot[{u v+(u+v)(1-u-v),u v(1-u-v)}/.v->u,{u,0,1},PlotStyle->Black]}]}]

Jac={{D[u v+ (u+v)(1-u-v),u],D[u v+ (u+v)(1-u-v),v]},{D[u v (1-u-v),u],D[u v (1-u-v),v]}}//Simplify;
Jac//MatrixForm
Solve[u-3 (u^2)+2 (u^3)-v+3 (u^2) v+3 (v^2)-3 u (v^2)-2 (v^3)==0,{u,v}]

5. Quadratic (Superscript[G, 0]) approximation

(Subscript[b, 2,0,0])=(Subscript[v, 0]); (Subscript[b, 0,2,0])=(Subscript[v, 1]); (Subscript[b, 0,0,2])=(Subscript[v, 2]);
(Subscript[b, 1,1,0])=(α/2)((Subscript[v, 0])+(Subscript[v, 1])); (Subscript[b, 0,1,1])=(α/2)((Subscript[v, 1])+(Subscript[v, 2])); (Subscript[b, 1,0,1])=(α/2)((Subscript[v, 0])+(Subscript[v, 2]));
p2[u_,v_,w_,α_]=Sum[(Subscript[B, i,j,2-i-j])(Subscript[b, i,j,2-i-j]),{i,0,2},{j,0,2-i}]//Simplify;
f[u_,v_,w_,α_]=p2[u,v,w,α].p2[u,v,w,α]-1;
g[u_,v_,w_,α_]=(Sqrt[1+f[u,v,w,α]])-1;

For[i=0,i≤2,i++,(Subscript[f, i])=Coefficient[SymmetricReduction[f[u,v,w,α],{u,v,w},{(Subscript[e, 1]),(Subscript[e, 2]),(Subscript[e, 3])}][[1]],α,i]]

{(Subscript[f, 2]),(Subscript[f, 1]),(Subscript[f, 0])}.{(α^2),α,1}

Collect[SymmetricReduction[f[u,v,w,α],{u,v,w},{(Subscript[e, 1]),(Subscript[e, 2]),(Subscript[e, 3])}][[1]],α]

({{}, {}})

(Subscript[f, 0])==(6 (c^2) (Subscript[e, 1]) (Subscript[e, 3])+(Subscript[e, 1]^4)-4 (Subscript[e, 1]^2) (Subscript[e, 2])+(4 -3 (c^2))(Subscript[e, 2]^2)-1)//Simplify
(Subscript[f, 1])==((4-3 (c^2))((Subscript[e, 1]^2) (Subscript[e, 2])- 2(Subscript[e, 2]^2) -3 (Subscript[e, 1])(Subscript[e, 3]))+12(1- (c^2))(Subscript[e, 1]) (Subscript[e, 3]))//Simplify
(Subscript[f, 2])==(12(1-(c^2))(Subscript[e, 1]) (Subscript[e, 3])+(4-3(c^2))((Subscript[e, 2]^2)-3 (Subscript[e, 1])(Subscript[e, 3])))//Simplify

Plot3D[(Subscript[e, 2]^2)-3 (Subscript[e, 1])(Subscript[e, 3])/.{(Subscript[e, 1])->u+v+(1-u-v),(Subscript[e, 2])->u v+(u+v)(1-u-v),(Subscript[e, 3])->u v(1-u-v)},{u,0,1},{v,0,1-u}]

Plot3D[(Subscript[e, 1]^2) (Subscript[e, 2])- 2(Subscript[e, 2]^2) -3 (Subscript[e, 1])(Subscript[e, 3])/.{(Subscript[e, 1])->u+v+(1-u-v),(Subscript[e, 2])->u v+(u+v)(1-u-v),(Subscript[e, 3])->u v(1-u-v)},{u,0,1},{v,0,1-u}]

αf=α/.Solve[f[(1/3),(1/3),(1/3),α]==-f[(1/2),(1/2),0,α],α][[1]]

αg=α/.Solve[Assuming[α>0,(g[(1/3),(1/3),(1/3),α]+g[(1/2),(1/2),0,α])//Simplify]==0,α][[1]]

αfT2=αf/.c->cT//Simplify;
errorT2s=f[(1/2),(1/2),0,αfT2]/.c->cT//Simplify;
ppT2s=p2[u,v,1-u-v,αfT2]/.c->cT;
(*Replacing variables u, v and w with symmetric polynomials (Subscript[e, 2]) and (Subscript[e, 3])*)ff=SymmetricReduction[f[u,v,w,αfT2]/.c->cT,{u,v,w},{1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
(*stacionary points in the interior of Ω*)
stacT2s=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacT2s=={});
(*stacionary points on the straight part of the boundary of Ω*)
stacBd1T2s=Select[Solve[D[f[u,0,1-u,αfT2]/.c->cT,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1T2s],i++,OK=OK&&(Abs[f[u,0,1-u,αfT2]/.c->cT]-errorT2s/.stacBd1T2s[[i]])≤0];
(*stacionary points on the curved parts of the boundary of Ω*)
stacBd2T2s=Select[Solve[D[f[u,u,1-2u,αfT2]/.c->cT,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2T2s],i++,OK=OK&&(Abs[f[u,u,1-2u,αfT2]/.c->cT]-errorT2s/.stacBd2T2s[[i]])≤0];
If[OK,Print["The best interpolant is induced by (Subscript[α, f])=",αfT2,"≐",N[αfT2]," with the simplified radial error ",errorT2s,"≐",N[errorT2s], ",
the max curvature ≐",N[GaussianCurvaturePoint[ppT2s,{u,v},{1/3,1/3}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppT2s,{u,v},{0,0}]]]]

αgT2=αg/.c->cT//Simplify;
errorT2=g[(1/2),(1/2),0,αgT2]/.c->cT//Simplify;
ppT2=p2[u,v,1-u-v,αgT2]/.c->cT;
(*Replacing variables u, v and w with symmetric polynomials (Subscript[e, 2]) and (Subscript[e, 3])*)ff=SymmetricReduction[f[u,v,w,αgT2]/.c->cT,{u,v,w},{1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
(*stacionary points in the interior of Ω*)
stacT2=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacT2=={});
(*stacionary points on the straight part of the boundary of Ω*)
stacBd1T2=Select[Solve[D[f[u,0,1-u,αgT2]/.c->cT,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1T2],i++,OK=OK&&Simplify[(Abs[g[u,0,1-u,αgT2]/.c->cT]-errorT2/.stacBd1T2[[i]])≤0]];
(*stacionary points on the curved parts of the boundary of Ω*)
stacBd2T2=Select[Solve[D[f[u,u,1-2u,αgT2]/.c->cT,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2T2],i++,OK=OK&&Simplify[(Abs[g[u,u,1-2u,αgT2]/.c->cT]-errorT2/.stacBd2T2[[i]])≤0]];
If[OK,Print["The best interpolant is induced by (Subscript[α, g])=",αgT2,"≐",N[αgT2],",
with the radial error ",errorT2,"≐",N[errorT2], "
the max curvature ≐",N[GaussianCurvaturePoint[ppT2,{u,v},{1/2,1/2}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppT2,{u,v},{0,0}]]]]

ppT2=p2[u,v,1-u-v,αgT2]/.c->cT;
Show[Graphics3D[{Blue,Opacity[.3],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesT[[i]].ppT2,{u,0,1},{v,0,1-u}],{i,4}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front]

αfO2=αf/.c->cO//Simplify;
errorO2s=f[(1/2),(1/2),0,αfO2]/.c->cO//Simplify;
ppO2s=p2[u,v,1-u-v,αfO2]/.c->cO;
(*Replacing variables u, v and w with symmetric polynomials (Subscript[e, 2]) and (Subscript[e, 3])*)ff=SymmetricReduction[f[u,v,w,αfO2]/.c->cO,{u,v,w},{1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
(*stacionary points in the interior of Ω*)
stacO2s=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacO2s=={});
(*stacionary points on the straight part of the boundary of Ω*)
stacBd1O2s=Select[Solve[D[f[u,0,1-u,αfO2]/.c->cO,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1O2s],i++,OK=OK&&(Abs[f[u,0,1-u,αfO2]/.c->cO]-errorO2s/.stacBd1O2s[[i]])≤0];
(*stacionary points on the curved parts of the boundary of Ω*)
stacBd2O2s=Select[Solve[D[f[u,u,1-2u,αfO2]/.c->cO,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2O2s],i++,OK=OK&&(Abs[f[u,u,1-2u,αfO2]/.c->cO]-errorO2s/.stacBd2O2s[[i]])≤0];
If[OK,Print["The best interpolant is induced by (Subscript[α, f])=",αfO2,"≐",N[αfO2]," with the simplified radial error ",errorO2s,"≐",N[errorO2s],",
and the max curvature ≐",N[GaussianCurvaturePoint[ppO2s,{u,v},{1/3,1/3}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppO2s,{u,v},{0,0}]]]]

αgO2=αg/.c->cO//Simplify;
errorO2=g[(1/2),(1/2),0,αgO2]/.c->cO//Simplify;
ppO2=p2[u,v,1-u-v,αgO2]/.c->cO;
(*Replacing variables u, v and w with symmetric polynomials (Subscript[e, 2]) and (Subscript[e, 3])*)ff=SymmetricReduction[f[u,v,w,αgO2]/.c->cO,{u,v,w},{1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
(*stacionary points in the interior of Ω*)
stacO2=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacO2=={});
(*stacionary points on the straight part of the boundary of Ω*)
stacBd1O2=Select[Solve[D[f[u,0,1-u,αgO2]/.c->cO,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1O2],i++,OK=OK&&Simplify[(Abs[g[u,0,1-u,αgO2]/.c->cO]-errorO2/.stacBd1O2[[i]])≤0]];
(*stacionary points on the curved parts of the boundary of Ω*)
stacBd2O2=Select[Solve[D[f[u,u,1-2u,αgO2]/.c->cO,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2O2],i++,OK=OK&&Simplify[(Abs[g[u,u,1-2u,αgO2]/.c->cO]-errorO2/.stacBd2O2[[i]])≤0]];
If[OK,Print["The best interpolant is induced by (Subscript[α, g])=",αgO2,"≐",N[αgO2]," with the radial error ",errorO2,"≐",N[errorO2],",
and the max curvature ≐",N[GaussianCurvaturePoint[ppO2,{u,v},{1/3,1/3}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppO2,{u,v},{0,0}]]]]

ppO2=p2[u,v,1-u-v,αgO2]/.c->cO;
Show[Graphics3D[{Blue,Opacity[.3],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesO[[i]].ppO2,{u,0,1},{v,0,1-u}],{i,8}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front]

αfI2=αf/.c->cI//Simplify;
errorI2s=f[(1/2),(1/2),0,αfI2]/.c->cI//Simplify;
ppI2s=p2[u,v,1-u-v,αfI2]/.c->cI;
(*Replacing variables u, v and w with symmetric polynomials (Subscript[e, 2]) and (Subscript[e, 3])*)ff=SymmetricReduction[f[u,v,w,αfI2]/.c->cI,{u,v,w},{1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
(*stacionary points in the interior of Ω*)
stacI2s=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacI2s=={});
(*stacionary points on the straight part of the boundary of Ω*)
stacBd1I2s=Select[Solve[D[f[u,0,1-u,αfI2]/.c->cI,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1I2s],i++,OK=OK&&Simplify[(Abs[f[u,0,1-u,αfI2]/.c->cI]-errorI2s/.stacBd1I2s[[i]])≤0]];
(*stacionary points on the curved parts of the boundary of Ω*)
stacBd2I2s=Select[Solve[D[f[u,u,1-2u,αfI2]/.c->cI,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2I2s],i++,OK=OK&&Simplify[(Abs[f[u,u,1-2u,αfI2]/.c->cI]-errorI2s/.stacBd2I2s[[i]])≤0]];
If[OK,Print["The best interpolant is induced by (Subscript[α, f])=",αfI2,"≐",N[αfI2]," with the simplified radial error ",errorI2s,"≐",N[errorI2s],",
and the max curvature ≐",N[GaussianCurvaturePoint[ppI2s,{u,v},{1/3,1/3}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppI2s,{u,v},{0,0}]]]]

αgI2=αg/.c->cI//Simplify;
errorI2=g[(1/2),(1/2),0,αgI2]/.c->cI//Simplify;
ppI2=p2[u,v,1-u-v,αgI2]/.c->cI;
(*Replacing variables u, v and w with symmetric polynomials (Subscript[e, 2]) and (Subscript[e, 3])*)ff=SymmetricReduction[f[u,v,w,αgI2]/.c->cI,{u,v,w},{1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
(*stacionary points in the interior of Ω*)
stacI2=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacI2=={});
(*stacionary points on the straight part of the boundary of Ω*)
stacBd1I2=Select[Solve[D[f[u,0,1-u,αgI2]/.c->cI,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1I2],i++,OK=OK&&Simplify[(Abs[g[u,0,1-u,αgI2]/.c->cI]-errorI2/.stacBd1I2[[i]])≤0]];
(*stacionary points on the curved parts of the boundary of Ω*)
stacBd2I2=Select[Solve[D[f[u,u,1-2u,αgI2]/.c->cI,u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2I2],i++,OK=OK&&FullSimplify[(Abs[g[u,u,1-2u,αgI2]/.c->cI]-errorI2/.stacBd2I2[[i]])≤0]];
If[OK,Print["The best interpolant is induced by (Subscript[α, g])=",αgI2,"≐",N[αgI2]," with the radial error ",errorI2,"≐",N[errorI2],",
and the max curvature ≐",N[GaussianCurvaturePoint[ppI2,{u,v},{1/3,1/3}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppI2,{u,v},{0,0}]]]]

ppI2=p2[u,v,1-u-v,αgI2]/.c->cI;
Show[Graphics3D[{Blue,Opacity[.3],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesI[[i]].ppI2,{u,0,1},{v,0,1-u}],{i,20}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front]

GraphicsRow[{Plot3D[g[u,v,1-u-v,αgT2]/.c->cT,{u,0,1},{v,0,1-u},ViewVector->{-1,0.5,3*errorT2}],Plot3D[g[u,v,1-u-v,αgO2]/.c->cO,{u,0,1},{v,0,1-u},ViewVector->{-1,0.5,3*errorO2}],Plot3D[g[u,v,1-u-v,αgI2]/.c->cI,{u,0,1},{v,0,1-u},ViewVector->{-1,0.5,3*errorI2}]}]

6. Cubic (Superscript[G, 1]) approximation

(Subscript[b, 3,0,0])=(Subscript[v, 0]); (Subscript[b, 0,3,0])=(Subscript[v, 1]); (Subscript[b, 0,0,3])=(Subscript[v, 2]);
(Subscript[b, 2,1,0])=α (Subscript[v, 0])+β (Subscript[v, 1]); (Subscript[b, 2,0,1])=α (Subscript[v, 0])+β (Subscript[v, 2]); (Subscript[b, 1,2,0])=α (Subscript[v, 1])+β (Subscript[v, 0]);
(Subscript[b, 0,2,1])=α (Subscript[v, 1])+β (Subscript[v, 2]); (Subscript[b, 1,0,2])=α (Subscript[v, 2])+β (Subscript[v, 0]); (Subscript[b, 0,1,2])=α (Subscript[v, 2])+β (Subscript[v, 1]);
(Subscript[b, 1,1,1])={0,0,γ};
p3[u_,v_,w_,α_,β_,γ_]=Sum[(Subscript[B, i,j,3-i-j])(Subscript[b, i,j,3-i-j]),{i,0,3},{j,0,3-i}]//Simplify;
f[u_,v_,w_,α_,β_,γ_]=p3[u,v,w,α,β,γ].p3[u,v,w,α,β,γ]-1;
g[u_,v_,w_,α_,β_,γ_]=(Sqrt[1+f[u,v,w,α,β,γ]])-1;

Solve[p3[0,0,1,α,β,γ].D[p3[u,v,1-u-v,α,β,γ],u]==0/.{u->0,v->0},α]

candidates=Solve[CoefficientList[Cross[(Subscript[b, 3,0,0]),(Subscript[b, 0,0,3])].Cross[D[p3[u,0,1-u,(1/2) (2-2 β+3 (c^2) β),β,γ],u],D[p3[u,v,1-u-v,(1/2) (2-2 β+3 (c^2) β),β,γ],v]/.v->0],u]==Table[0,{i,5}],{β,γ}]//Simplify;
candidates=Permute[ candidates, Cycles[{{1,2}}]];

For[k=1,k≤3,k++,{(Subscript[α, k]),(Subscript[β, k]),(Subscript[γ, k])}={(1/2) (2-2 β+3 (c^2) β),β,γ}/.candidates[[k]]//Simplify;Print["k=",k,": {(Subscript[α, k]),(Subscript[β, k]),(Subscript[γ, k])}=",{(Subscript[α, k]),(Subscript[β, k]),(Subscript[γ, k])}]]

(Subscript[a, 0])=24-324 (Subscript[e, 2]^2)+324 (Subscript[e, 3]);
(Subscript[a, 1])=-39+567 (Subscript[e, 2]^2)-324 (Subscript[e, 3])-972 (Subscript[e, 2]) (Subscript[e, 3]);
(Subscript[a, 2])=16-243 (Subscript[e, 2]^2)+972 (Subscript[e, 2]) (Subscript[e, 3])-729 (Subscript[e, 3]^2);

f[(1/3),(1/3),(1/3),(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]-f[u,v,w,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]-
((c^4)/(9 (((4-3 (c^2)))^2) (1-(c^2))))((Subscript[a, 2])((( (c^2)+((Subscript[a, 1])/(2(Subscript[a, 2])))))^2)-((Subscript[a, 1]^2)/(4(Subscript[a, 2])))+(Subscript[a, 0]))/.{(Subscript[e, 2])->u v+ u w+v w,(Subscript[e, 3])->u v w}/.w->1-u-v//Simplify

(*stacionary points in the interior of Ω*)
staca2=Solve[D[(Subscript[a, 2]),(Subscript[e, 2])]==0&&D[(Subscript[a, 2]),(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(staca2=={});
(*stacionary points on the straight part of the boundary of Ω*)
staca2Bd1=Select[Solve[D[(Subscript[a, 2])/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[staca2Bd1],i++,OK=OK&&((Subscript[a, 2])/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0}/.staca2Bd1[[i]])≥0];
(*stacionary points on the curved parts of the boundary of Ω*)
staca2Bd2=Select[Solve[D[(Subscript[a, 2])/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[staca2Bd2],i++,OK=OK&&((Subscript[a, 2])/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)}/.staca2Bd2[[i]])≥0];
OK

(*stacionary points in the interior of Ω*)
stacv=Solve[D[-2(Subscript[a, 2])-(Subscript[a, 1]),(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacv=={});
(*stacionary points on the straight part of the boundary of Ω*)
staca2Bd1=Select[Solve[D[-2(Subscript[a, 2])-(Subscript[a, 1])/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[staca2Bd1],i++,OK=OK&&(-2(Subscript[a, 2])-(Subscript[a, 1])/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0}/.staca2Bd1[[i]])≥0];
(*stacionary points on the curved parts of the boundary of Ω*)
staca2Bd2=Select[Solve[D[-2(Subscript[a, 2])-(Subscript[a, 1])/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[staca2Bd2],i++,OK=OK&&(-2(Subscript[a, 2])-(Subscript[a, 1])/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)}/.staca2Bd2[[i]])≥0];
OK

((Subscript[a, 2])((( (c^2)+((Subscript[a, 1])/(2(Subscript[a, 2])))))^2)-((Subscript[a, 1]^2)/(4(Subscript[a, 2])))+(Subscript[a, 0]))/.c->1//Simplify

OK=True;
(*1-729 (Subscript[e, 3]^2) has relevant stacionary points only on the curved parts of the boundary of Ω*)
stacb=Select[Solve[D[1-729 (Subscript[e, 3]^2)/.{(Subscript[e, 3])->(1-2 u) (u^2)},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacb],i++,OK=OK&&(1-729 (Subscript[e, 3]^2)/.{(Subscript[e, 3])->(1-2 u) (u^2)}/.stacb[[i]])≥0];
OK

Cross[D[p3[u,v,1-u-v,(Subscript[α, 2]),(Subscript[β, 2]),(Subscript[γ, 2])],u],D[p3[u,v,1-u-v,(Subscript[α, 2]),(Subscript[β, 2]),(Subscript[γ, 2])],v]]/.{u->0,v->0}//Simplify
Cross[D[p3[u,v,1-u-v,(Subscript[α, 3]),(Subscript[β, 3]),(Subscript[γ, 3])],u],D[p3[u,v,1-u-v,(Subscript[α, 3]),(Subscript[β, 3]),(Subscript[γ, 3])],v]]/.{u->(1/2),v->(1/2)}//Simplify

errorT3s=f[(1/3),(1/3),(1/3),(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cT//Simplify;
errorT3=g[(1/3),(1/3),(1/3),(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cT//Simplify;
ppT3=p3[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cT;
Print["The best interpolant is induced by 
{(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])}=",Simplify[{(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])}/.c->cT],
"
with the radial error ",errorT3,"≐",N[errorT3]," with the simplified radial error ",errorT3s,"≐",N[errorT3s],"
and the max curvature ", GaussianCurvaturePoint[ppT3,{u,v},{1/3,1/3}],"≐",N[GaussianCurvaturePoint[ppT3,{u,v},{1/3,1/3}]]," and the min curvature ",GaussianCurvaturePoint[ppT3,{u,v},{0,0}],"≐",N[GaussianCurvaturePoint[ppT3,{u,v},{0,0}]]]

ppT3=p3[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cT;
Show[Graphics3D[{Blue,Opacity[.3],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesT[[i]].ppT3,{u,0,1},{v,0,1-u}],{i,4}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front]

errorO3s=f[(1/3),(1/3),(1/3),(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cO//Simplify;
errorO3=g[(1/3),(1/3),(1/3),(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cO//Simplify;
ppO3=p3[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cO;
Print["The best interpolant is induced by 
{(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])}=",Simplify[{(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])}/.c->cO],
"
with the radial error ",errorO3,"≐",N[errorO3]," with the simplified radial error ",errorO3s,"≐",N[errorO3s],"
and the max curvature ", GaussianCurvaturePoint[ppO3,{u,v},{1/3,1/3}],"≐",N[GaussianCurvaturePoint[ppO3,{u,v},{1/3,1/3}]]," and the min curvature ",GaussianCurvaturePoint[ppO3,{u,v},{0,0}],"≐",N[GaussianCurvaturePoint[ppO3,{u,v},{0,0}]]]

ppO3=p3[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cO;
Show[Graphics3D[{Blue,Opacity[.3],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesO[[i]].ppO3,{u,0,1},{v,0,1-u}],{i,8}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front]

errorI3s=f[(1/3),(1/3),(1/3),(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cI//Simplify;
errorI3=g[(1/3),(1/3),(1/3),(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cI//Simplify;
ppI3=p3[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cI;
Print["The best interpolant is induced by 
{(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])}=",Simplify[{(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])}/.c->cI],
"
with the radial error ",errorI3,"≐",N[errorI3]," with the simplified radial error ",errorI3s,"≐",N[errorI3s],"
and the max curvature ", GaussianCurvaturePoint[ppI3,{u,v},{1/3,1/3}],"≐",N[GaussianCurvaturePoint[ppI3,{u,v},{1/3,1/3}]]," and the min curvature ",GaussianCurvaturePoint[ppI3,{u,v},{0,0}],"≐",N[GaussianCurvaturePoint[ppI3,{u,v},{0,0}]]]

ppI3=p3[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cI;
Show[Graphics3D[{Blue,Opacity[.3],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesI[[i]].ppI3,{u,0,1},{v,0,1-u}],{i,20}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front]

GraphicsRow[{Plot3D[g[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cT,{u,0,1},{v,0,1-u},ViewVector->{-1,0.5,3*errorT3}],Plot3D[g[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cO,{u,0,1},{v,0,1-u},ViewVector->{-1,0.5,3*errorO3}],Plot3D[g[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.c->cI,{u,0,1},{v,0,1-u},ViewVector->{-1,0.5,3*errorI3}]},ImageSize->Large]

ppT2a=p3[u,v,1-u-v,(Subscript[α, 2]),(Subscript[β, 2]),(Subscript[γ, 2])]/.c->cT;
ppT2b=p3[u,v,1-u-v,(Subscript[α, 3]),(Subscript[β, 3]),(Subscript[γ, 3])]/.c->cT;
GraphicsRow[{Show[Graphics3D[{Blue,Opacity[.1],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesT[[i]].ppT2a,{u,0,1},{v,0,1-u}],{i,4}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front],
Show[Graphics3D[{Blue,Opacity[.1],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesT[[i]].ppT2b,{u,0,1},{v,0,1-u}],{i,4}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front]}]
ppO2a=p3[u,v,1-u-v,(Subscript[α, 2]),(Subscript[β, 2]),(Subscript[γ, 2])]/.c->cO;
ppO2b=p3[u,v,1-u-v,(Subscript[α, 3]),(Subscript[β, 3]),(Subscript[γ, 3])]/.c->cO;
GraphicsRow[{Show[Graphics3D[{Blue,Opacity[.1],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesO[[i]].ppO2a,{u,0,1},{v,0,1-u}],{i,8}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front],
Show[Graphics3D[{Blue,Opacity[.1],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesO[[i]].ppO2b,{u,0,1},{v,0,1-u}],{i,8}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front]}]
ppI2a=p3[u,v,1-u-v,(Subscript[α, 2]),(Subscript[β, 2]),(Subscript[γ, 2])]/.c->cI;
ppI2b=p3[u,v,1-u-v,(Subscript[α, 3]),(Subscript[β, 3]),(Subscript[γ, 3])]/.c->cI;
GraphicsRow[{Show[Graphics3D[{Blue,Opacity[.1],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesI[[i]].ppI2a,{u,0,1},{v,0,1-u}],{i,20}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front],
Show[Graphics3D[{Blue,Opacity[.1],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesI[[i]].ppI2b,{u,0,1},{v,0,1-u}],{i,20}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front]}]

7. Cubic (Superscript[G, 2]) approximation

RR=Assuming[0<c<1,ReflectionMatrix[Cross[(Subscript[b, 0,3,0]),(Subscript[b, 0,0,3])]]//Simplify];
RR//MatrixForm

(Subscript[γ, 1,v])[t_]=p3[-t,v,1-(-t)-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])];
(Subscript[γ, 2,v])[t_]=RR.p3[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1])]/.{u-> t+((6 (c^2) (4-5 (c^2)+6 (c^4) (1-v) v))/((4-3 (c^2)) (4-4 (c^2)+3 (c^4) (1-v) v))) (t^2), v->v+((6 (c^2) (1-v)-4)/(4-3 (c^2))) t+((6 (c^4) (2+9 (c^4) (((1-v))^2) v-3(c^2) (1+ v-2 (v^2))))/((((4-3 (c^2)))^2) (4-4 (c^2)+3 (c^4) (1-v) v)))(t^2)};

(Subscript[γ, 1,v])[0]==(Subscript[γ, 2,v])[0]//Simplify
D[(Subscript[γ, 1,v])[t],t]==D[(Subscript[γ, 2,v])[t],t]/.{t->0}//Simplify
D[(Subscript[γ, 1,v])[t],{t,2}]==D[(Subscript[γ, 2,v])[t],{t,2}]/.{t->0}//Simplify

8. Quartic (Superscript[G, 1]) approximation

(Subscript[b, 4,0,0])=(Subscript[v, 0]); (Subscript[b, 0,4,0])=(Subscript[v, 1]); (Subscript[b, 0,0,4])=(Subscript[v, 2]);
(Subscript[b, 3,1,0])=α (Subscript[v, 0])+β (Subscript[v, 1]); (Subscript[b, 3,0,1])=α (Subscript[v, 0])+β (Subscript[v, 2]); (Subscript[b, 1,3,0])=α (Subscript[v, 1])+β (Subscript[v, 0]);
(Subscript[b, 0,3,1])=α (Subscript[v, 1])+β (Subscript[v, 2]); (Subscript[b, 1,0,3])=α (Subscript[v, 2])+β (Subscript[v, 0]); (Subscript[b, 0,1,3])=α (Subscript[v, 2])+β (Subscript[v, 1]);
(Subscript[b, 2,2,0])=γ((Subscript[v, 0])+(Subscript[v, 1]));(Subscript[b, 2,0,2])=γ((Subscript[v, 0])+(Subscript[v, 2]));(Subscript[b, 0,2,2])=γ((Subscript[v, 1])+(Subscript[v, 2]));
(Subscript[b, 2,1,1])=ζ (Subscript[v, 0])+ξ((Subscript[v, 1])+(Subscript[v, 2]));(Subscript[b, 1,2,1])=ζ (Subscript[v, 1])+ξ((Subscript[v, 0])+(Subscript[v, 2]));(Subscript[b, 1,1,2])=ζ (Subscript[v, 2])+ξ((Subscript[v, 0])+(Subscript[v, 1]));
p4[u_,v_,w_,α_,β_,γ_,ζ_,ξ_]=Sum[(Subscript[B, i,j,4-i-j])(Subscript[b, i,j,4-i-j]),{i,0,4},{j,0,4-i}]//Simplify;
f[u_,v_,w_,α_,β_,γ_,ζ_,ξ_]=p4[u,v,w,α,β,γ,ζ,ξ].p4[u,v,w,α,β,γ,ζ,ξ]-1;
g[u_,v_,w_,α_,β_,γ_,ζ_,ξ_]=(Sqrt[1+f[u,v,w,α,β,γ,ζ,ξ]])-1;

Solve[p4[0,0,1,α,β,γ,ζ,ξ].D[p4[u,v,1-u-v,α,β,γ,ζ,ξ],u]==0/.{u->0,v->0},α]

candidates=Solve[CoefficientList[Cross[(Subscript[b, 4,0,0]),(Subscript[b, 0,0,4])].Cross[D[p4[u,0,1-u,(1/2) (2-2 β+3 (c^2) β),β,γ,ζ,ξ],u],D[p4[u,v,1-u-v,(1/2) (2-2 β+3 (c^2) β),β,γ,ζ,ξ],v]/.v->0],u]==Table[0,{i,7}],{β,γ,ζ,ξ}]//Simplify;
Length[candidates]

D[p4[u,v,1-u-v,(1/2) (2-2 β+3 (c^2) β),β,γ,ζ,ξ],u]/.candidates[[1]]/.{u->0,v->0}//Simplify
D[p4[u,v,1-u-v,(1/2) (2-2 β+3 (c^2) β),β,γ,ζ,ξ],u]/.candidates[[2]]/.{u->0,v->0}//Simplify
D[p4[u,v,1-u-v,(1/2) (2-2 β+3 (c^2) β),β,γ,ζ,ξ],u]/.candidates[[3]]/.{u->(1/2),v->0}//Simplify
D[p4[u,v,1-u-v,(1/2) (2-2 β+3 (c^2) β),β,γ,ζ,ξ],u]/.candidates[[4]]/.{u->(1/2),v->0}//Simplify

For[k=1,k≤2,k++,{(Subscript[α, k]),(Subscript[β, k]),(Subscript[γ, k]),(Subscript[ζ, k]),(Subscript[ξ, k])}={(1/2) (2-2 β+3 (c^2) β),β,γ,ζ,ξ}/.candidates[[k+4]]/.Solve[(γ/.candidates[[k+4]])==γ,β][[1]]//Simplify;Print["k=",k,": {(Subscript[α, k]),(Subscript[β, k]),(Subscript[γ, k]),(Subscript[ζ, k]),(Subscript[ξ, k])}=",Collect[{(Subscript[α, k]),(Subscript[β, k]),(Subscript[γ, k]),(Subscript[ζ, k]),(Subscript[ξ, k])},γ]]]

For[k=1,k≤2,k++,(Subscript[f, k])[u_,v_,w_]=f[u,v,w,(Subscript[α, k]),(Subscript[β, k]),(Subscript[γ, k]),(Subscript[ζ, k]),(Subscript[ξ, k])]]
For[k=1,k≤2,k++,(Subscript[g, k])[u_,v_,w_]=(Sqrt[1+(Subscript[f, k])[u,v,w]])-1]

8.1. Tetrahedron

ff=SymmetricReduction[Coefficient[(Subscript[f, 1])[u,v,w]/.c->cT,γ,2],{u,v,w},{ 1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
staclc=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(staclc=={});
stacBd1lc=Select[Solve[D[ff/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1lc],i++,OK=OK&&(ff/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0}/.stacBd1lc[[i]])≥0];
stacBd2lc=Select[Solve[D[ff/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2lc],i++,OK=OK&&(ff/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)}/.stacBd2lc[[i]])≥0];
OK

ff=SymmetricReduction[Coefficient[(Subscript[f, 1])[u,v,w]/.c->cT,γ,2]+Coefficient[(Subscript[f, 1])[u,v,w]/.c->cT,γ,1],{u,v,w},{ 1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
stacv=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacv=={});
stacBd1v=Select[Solve[D[ff/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1v],i++,OK=OK&&(ff/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0}/.stacBd1v[[i]])≥0];
stacBd2v=Select[Solve[D[ff/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2v],i++,OK=OK&&(ff/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)}/.stacBd2v[[i]])≥0];
OK

γfT4=γ/.Solve[(Subscript[f, 1])[(1/2),(1/2),0]==-(Subscript[f, 1])[(1/3),(1/3),(1/3)]/.c->cT,γ,PositiveReals][[1]];
Plot3D[(Subscript[f, 1])[u,v,1-u-v]/.{c->cT,γ->γfT4},{u,0,1},{v,0,1-u}]
errorT4s=(Subscript[f, 1])[(1/3),(1/3),(1/3)]/.{c->cT,γ->γfT4}//Simplify;
ppT4s=p4[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1]),(Subscript[ζ, 1]),(Subscript[ξ, 1])]/.{c->cT,γ->γfT4};
(*Replacing variables u, v and w with symmetric polynomials (Subscript[e, 2]) and (Subscript[e, 3])*)ff=SymmetricReduction[(Subscript[f, 1])[u,v,w]/.{c->cT,γ->γfT4},{u,v,w},{1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
(*stacionary points in the interior of Ω*)
stacT4s=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacT4s=={});
(*stacionary points on the straight part of the boundary of Ω*)
stacBd1T4s=Select[Solve[D[(Subscript[f, 1])[u,0,1-u]/.{c->cT,γ->γfT4},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1T4s],i++,OK=OK&&Simplify[(Abs[(Subscript[f, 1])[u,0,1-u]/.{c->cT,γ->γfT4}]-errorT4s/.stacBd1T4s[[i]])≤0]];
(*stacionary points on the curved parts of the boundary of Ω*)
stacBd2T4s=Select[Solve[D[(Subscript[f, 1])[u,u,1-2u]/.{c->cT,γ->γfT4},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2T4s],i++,OK=OK&&Simplify[(Abs[(Subscript[f, 1])[u,u,1-2u]/.{c->cT,γ->γfT4}]-errorT4s/.stacBd2T4s[[i]])≤0]];
If[OK,Print["The best interpolant is induced by (Subscript[γ, f])=",γfT4,"≐",N[γfT4],",
with the simplified radial error ",errorT4s,"≐",N[errorT4s], ",
the max curvature ≐",N[GaussianCurvaturePoint[ppT4s,{u,v},{(1/3),(1/3)}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppT4s,{u,v},{(1/2),(1/2)}]]]]

(*The above interpolant ppT4s is better then (Subscript[f, 2])[.,.] if γ>((22229-216 (Sqrt[2291]))/7602)*)
(Subscript[f, 2])[(1/2),(1/2),0]/.c->cT//Simplify
Solve[(Subscript[f, 2])[(1/2),(1/2),0]==((Subscript[f, 1])[(1/2),(1/2),0]/.γ->γfT4)/.c->cT,γ]

(*The above interpolant ppT4s is better then (Subscript[f, 2])[.,.] if γ<((22229-216 (Sqrt[2291]))/7602)*)
(Subscript[f, 2])[(1/3),(1/3),(1/3)]/.c->cT//Simplify
Solve[(Subscript[f, 2])[(1/3),(1/3),(1/3)]==((Subscript[f, 1])[(1/3),(1/3),(1/3)]/.γ->γfT4)/.c->cT,γ]

(*The above interpolant ppT4s is better then (Subscript[f, 2])[.,.] if γ=((22229-216 (Sqrt[2291]))/7602)*)
((Subscript[f, 1])[(1/2),(1/2),0]/.{c->cT,γ->γfT4})-((Subscript[f, 2])[u,v,1-u-v]/.{c->cT,u->(1/2)-(1/2)(Sqrt[((27(Sqrt[2291])-979)/2649)]),v->0,γ->((22229-216 (Sqrt[2291]))/7602)})≥0

γgT4=γ/.Solve[(Subscript[g, 1])[(1/2),(1/2),0]==-(Subscript[g, 1])[(1/3),(1/3),(1/3)]/.c->cT,γ,PositiveReals][[1]];
Plot3D[(Subscript[g, 1])[u,v,1-u-v]/.{c->cT,γ->γgT4},{u,0,1},{v,0,1-u}]
errorT4=(Subscript[g, 1])[(1/3),(1/3),(1/3)]/.{c->cT,γ->γgT4}//Simplify;
ppT4=p4[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1]),(Subscript[ζ, 1]),(Subscript[ξ, 1])]/.{c->cT,γ->γgT4};
(*Replacing variables u, v and w with symmetric polynomials (Subscript[e, 2]) and (Subscript[e, 3])*)ff=SymmetricReduction[(Subscript[f, 1])[u,v,w]/.{c->cT,γ->γgT4},{u,v,w},{1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
(*stacionary points in the interior of Ω*)
stacT4=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacT4=={});
(*stacionary points on the straight part of the boundary of Ω*)
stacBd1T4=Select[Solve[D[(Subscript[f, 1])[u,0,1-u]/.{c->cT,γ->γgT4},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1T4],i++,OK=OK&&Simplify[(Abs[(Subscript[g, 1])[u,0,1-u]/.{c->cT,γ->γgT4}]-errorT4/.stacBd1T4[[i]])≤0]];
(*stacionary points on the curved parts of the boundary of Ω*)
stacBd2T4=Select[Solve[D[(Subscript[f, 1])[u,u,1-2u]/.{c->cT,γ->γgT4},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2T4],i++,OK=OK&&Simplify[(Abs[(Subscript[g, 1])[u,u,1-2u]/.{c->cT,γ->γfT4}]-errorT4/.stacBd2T4[[i]])≤0]];
If[OK,Print["The best interpolant is induced by (Subscript[γ, g])=",γgT4,"≐",N[γgT4],",
with the simplified radial error ",errorT4,"≐",N[errorT4], ",
the max curvature ≐",N[GaussianCurvaturePoint[ppT4,{u,v},{(1/3),(1/3)}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppT4,{u,v},{(1/2),(1/2)}]]]]

(*The above interpolant ppT4 is better then (Subscript[f, 2])[.,.] if γ>((-2773+5832 (Sqrt[3]))/4686)*)
(Subscript[g, 2])[(1/2),(1/2),0]/.c->cT//Simplify
Solve[(Subscript[g, 2])[(1/2),(1/2),0]==((Subscript[g, 1])[(1/2),(1/2),0]/.γ->γgT4)/.c->cT,γ]

(*The above interpolant ppT4s is better then (Subscript[f, 2])[.,.] if γ<((-2773+5832 (Sqrt[3]))/4686)*)
(Subscript[g, 2])[(1/3),(1/3),(1/3)]/.c->cT//Simplify
Solve[(Subscript[g, 2])[(1/3),(1/3),(1/3)]==((Subscript[g, 1])[(1/3),(1/3),(1/3)]/.γ->γgT4)/.c->cT,γ]

(*The above interpolant ppT4s is better then (Subscript[f, 2])[.,.] if γ=((-2773+5832 (Sqrt[3]))/4686)*)
((Subscript[g, 1])[(1/2),(1/2),0]/.{c->cT,γ->γgT4})-((Subscript[g, 2])[u,v,1-u-v]/.{c->cT,u->1/3,v->0,γ->((-2773+5832 (Sqrt[3]))/4686)})≥0//Simplify

ppT4=p4[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1]),(Subscript[ζ, 1]),(Subscript[ξ, 1])]/.{c->cT,γ->γgT4};
Show[Graphics3D[{Blue,Opacity[.3],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesT[[i]].ppT4,{u,0,1},{v,0,1-u}],{i,4}],PlotRange->All,Boxed->False,Axis->False]

8.2. Octahedron

ff=SymmetricReduction[Coefficient[(Subscript[f, 1])[u,v,1-u-v]/.c->cO,γ,2],{u,v,w},{ 1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
staclc=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(staclc=={});
stacBd1lc=Select[Solve[D[ff/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1lc],i++,OK=OK&&(ff/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0}/.stacBd1lc[[i]])≥0];
stacBd2lc=Select[Solve[D[ff/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2lc],i++,OK=OK&&(ff/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)}/.stacBd2lc[[i]])≥0];
OK

Coefficient[(Subscript[f, 1])[(1/2),(1/2),0]/.c->cO,γ,1]/(2Coefficient[(Subscript[f, 1])[(1/2),(1/2),0]/.c->cO,γ,2])//Simplify

Coefficient[(Subscript[f, 1])[(1/3),(1/3),(1/3)]/.c->cO,γ,1]/(2Coefficient[(Subscript[f, 1])[(1/3),(1/3),(1/3)]/.c->cO,γ,2])//Simplify

γfO4=γ/.Solve[(Subscript[f, 1])[(1/2),(1/2),0]==-(Subscript[f, 1])[(1/3),(1/3),(1/3)]/.c->cO,γ,PositiveReals][[1]];
Plot3D[(Subscript[f, 1])[u,v,1-u-v]/.{c->cO,γ->γfO4},{u,0,1},{v,0,1-u}]
errorO4s=(Subscript[f, 1])[(1/3),(1/3),(1/3)]/.{c->cO,γ->γfO4}//Simplify;
ppO4s=p4[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1]),(Subscript[ζ, 1]),(Subscript[ξ, 1])]/.{c->cO,γ->γfO4};
(*Replacing variables u, v and w with symmetric polynomials (Subscript[e, 2]) and (Subscript[e, 3])*)ff=SymmetricReduction[(Subscript[f, 1])[u,v,w]/.{c->cO,γ->γfO4},{u,v,w},{1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
(*stacionary points in the interior of Ω*)
stacO4s=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacO4s=={});
(*stacionary points on the straight part of the boundary of Ω*)
stacBd1O4s=Select[Solve[D[(Subscript[f, 1])[u,0,1-u]/.{c->cO,γ->γfO4},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1O4s],i++,OK=OK&&Simplify[(Abs[(Subscript[f, 1])[u,0,1-u]/.{c->cO,γ->γfO4}]-errorO4s/.stacBd1O4s[[i]])≤0]];
(*stacionary points on the curved parts of the boundary of Ω*)
stacBd2O4s=Select[Solve[D[(Subscript[f, 1])[u,u,1-2u]/.{c->cO,γ->γfO4},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2O4s],i++,OK=OK&&Simplify[(Abs[(Subscript[f, 1])[u,u,1-2u]/.{c->cO,γ->γfO4}]-errorO4s/.stacBd2O4s[[i]])≤0]];
If[OK,Print["The best interpolant is induced by (Subscript[γ, f])=",γfO4,"≐",N[γfO4],",
with the simplified radial error ",errorO4s,"≐",N[errorO4s], ",
the max curvature ≐",N[GaussianCurvaturePoint[ppO4s,{u,v},{(1/3),(1/3)}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppO4s,{u,v},{(1/2),(1/2)}]]]]

(*The above interpolant ppO4s is better then (Subscript[f, 2])[.,.] if γ>((9 (409-8 (Sqrt[974])))/1510)*)
(Subscript[f, 2])[(1/2),(1/2),0]/.c->cO//Simplify
Solve[(Subscript[f, 2])[(1/2),(1/2),0]==((Subscript[f, 1])[(1/2),(1/2),0]/.γ->γfO4)/.c->cO,γ]

(*The above interpolant ppO4s is better then (Subscript[f, 2])[.,.] if γ<((9 (409-8 (Sqrt[974])))/1510)*)
(Subscript[f, 2])[(1/3),(1/3),(1/3)]/.c->cO//Simplify
Solve[(Subscript[f, 2])[(1/3),(1/3),(1/3)]==((Subscript[f, 1])[(1/3),(1/3),(1/3)]/.γ->γfO4)/.c->cO,γ]

(*The above interpolant ppO4s is better then (Subscript[f, 2])[.,.] if γ=((9 (409-8 (Sqrt[974])))/1510)*)
((Subscript[f, 1])[(1/2),(1/2),0]/.{c->cO,γ->γfO4})-((Subscript[f, 2])[u,v,1-u-v]/.{c->cO,u->(1/2)-(1/2)(Sqrt[((27(Sqrt[974])-619)/3666)]),v->0,γ->((9 (409-8 (Sqrt[974])))/1510)})≥0

γgO4=γ/.Solve[(Subscript[g, 1])[(1/2),(1/2),0]==-(Subscript[g, 1])[(1/3),(1/3),(1/3)]/.c->cO,γ,PositiveReals][[1]];
Plot3D[(Subscript[f, 1])[u,v,1-u-v]/.{c->cO,γ->γgO4},{u,0,1},{v,0,1-u}]
errorO4=(Subscript[g, 1])[(1/3),(1/3),(1/3)]/.{c->cO,γ->γgO4}//Simplify;
ppO4=p4[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1]),(Subscript[ζ, 1]),(Subscript[ξ, 1])]/.{c->cO,γ->γgO4};
(*Replacing variables u, v and w with symmetric polynomials (Subscript[e, 2]) and (Subscript[e, 3])*)ff=SymmetricReduction[(Subscript[f, 1])[u,v,w]/.{c->cO,γ->γgO4},{u,v,w},{1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//FullSimplify;
(*stacionary points in the interior of Ω*)
stacO4=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacO4=={});
(*stacionary points on the straight part of the boundary of Ω*)
stacBd1O4=Select[Solve[D[(Subscript[f, 1])[u,0,1-u]/.{c->cO,γ->γgO4},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1O4],i++,OK=OK&&Simplify[(Abs[(Subscript[g, 1])[u,0,1-u]/.{c->cO,γ->γgO4}]-errorO4/.stacBd1O4[[i]])≤0]];
(*stacionary points on the curved parts of the boundary of Ω*)
stacBd2O4=Select[Solve[D[(Subscript[f, 1])[u,u,1-2u]/.{c->cO,γ->γgO4},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2O4],i++,OK=OK&&Simplify[(Abs[(Subscript[g, 1])[u,u,1-2u]/.{c->cO,γ->γgO4}]-errorO4/.stacBd2O4[[i]])≤0]];
If[OK,Print["The best interpolant is induced by (Subscript[γ, g])=",γgO4,"≐",N[γgO4],",
with the simplified radial error ",errorO4,"≐",N[errorO4], ",
the max curvature ≐",N[GaussianCurvaturePoint[ppO4,{u,v},{(1/3),(1/3)}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppO4,{u,v},{(1/2),(1/2)}]]]]

(*The above interpolant ppO4 is better then (Subscript[f, 2])[.,.] if γ>(3/538) (309-512 (Sqrt[3])+8 (Sqrt[6126+1512 (Sqrt[3])]))*)
(Subscript[g, 2])[(1/2),(1/2),0]/.c->cO//Simplify
Solve[(Subscript[g, 2])[(1/2),(1/2),0]==((Subscript[g, 1])[(1/2),(1/2),0]/.γ->γgO4)/.c->cO,γ]
Simplify[ToRadicals[γ/.%[[1]]]]

(*The above interpolant ppO4s is better then (Subscript[f, 2])[.,.] if γ<(3/538) (309-512 (Sqrt[3])+8 (Sqrt[6126+1512 (Sqrt[3])]))*)
(Subscript[g, 2])[(1/3),(1/3),(1/3)]/.c->cO//Simplify
Solve[(Subscript[g, 2])[(1/3),(1/3),(1/3)]==((Subscript[g, 1])[(1/3),(1/3),(1/3)]/.γ->γgO4)/.c->cO,γ]
Simplify[ToRadicals[γ/.%[[1]]]]

(*The above interpolant ppO4 is better then (Subscript[f, 2])[.,.] if γ=(3/538) (309-512 (Sqrt[3])+8 (Sqrt[6126+1512 (Sqrt[3])]))*)
((Subscript[g, 1])[(1/2),(1/2),0]/.{c->cO,γ->γgO4})-((Subscript[g, 2])[u,v,1-u-v]/.{c->cO,u->1/3,v->0,γ->(3/538) (309-512 (Sqrt[3])+8 (Sqrt[6126+1512 (Sqrt[3])]))})≥0//Simplify

ppO4=p4[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1]),(Subscript[ζ, 1]),(Subscript[ξ, 1])]/.{c->cO,γ->γgO4};
Show[Graphics3D[{Blue,Opacity[.3],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesO[[i]].ppO4,{u,0,1},{v,0,1-u}],{i,8}],PlotRange->All,Boxed->False,Axis->False]

8.3. Icosahedron

(* The method from the previous two cases does not work in the icosahedron case *)
γ0=((-1931-1100 (Sqrt[5])+18 (Sqrt[6 (556615+248877 (Sqrt[5]))]))/(6(5771+2508 (Sqrt[5]))));
γγ=γ/.Solve[(Subscript[f, 1])[(1/2),(1/2),0]==-(Subscript[f, 1])[(1/3),(1/3),(1/3)]/.c->cI,γ,PositiveReals][[1]]//Simplify;
γγ==γ0//FullSimplify
N[γ0]
Plot3D[(Subscript[f, 1])[u,v,1-u-v]/.{c->cI,γ->γγ},{u,0,1},{v,0,1-u}]

ff=SymmetricReduction[Coefficient[(Subscript[f, 1])[u,v,w]/.c->cI,γ,2],{u,v,w},{ 1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
staclc=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(staclc=={});
stacBd1lc=Select[Solve[D[ff/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1lc],i++,OK=OK&&(ff/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0}/.stacBd1lc[[i]])≥0];
stacBd2lc=Select[Solve[D[ff/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2lc],i++,OK=OK&&(ff/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)}/.stacBd2lc[[i]])≥0];
OK

ff=SymmetricReduction[Coefficient[(Subscript[f, 1])[u,v,w]/.c->cI,γ,2]+Coefficient[(Subscript[f, 1])[u,v,w]/.c->cI,γ,1],{u,v,w},{ 1, (Subscript[e, 2]), (Subscript[e, 3])}][[1]]//Simplify;
stacv=Solve[D[ff,(Subscript[e, 2])]==0&&D[ff,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])},PositiveReals];
OK=(stacv=={});
stacBd1v=Select[Solve[D[ff/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1&];
For[i=1,i≤Length[stacBd1v],i++,OK=OK&&(ff/.{(Subscript[e, 2])->(1-u) u,(Subscript[e, 3])->0}/.stacBd1v[[i]])≥0];
stacBd2v=Select[Solve[D[ff/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)},u]==0,{u},PositiveReals]//Simplify,0<(u/.#)<1/2&];
For[i=1,i≤Length[stacBd2v],i++,OK=OK&&(ff/.{(Subscript[e, 2])->2 (1-2 u) u+(u^2),(Subscript[e, 3])->(1-2 u) (u^2)}/.stacBd2v[[i]])≥0];
OK

ϵ=(10^(-4));
γl=(1/2);
γr=N[γγ];
While[γr-γl>ϵ,
γn=(1/2)(γl+γr);
fn=N[(Subscript[f, 1])[u,v,1-u-v]/.{c->cI,γ->γn}];
Mγ=0;mγ=0;
ExtremsLocal=Select[NSolve[D[fn,u]==0&&D[fn,v]==0,{u,v},PositiveReals],(u/.#)>0&&(v/.#)>0&&(u+v/.#)<1&];
ExtremsBoundary=Select[Solve[ D[fn/.v->0,u]==0,u,PositiveReals],(u/.#)≥0&&(u/.#)≤ (1/2)&];
For[i=1,i≤Length[ExtremsLocal],i++,Mγ=Max[Mγ,(fn/.ExtremsLocal[[i]])];mγ=Min[mγ,(fn/.ExtremsLocal[[i]])]];
For[i=1,i≤Length[ExtremsBoundary],i++,Mγ=Max[Mγ,(fn/.ExtremsLocal[[i]])];mγ=Min[mγ,(fn/.ExtremsLocal[[i]])]];
If[Mγ+mγ>0,γr=γn,γl=γn]]
Print["The parameter (Subscript[γ, f]) is on the interval [",γl,",",γr,"]. 
The max of the simplified radial error function is ",Mγ," and the min is ",mγ,"."]

ϵ=(10^(-4));
γl=(1/2);
γr=N[γγ];
While[γr-γl>ϵ,
γn=(1/2)(γl+γr);
fn=N[(Subscript[f, 1])[u,v,1-u-v]/.{c->cI,γ->γn}];
gn=N[(Subscript[g, 1])[u,v,1-u-v]/.{c->cI,γ->γn}];
Mγ=0;mγ=0;
ExtremsLocal=Select[NSolve[D[fn,u]==0&&D[fn,v]==0,{u,v},PositiveReals],(u/.#)>0&&(v/.#)>0&&(u+v/.#)<1&];
ExtremsBoundary=Select[Solve[ D[fn/.v->0,u]==0,u,PositiveReals],(u/.#)≥0&&(u/.#)≤ (1/2)&];
For[i=1,i≤Length[ExtremsLocal],i++,Mγ=Max[Mγ,(gn/.ExtremsLocal[[i]])];mγ=Min[mγ,(gn/.ExtremsLocal[[i]])]];
For[i=1,i≤Length[ExtremsBoundary],i++,Mγ=Max[Mγ,(gn/.ExtremsLocal[[i]])];mγ=Min[mγ,(gn/.ExtremsLocal[[i]])]];
If[Mγ+mγ>0,γr=γn,γl=γn]]
Print["The parameter (Subscript[γ, g]) is on the interval [",γl,",",γr,"]. 
The max of the radial error function is ",Mγ," and the min is ",mγ,"."]

g=SymmetricReduction[(Subscript[f, 1])[u,v,w],{u,v,w},{1,(Subscript[e, 2]),(Subscript[e, 3])}][[1]]//Simplify

ext=Solve[D[g,(Subscript[e, 2])]==0&&D[g,(Subscript[e, 3])]==0,{(Subscript[e, 2]),(Subscript[e, 3])}];
Length[ext]

ext[[1]]//Simplify

(36 (c^2) (-1+(c^2)) (((4-8 γ+(c^2) (-1+6 γ)))^2))(Subscript[e, 2])/.ext[[2]]/.c->cI//Simplify
(36 (c^2) (-1+(c^2)) (((4-8 γ+(c^2) (-1+6 γ)))^2))(Subscript[e, 2])/.ext[[3]]/.c->cI//Simplify

Coefficient[-5 (539+227 (Sqrt[5]))+4 (2125+913 (Sqrt[5])) γ-68 (99+43 (Sqrt[5])) (γ^2),γ,2]<0
Solve[-5 (539+227 (Sqrt[5]))+4 (2125+913 (Sqrt[5])) γ-68 (99+43 (Sqrt[5])) (γ^2)==0,γ,Reals]

(Subscript[γ, 1])=γ/.Solve[(Subscript[f, 1])[(1/3),(1/3),(1/3)]==0/.c->cI,γ,PositiveReals][[1]]//FullSimplify
(Subscript[γ, 2])=γ/.Solve[(Subscript[f, 1])[(1/2),(1/2),0]==0/.c->cI,γ,PositiveReals][[1]]

{N[(Subscript[γ, 1]),10],N[(Subscript[γ, 2]),10]}

GraphicsRow[{Plot[(Subscript[f, 1])[u,v,1-u-v]/.{c->cI,γ->(Subscript[γ, 1]),v->0},{u,0,1}],Plot[(Subscript[f, 1])[u,v,1-u-v]/.{c->cI,γ->(Subscript[γ, 1]),v->u},{u,0,1/2}]}]
GraphicsRow[{Plot[(Subscript[g, 1])[u,v,1-u-v]/.{c->cI,γ->(Subscript[γ, 1]),v->0},{u,0,1}],Plot[(Subscript[g, 1])[u,v,1-u-v]/.{c->cI,γ->(Subscript[γ, 1]),v->u},{u,0,1/2}]}]

GraphicsRow[{Plot[(Subscript[f, 1])[u,v,1-u-v]/.{c->cI,γ->(Subscript[γ, 2]),v->0},{u,0,1}],Plot[(Subscript[f, 1])[u,v,1-u-v]/.{c->cI,γ->(Subscript[γ, 2]),v->u},{u,0,1/2}]}]
GraphicsRow[{Plot[(Subscript[g, 1])[u,v,1-u-v]/.{c->cI,γ->(Subscript[γ, 2]),v->0},{u,0,1}],Plot[(Subscript[g, 1])[u,v,1-u-v]/.{c->cI,γ->(Subscript[γ, 2]),v->u},{u,0,1/2}]}]

({{}, {}})

Clear[h];h[u_,t_]:=(Subscript[f, 1])[u,u,1-2u]-(Subscript[f, 1])[(1/2),(1/2),0]/.{c->cI,γ->((Subscript[γ, 2])-(Subscript[γ, 1]))t+(Subscript[γ, 1])}//Simplify;
Select[Solve[Coefficient[h[u,t],t,0]==0,u,Reals],0<(u/.#)<(1/2)&]
Select[Solve[Coefficient[h[u,t],t,0]+Coefficient[h[u,t],t,1]==0,u,Reals],0<(u/.#)<(1/2)&]
Select[Solve[Coefficient[h[u,t],t,0]+Coefficient[h[u,t],t,1]+Coefficient[h[u,t],t,2]==0,u,Reals],0<(u/.#)<(1/2)&]

boundary=Solve[D[(Subscript[f, 1])[u,0,1-u]/.{c->cI},u]==0&&D[(Subscript[f, 1])[u,0,1-u]/.{c->cI},{u,2}]<0&&(Subscript[f, 1])[u,0,1-u]+(Subscript[f, 1])[(1/2),(1/2),0]==0/.{c->cI},{u,γ},PositiveReals]

(Subscript[f, 1])[(1/3),(1/3),(1/3) ]+(Subscript[f, 1])[(1/2),0,(1/2) ]>0/.{c->cI}/.boundary[[1]]//Simplify
(Subscript[f, 1])[(1/3),(1/3),(1/3) ]+(Subscript[f, 1])[(1/2),0,(1/2) ]>0/.{c->cI}/.boundary[[2]]//Simplify

median=Solve[D[(Subscript[f, 1])[u,u,1-2u]/.{c->cI},u]==0&&D[(Subscript[f, 1])[u,u,1-2u]/.{c->cI},{u,2}]<0&&(Subscript[f, 1])[u,u,1-2u]+(Subscript[f, 1])[(1/2),(1/2),0]==0/.{c->cI},{u,γ},PositiveReals]

γfI4=γ/.median[[1]];
errorI4s=(Subscript[f, 1])[(1/3),(1/3),(1/3)]/.{c->cI,γ->γfI4}//Simplify;
ppI4s=Simplify[N[p4[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1]),(Subscript[ζ, 1]),(Subscript[ξ, 1])]/.{c->cI,γ->γγI4}]];
Print["The best interpolant is induced by (Subscript[γ, f])=",γfI4,"≐",N[γfI4],",
with the simplified radial error ",errorI4s,"≐",N[errorI4s], ",
the max curvature ≐",N[GaussianCurvaturePoint[ppI4s,{u,v},{0.13728818302662513`,0.13728818302662513`}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppI4s,{u,v},{0,0}]]]

(*The above interpolant ppI4s is better then (Subscript[f, 2])[.,.] if γ>0.6475484005265687`*)res1=Solve[((Subscript[f, 1])[(1/2),(1/2),0]/.{c->cI,γ->γfI4})-((Subscript[f, 2])[(1/3),(2/3),0]/.{c->cI})==0,γ]//Simplify;
N[res1]
γu=γ/.res1[[1]];

(*The above interpolant ppI4s is better then (Subscript[f, 2])[.,.] if γ<0.6479216169802989`*)res2=Solve[((Subscript[f, 1])[(1/3),(1/3),(1/3)]/.{c->cI,γ->γfI4})-((Subscript[f, 2])[(1/3),(1/3),(1/3)]/.{c->cI})==0,γ]//Simplify;
N[res2]
γd=γ/.res2[[2]];

γd-γu>0//Simplify

boundary=Select[Solve[D[(Subscript[f, 1])[u,0,1-u]/.{c->cI},u]==0&&D[(Subscript[f, 1])[u,0,1-u]/.{c->cI},{u,2}]<0&&(Subscript[g, 1])[u,0,1-u]+(Subscript[g, 1])[(1/2),(1/2),0]==0/.{c->cI},{u,γ},PositiveReals],0<(u/.#)<(1/2)&]

(Subscript[g, 1])[(1/3),(1/3),(1/3) ]+(Subscript[g, 1])[(1/2),0,(1/2) ]>0/.{c->cI}/.boundary[[1]]//Simplify

median=Solve[D[(Subscript[f, 1])[u,u,1-2u]/.{c->cI},u]==0&&D[(Subscript[f, 1])[u,u,1-2u]/.{c->cI},{u,2}]<0&&(Subscript[g, 1])[u,u,1-2u]+(Subscript[g, 1])[(1/2),(1/2),0]==0/.{c->cI},{u,γ},PositiveReals]

γfI4=γ/.median[[1]];
errorI4s=(Subscript[f, 1])[(1/3),(1/3),(1/3)]/.{c->cI,γ->γfI4}//Simplify;
ppI4s=Simplify[N[p4[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1]),(Subscript[ζ, 1]),(Subscript[ξ, 1])]/.{c->cI,γ->γγI4}]];
Print["The best interpolant is induced by (Subscript[γ, f])=",γfI4,"≐",N[γfI4],",
with the simplified radial error ",errorI4s,"≐",N[errorI4s], ",
the max curvature ≐",N[GaussianCurvaturePoint[ppI4s,{u,v},{0.13728818302662513`,0.13728818302662513`}]]," and the min curvature ≐",N[GaussianCurvaturePoint[ppI4s,{u,v},{0,0}]]]

(*The above interpolant ppI4s is better then (Subscript[f, 2])[.,.] if γ>0.6475484005265687`*)res1=Solve[((Subscript[f, 1])[(1/2),(1/2),0]/.{c->cI,γ->γfI4})-((Subscript[f, 2])[(1/3),(2/3),0]/.{c->cI})==0,γ]//Simplify;
N[res1]
γu=γ/.res1[[1]];

(*The above interpolant ppI4s is better then (Subscript[f, 2])[.,.] if γ<0.6479216169802989`*)res2=Solve[((Subscript[f, 1])[(1/3),(1/3),(1/3)]/.{c->cI,γ->γfI4})-((Subscript[f, 2])[(1/3),(1/3),(1/3)]/.{c->cI})==0,γ]//Simplify;
N[res2]
γd=γ/.res2[[2]];

γd-γu>0//Simplify

ppI4=Simplify[N[p4[u,v,1-u-v,(Subscript[α, 1]),(Subscript[β, 1]),(Subscript[γ, 1]),(Subscript[ζ, 1]),(Subscript[ξ, 1])]/.{c->cI,γ->γfI4}]];
Show[Graphics3D[{Blue,Opacity[.3],Sphere[{0,0,0},1],Boxed->False}],
Table[ParametricPlot3D[SymmetriesI[[i]].ppI4,{u,0,1},{v,0,1-u}],{i,20}],PlotRange->All,Boxed->False,Axis->False,ViewPoint->Front]