(*Example from 591, March 7, 2019*)

(*Example of the algebraic boundary of a convex cone and its dual*)
f = (x+z)(x-z)^2-y^2z;
f1 = f/.{z->1};
K = RegionPlot[f1≥0&&1≥x,{x,-2,2},{y,-2,2}]
Kb =Show[K,ContourPlot[f1==0,{x,-2,2},{y,-2,2},ContourStyle->{Thick,Black}]]

(*Finding the dual variety to the algebraic boundary of K*)
g = (GroebnerBasis[{f,a-D[f,x],b-D[f,y],c-D[f,z]},{a,b,c},{x,y,z}])[[1]]

(*The algebraic boundary of the dual cone is defined by g and the hyperplane a+c=0 (which corresponds to the singular point on the boundary of K)*)
Kd = RegionPlot[{(g/.{c->1})<=0||(a≥-1&&a≤0&&a^2≥2b^2)},{a,-2,2},{b,-2,2}];
Kdb = Show[Kd,ContourPlot[{(g/.{c->1})==0,a+1==0},{a,-2,2},{b,-2,2},ContourStyle->Directive[{Thick,Black}]]]

(* To optimize ax+by over the intersection of K with the plane z=1, we solve the polynomial g for the variable c *)
sol = c/.Simplify[Solve[g==0,c]]

s = {-((2 (a^3)-9 a (b^2)+2 (Sqrt[((((a^2)+6 (b^2)))^3)]))/(27 (b^2))),((-2 (a^3)+9 a (b^2)+2 (Sqrt[((((a^2)+6 (b^2)))^3)]))/(27 (b^2)))};

Manipulate[
{Show[K,ContourPlot[a x + b y==-(((-2 (a^3)+9 a (b^2)+2 (Sqrt[(a^6)+18 (a^4) (b^2)+108 (a^2) (b^4)+216 (b^6)]))/(27 (b^2)))),{x,-2,2},{y,-2,2},ContourStyle->Red]],
Show[Kdb,ParametricPlot[{t{a/(((-2 (a^3)+9 a (b^2)+2 (Sqrt[(a^6)+18 (a^4) (b^2)+108 (a^2) (b^4)+216 (b^6)]))/(27 (b^2)))),b/(((-2 (a^3)+9 a (b^2)+2 (Sqrt[(a^6)+18 (a^4) (b^2)+108 (a^2) (b^4)+216 (b^6)]))/(27 (b^2))))}} ,{t,-10,10},PlotStyle->Directive[Dashed]],ListPlot[{{a/(((-2 (a^3)+9 a (b^2)+2 (Sqrt[(a^6)+18 (a^4) (b^2)+108 (a^2) (b^4)+216 (b^6)]))/(27 (b^2)))),b/(((-2 (a^3)+9 a (b^2)+2 (Sqrt[(a^6)+18 (a^4) (b^2)+108 (a^2) (b^4)+216 (b^6)]))/(27 (b^2))))}} ,PlotStyle->Directive[Red,PointSize[.02]]],ListPlot[{{a,b}} ,PlotStyle->Directive[Blue,PointSize[.02]]]]},{a,-2,2},{b,-2,2}]

(*We can use the same techniques to study a four-dimensional cone, which we view as the cone over a three-dimensional convex body*)
M = {{t,x,y},{x,t,z},{y,z,t}};
MatrixForm[M]

F = Det[M]
ContourPlot3D[(F/.{t->1})==0,{x,-1,1},{y,-1,1},{z,-1,1}]

(* To find the dual variety of the boundary of K, we use a Grobner basis calculation*)
g = (GroebnerBasis[{F,a-D[F,x],b-D[F,y],c-D[F,z],d-D[F,t]},{a,b,c,d},{x,y,z,t}])[[1]]

S = ContourPlot3D[(g/.{d->1})==0,{a,-1,1},{b,-1,1},{c,-1,1},PlotPoints->60]

(* The algebraic boundary of the dual cone is given by this quartic surface and four planes (corresponding to the four singular points on the boundary of K)*)
Show[S,ContourPlot3D[((a+b+c+d)/.{d->1})==0,{a,-1,1},{b,-1,1},{c,-1,1},PlotPoints->60,ContourStyle->{Opacity[.5],Red},Mesh->None]]

Show[S,ContourPlot3D[((a+b+c+d)/.{d->1})==0,{a,-1,1},{b,-1,1},{c,-1,1},PlotPoints->60,ContourStyle->{Opacity[.5],Red},Mesh->None],ContourPlot3D[((-a-b+c+d)/.{d->1})==0,{a,-1,1},{b,-1,1},{c,-1,1},PlotPoints->60,ContourStyle->{Opacity[.5],Red},Mesh->None],ContourPlot3D[((-a+b-c+d)/.{d->1})==0,{a,-1,1},{b,-1,1},{c,-1,1},PlotPoints->60,ContourStyle->{Opacity[.5],Red},Mesh->None],ContourPlot3D[((a-b-c+d)/.{d->1})==0,{a,-1,1},{b,-1,1},{c,-1,1},PlotPoints->60,ContourStyle->{Opacity[.5],Red},Mesh->None]]