(* example of creating degree 9 octahedrally symmetric body of constant width -- make sure to run rest of notebook before running this part*)

d=9;
invs = Invariants[[1]];
basis = SubspaceBasisFromInvariants[invs, d]/. {x^2 + y^2 + z^2 -> 1};
obasis = Orthogonalize[basis,multipolyinner]

coeffs =Join[Table[0, {i,1, Length[obasis]-1}],{1}]
f = obasis.coeffs

plotpoints = 40; 
e = MaxEpsilon[f]  
out =SimpleSurface[f,e, plotpoints]

(* Invariants and subspace finding functions *)

Invariants = {
(* oct *) {(* prim *){{(x^2)+(y^2)+(z^2),2},{(x^4)+(y^4)+(z^4),4},{(x^6)+(y^6)+(z^6),6}}, (* sec *) {{(x^5) (y^3) z-(x^3) (y^5) z-(x^5) y (z^3)+x (y^5) (z^3)+(x^3) y (z^5)-x (y^3) (z^5),9}}},
(* ico *){(* prim *){{(x^2)+(1/2) (1-(Sqrt[5])) x y+(y^2)+(1/2) (1-(Sqrt[5])) x z+(1/2) (1-(Sqrt[5])) y z+(z^2),2},{(x^6)+(1/2) (3-3 (Sqrt[5])) (x^5) y+(1/2) (15-3 (Sqrt[5])) (x^4) (y^2)+(5-4 (Sqrt[5])) (x^3) (y^3)+(1/2) (15-3 (Sqrt[5])) (x^2) (y^4)+(1/2) (3-3 (Sqrt[5])) x (y^5)+(y^6)+(1/2) (3-3 (Sqrt[5])) (x^5) z+(1/4) (45-15 (Sqrt[5])) (x^4) y z+(15-9 (Sqrt[5])) (x^3) (y^2) z+(15-9 (Sqrt[5])) (x^2) (y^3) z+(1/4) (45-15 (Sqrt[5])) x (y^4) z+(1/2) (3-3 (Sqrt[5])) (y^5) z+(1/2) (15-3 (Sqrt[5])) (x^4) (z^2)+(15-9 (Sqrt[5])) (x^3) y (z^2)+(1/4) (135-27 (Sqrt[5])) (x^2) (y^2) (z^2)+(15-9 (Sqrt[5])) x (y^3) (z^2)+(1/2) (15-3 (Sqrt[5])) (y^4) (z^2)+(5-4 (Sqrt[5])) (x^3) (z^3)+(15-9 (Sqrt[5])) (x^2) y (z^3)+(15-9 (Sqrt[5])) x (y^2) (z^3)+(5-4 (Sqrt[5])) (y^3) (z^3)+(1/2) (15-3 (Sqrt[5])) (x^2) (z^4)+(1/4) (45-15 (Sqrt[5])) x y (z^4)+(1/2) (15-3 (Sqrt[5])) (y^2) (z^4)+(1/2) (3-3 (Sqrt[5])) x (z^5)+(1/2) (3-3 (Sqrt[5])) y (z^5)+(z^6),6},{(x^10)+(1/2) (5-5 (Sqrt[5])) (x^9) y+(1/20) (405-99 (Sqrt[5])) (x^8) (y^2)+(1/10) (735-231 (Sqrt[5])) (x^6) (y^4)+(1/5) (315-189 (Sqrt[5])) (x^5) (y^5)+(1/10) (735-231 (Sqrt[5])) (x^4) (y^6)+(1/5) (150-102 (Sqrt[5])) (x^3) (y^7)+(1/20) (405-99 (Sqrt[5])) (x^2) (y^8)+(1/2) (5-5 (Sqrt[5])) x (y^9)+(y^10)+(1/2) (5-5 (Sqrt[5])) (x^9) z+(36-9 (Sqrt[5])) (x^8) y z+(1/5) (360-252 (Sqrt[5])) (x^7) (y^2) z+(1/5) (1050-378 (Sqrt[5])) (x^6) (y^3) z+(252-126 (Sqrt[5])) (x^5) (y^4) z+(252-126 (Sqrt[5])) (x^4) (y^5) z+(1/5) (1050-378 (Sqrt[5])) (x^3) (y^6) z+(1/5) (360-252 (Sqrt[5])) (x^2) (y^7) z+(36-9 (Sqrt[5])) x (y^8) z+(1/2) (5-5 (Sqrt[5])) (y^9) z+(1/5) (360-252 (Sqrt[5])) (x^7) y (z^2)+(1/5) (1575-441 (Sqrt[5])) (x^6) (y^2) (z^2)+(1/5) (1890-1134 (Sqrt[5])) (x^5) (y^3) (z^2)+(1/100) (150-102 (Sqrt[5])) (405-99 (Sqrt[5])) (x^15) (y^3) (z^2)+(630-189 (Sqrt[5])) (x^4) (y^4) (z^2)+(1/5) (1890-1134 (Sqrt[5])) (x^3) (y^5) (z^2)+(1/5) (1575-441 (Sqrt[5])) (x^2) (y^6) (z^2)+(1/5) (360-252 (Sqrt[5])) x (y^7) (z^2)+(1/20) (405-99 (Sqrt[5])) (y^8) (z^2)+(1/5) (150-102 (Sqrt[5])) (x^7) (z^3)+(1/5) (1890-1134 (Sqrt[5])) (x^5) (y^2) (z^3)+(630-294 (Sqrt[5])) (x^4) (y^3) (z^3)+(630-294 (Sqrt[5])) (x^3) (y^4) (z^3)+(1/5) (1890-1134 (Sqrt[5])) (x^2) (y^5) (z^3)+(1/5) (1050-378 (Sqrt[5])) x (y^6) (z^3)+(1/5) (150-102 (Sqrt[5])) (y^7) (z^3)+(252-126 (Sqrt[5])) (x^5) y (z^4)+(630-189 (Sqrt[5])) (x^4) (y^2) (z^4)+(630-294 (Sqrt[5])) (x^3) (y^3) (z^4)+(630-189 (Sqrt[5])) (x^2) (y^4) (z^4)+(252-126 (Sqrt[5])) x (y^5) (z^4)+(1/10) (735-231 (Sqrt[5])) (y^6) (z^4)+(1/5) (315-189 (Sqrt[5])) (x^5) (z^5)+(252-126 (Sqrt[5])) (x^4) y (z^5)+(1/5) (1890-1134 (Sqrt[5])) (x^3) (y^2) (z^5)+(1/5) (1890-1134 (Sqrt[5])) (x^2) (y^3) (z^5)+(252-126 (Sqrt[5])) x (y^4) (z^5)+(1/5) (315-189 (Sqrt[5])) (y^5) (z^5)+(1/10) (735-231 (Sqrt[5])) (x^4) (z^6)+(1/5) (1050-378 (Sqrt[5])) (x^3) y (z^6)+(1/5) (1575-441 (Sqrt[5])) (x^2) (y^2) (z^6)+(1/5) (1050-378 (Sqrt[5])) x (y^3) (z^6)+(1/10) (735-231 (Sqrt[5])) (y^4) (z^6)+(1/5) (150-102 (Sqrt[5])) (x^3) (z^7)+(1/5) (360-252 (Sqrt[5])) (x^2) y (z^7)+(1/50) (1050-378 (Sqrt[5])) (735-231 (Sqrt[5])) (x^12) y (z^7)+(1/5) (360-252 (Sqrt[5])) x (y^2) (z^7)+(1/5) (150-102 (Sqrt[5])) (y^3) (z^7)+(1/20) (405-99 (Sqrt[5])) (x^2) (z^8)+(36-9 (Sqrt[5])) x y (z^8)+(1/20) (405-99 (Sqrt[5])) (y^2) (z^8)+(1/2) (5-5 (Sqrt[5])) x (z^9)+(1/2) (5-5 (Sqrt[5])) y (z^9)+(z^10),10}}, (* sec *) {{(x^13) (y^2)+(1/4) (13-13 (Sqrt[5])) (x^12) (y^3)+(1/2) (41-23 (Sqrt[5])) (x^11) (y^4)+(1/2) (143-33 (Sqrt[5])) (x^10) (y^5)+(1/4) (237-141 (Sqrt[5])) (x^9) (y^6)+(6-27 (Sqrt[5])) (x^8) (y^7)+(-6+27 (Sqrt[5])) (x^7) (y^8)+(1/4) (-237+141 (Sqrt[5])) (x^6) (y^9)+(1/2) (-143+33 (Sqrt[5])) (x^5) (y^10)+(1/2) (-41+23 (Sqrt[5])) (x^4) (y^11)+(1/4) (-13+13 (Sqrt[5])) (x^3) (y^12)-(x^2) (y^13)+(1/4) (13-13 (Sqrt[5])) (x^12) (y^2) z+(63-15 (Sqrt[5])) (x^11) (y^3) z+(165-88 (Sqrt[5])) (x^10) (y^4) z+(1/2) (555-371 (Sqrt[5])) (x^9) (y^5) z+(1/4) (1827-315 (Sqrt[5])) (x^8) (y^6) z+(1/4) (-1827+315 (Sqrt[5])) (x^6) (y^8) z+(1/2) (-555+371 (Sqrt[5])) (x^5) (y^9) z+(-165+88 (Sqrt[5])) (x^4) (y^10) z+(-63+15 (Sqrt[5])) (x^3) (y^11) z+(1/4) (-13+13 (Sqrt[5])) (x^2) (y^12) z-(x^13) (z^2)+(1/4) (-13+13 (Sqrt[5])) (x^12) y (z^2)+(1/2) (187-143 (Sqrt[5])) (x^10) (y^3) (z^2)+(1/4) (2465-645 (Sqrt[5])) (x^9) (y^4) (z^2)+(702-378 (Sqrt[5])) (x^8) (y^5) (z^2)+(177-285 (Sqrt[5])) (x^7) (y^6) (z^2)+(-177+285 (Sqrt[5])) (x^6) (y^7) (z^2)+(-702+378 (Sqrt[5])) (x^5) (y^8) (z^2)+(1/4) (-2465+645 (Sqrt[5])) (x^4) (y^9) (z^2)+(1/2) (-187+143 (Sqrt[5])) (x^3) (y^10) (z^2)+(1/4) (13-13 (Sqrt[5])) x (y^12) (z^2)+(y^13) (z^2)+(1/4) (-13+13 (Sqrt[5])) (x^12) (z^3)+(-63+15 (Sqrt[5])) (x^11) y (z^3)+(1/2) (-187+143 (Sqrt[5])) (x^10) (y^2) (z^3)+(1/4) (1005-1305 (Sqrt[5])) (x^8) (y^4) (z^3)+(1086-126 (Sqrt[5])) (x^7) (y^5) (z^3)+(-1086+126 (Sqrt[5])) (x^5) (y^7) (z^3)+(1/4) (-1005+1305 (Sqrt[5])) (x^4) (y^8) (z^3)+(1/2) (187-143 (Sqrt[5])) (x^2) (y^10) (z^3)+(63-15 (Sqrt[5])) x (y^11) (z^3)+(1/4) (13-13 (Sqrt[5])) (y^12) (z^3)+(1/2) (-41+23 (Sqrt[5])) (x^11) (z^4)+(-165+88 (Sqrt[5])) (x^10) y (z^4)+(1/4) (-2465+645 (Sqrt[5])) (x^9) (y^2) (z^4)+(1/4) (-1005+1305 (Sqrt[5])) (x^8) (y^3) (z^4)+(1/2) (-441-483 (Sqrt[5])) (x^6) (y^5) (z^4)+(1/2) (441+483 (Sqrt[5])) (x^5) (y^6) (z^4)+(1/4) (1005-1305 (Sqrt[5])) (x^3) (y^8) (z^4)+(1/4) (2465-645 (Sqrt[5])) (x^2) (y^9) (z^4)+(165-88 (Sqrt[5])) x (y^10) (z^4)+(1/2) (41-23 (Sqrt[5])) (y^11) (z^4)+(1/2) (-143+33 (Sqrt[5])) (x^10) (z^5)+(1/2) (-555+371 (Sqrt[5])) (x^9) y (z^5)+(-702+378 (Sqrt[5])) (x^8) (y^2) (z^5)+(-1086+126 (Sqrt[5])) (x^7) (y^3) (z^5)+(1/2) (441+483 (Sqrt[5])) (x^6) (y^4) (z^5)+(1/2) (-441-483 (Sqrt[5])) (x^4) (y^6) (z^5)+(1086-126 (Sqrt[5])) (x^3) (y^7) (z^5)+(702-378 (Sqrt[5])) (x^2) (y^8) (z^5)+(1/2) (555-371 (Sqrt[5])) x (y^9) (z^5)+(1/2) (143-33 (Sqrt[5])) (y^10) (z^5)+(1/4) (-237+141 (Sqrt[5])) (x^9) (z^6)+(1/4) (-1827+315 (Sqrt[5])) (x^8) y (z^6)+(-177+285 (Sqrt[5])) (x^7) (y^2) (z^6)+(1/2) (-441-483 (Sqrt[5])) (x^5) (y^4) (z^6)+(1/2) (441+483 (Sqrt[5])) (x^4) (y^5) (z^6)+(177-285 (Sqrt[5])) (x^2) (y^7) (z^6)+(1/4) (1827-315 (Sqrt[5])) x (y^8) (z^6)+(1/4) (237-141 (Sqrt[5])) (y^9) (z^6)+(-6+27 (Sqrt[5])) (x^8) (z^7)+(177-285 (Sqrt[5])) (x^6) (y^2) (z^7)+(1086-126 (Sqrt[5])) (x^5) (y^3) (z^7)+(-1086+126 (Sqrt[5])) (x^3) (y^5) (z^7)+(-177+285 (Sqrt[5])) (x^2) (y^6) (z^7)+(6-27 (Sqrt[5])) (y^8) (z^7)+(6-27 (Sqrt[5])) (x^7) (z^8)+(1/4) (1827-315 (Sqrt[5])) (x^6) y (z^8)+(702-378 (Sqrt[5])) (x^5) (y^2) (z^8)+(1/4) (1005-1305 (Sqrt[5])) (x^4) (y^3) (z^8)+(1/4) (-1005+1305 (Sqrt[5])) (x^3) (y^4) (z^8)+(-702+378 (Sqrt[5])) (x^2) (y^5) (z^8)+(1/4) (-1827+315 (Sqrt[5])) x (y^6) (z^8)+(-6+27 (Sqrt[5])) (y^7) (z^8)+(1/4) (237-141 (Sqrt[5])) (x^6) (z^9)+(1/2) (555-371 (Sqrt[5])) (x^5) y (z^9)+(1/4) (2465-645 (Sqrt[5])) (x^4) (y^2) (z^9)+(1/4) (-2465+645 (Sqrt[5])) (x^2) (y^4) (z^9)+(1/2) (-555+371 (Sqrt[5])) x (y^5) (z^9)+(1/4) (-237+141 (Sqrt[5])) (y^6) (z^9)+(1/2) (143-33 (Sqrt[5])) (x^5) (z^10)+(165-88 (Sqrt[5])) (x^4) y (z^10)+(1/2) (187-143 (Sqrt[5])) (x^3) (y^2) (z^10)+(1/2) (-187+143 (Sqrt[5])) (x^2) (y^3) (z^10)+(-165+88 (Sqrt[5])) x (y^4) (z^10)+(1/2) (-143+33 (Sqrt[5])) (y^5) (z^10)+(1/2) (41-23 (Sqrt[5])) (x^4) (z^11)+(63-15 (Sqrt[5])) (x^3) y (z^11)+(-63+15 (Sqrt[5])) x (y^3) (z^11)+(1/2) (-41+23 (Sqrt[5])) (y^4) (z^11)+(1/4) (13-13 (Sqrt[5])) (x^3) (z^12)+(1/4) (13-13 (Sqrt[5])) (x^2) y (z^12)+(1/4) (-13+13 (Sqrt[5])) x (y^2) (z^12)+(1/4) (-13+13 (Sqrt[5])) (y^3) (z^12)+(x^2) (z^13)-(y^2) (z^13),15}}},
(* c4 *) {(* prim *){{z,1},{x^2+y^2,2}, {x^4+y^4, 4}}, (* sec *) {{x^3*y-x*y^3,4}}},
(* tet *) {(* prim *){{x^2 + y^2 + z^2,2},{x*y*z,3},{x^4 + y^4 + z^4, 4}}, (* sec *) {}}
};

SubspaceBasisFromInvariants[invariants_, d_] := (
primary = invariants[[1]];
secondary =invariants[[2]];
primarylist =Table[primary[[i,2]], {i,1,Length[primary]}];
secondarylist =Table[secondary[[i,2]], {i,1,Length[secondary]}];
out = {};
multilist = Join[Flatten[Table[Table[primarylist[[k]], {i, 1, d/primarylist[[k]]}], {k,1,Length[primarylist]}]], secondarylist];
basisindices = subs[multilist,d];
basisdegrees = DeleteDuplicates[Table[Table[multilist[[basisindices[[j, i]]]], {i, 1, Length[basisindices[[j]]]}], {j,1,Length[basisindices]}]];
basis  = {};
Do[b=1;
Do[b *= Select[Join[primary,secondary], #[[2]] == basisdegrees[[j,i]]& ][[1,1]], {i,1,Length[basisdegrees[[j]]]}];
basis = Append[basis, b];, {j,1,Length[basisdegrees]}];
basis
)

(* helper function for SubspaceBasisFromInvariants *)
subs[list_,n_]:=If[#1==n,#2,Nothing]&@@@Fold[Function[{x,a},Join[x,(If[#1+a[[2]]>n,Nothing,{#1+a[[2]],Append[#2,a[[1]]]}]&@@@x)]],{{0,{}}},Transpose@{Range@Length@list,list}]

(* writes a vector in terms of a basis *)
VectorInTermsofBasis[u_, basis_]:=(ClearAll[a]; ClearAll[b];Flatten[Table[b[i], {i, Length[basis]}] /.SolveAlways[u == Table[b[i], {i, Length[basis]}].basis, {x,y,z}]]);

(* inner product for orthogonalization *)

multipolyinner[f_,g_]:=(∫_({x,y,z}∈Sphere[{0,0,0}]))f*g

(* use this to test *)
innerprods[vs_, inner_] := Table[Table[inner[vs[[i]],vs[[j]]], {i,1,Length[vs]}], {j,1,Length[vs]}];

UnitBallQuotient[basis_] := ( q= basis /. {(x^2 + y^2 + z^2) -> 1};DeleteDuplicates[q]);

(* Functions to generate surfaces from support functions *)

(* support function for our cw body is 1 + ϵ*f(u) where u is a unit vector in S^n and f is an odd spherical harmonic invariant under G *)
SimpleFillmore[f_, ϵ_] := 1 + ϵ*f /. {x->#[[1]], y->#[[2]], z->#[[3]]} &;

(* parametrizes by unit vector the surface defined by h *)
SurfaceFromSupport[f_, ϵ_, u_]:=(
F[p_]:= f/. {x->p[[1]],y->p[[2]],z->p[[3]]};
h[v_] := 1 + ϵ*F[v];
H[v_] := Norm[v]*h[v/Norm[v]] /. Abs->Identity;
xs = Table[x[i], {i, 1, Length[u]}];
Grad[H[xs],xs] /. Table[xs[[i]]-> u[[i]], {i,1,Length[xs]}]
);

SurfaceFromSupport2D[f_, u_]:=(
F[p_]:= f/. {x->p[[1]],y->p[[2]]};
H[v_] := Norm[v]*F[v/Norm[v]] /. Abs->Identity;
xs = Table[x[i], {i, 1, Length[u]}];
Grad[H[xs],xs] /. Table[xs[[i]]-> u[[i]], {i,1,Length[xs]}]
);

SupportPlanes[f_, u_]:=(
F[p_]:= f/. {x->p[[1]],y->p[[2]]};
H[v_] := Norm[v]*F[v/Norm[v]] /. Abs->Identity;
xs = Table[x[i], {i, 1, Length[u]}];
pt =Grad[H[xs],xs] /. Table[xs[[i]]-> u[[i]], {i,1,Length[xs]}];
pt + s*{-u[[2]], u[[1]]}
);

(* spherical coordinate transforms in 2 and 3 dimensions *)
unitvs2[t_] := {Cos[t], Sin[t]};
unitvs3[t_, p_] := {Cos[t]*Sin[p], Sin[t]*Sin[p], Cos[p]};

(* wrapper to draw socw *)
SimpleSurface[f_, ϵ_, n_] :=ParametricPlot3D[SurfaceFromSupport[f, ϵ, unitvs3[t,p]], {t,0,2Pi}, {p,0,Pi},Mesh->{0}, PlotPoints->n]
UnwrappedSurface[f_, ϵ_, n_] :=Plot3D[Norm[SurfaceFromSupport[f, ϵ, unitvs3[t,p]]],{t,0, 2*Pi}, {p,0,Pi},Mesh->{40}];
SimpleSurfaceWidth[f_,ϵ_]:=Plot3D[SimpleFillmore[f, ϵ][unitvs3[t, p]] + SimpleFillmore[f, ϵ][-unitvs3[t, p]], {t, 0,2Pi},  {p,0,Pi}];

(* approximation of bounds on epsilon *)

MaxEpsilon[f_] :=(
ClearAll[ϵ];
F[p_]:= f/. {x->p[[1]],y->p[[2]],z->p[[3]]};
h[u_] := 1 + ϵ*F[u];
H[v_] := Norm[v]*h[v/Norm[v]] /. Abs->Identity;
epsi[{u1_,u2_,u3_},{v1_,v2_,v3_}] =(Abs[ϵ] /. Solve[H[{u1,u2,u3}+{v1,v2,v3}] - H[{u1,u2,u3}] - H[{v1,v2,v3}]==0,ϵ])[[1]];
NMinimize[{epsi[unitvs3[t1,p1], unitvs3[t2,p2]], 0<t1<Pi/2, 0 < t2 < Pi/2, Pi/2 < p1<Pi, Pi/2<p2<Pi}, {t1,t2,p1,p2}][[1]]
);

ConvexCheck[f_, q_] :=(
F[p_]:= f/. {x->p[[1]],y->p[[2]],z->p[[3]]};
h[u_] := 1 + q*F[u];
H[v_] := Norm[v]*h[v/Norm[v]] /. Abs->Identity;
Maximize[H[{u1,u2,u3} +{v1,v2,v3}] -H[{u1,u2,u3}] - H[{v1,v2,v3}], {u1,u2,u3,v1,v2,v3}]
);