How to find quadrilateral with maximum area?

Let the lengths of the sides of a quadrilateral equal a,b,c,d. What is the maximum area of such quadrilaterals? Wiki says that the quadrilateral with its maximum area must be cyclic and presents Brahmagupta's formula as an answer.

I am able to do that in Mathematica only for concrete values of a,b,c,d in such a way. The command

Maximize[{(a*d*Sin[x] + b*c*Sin[y])/2,
a^2 + b^2 - 2*a*d*Cos[x] == b^2 + c^2 - 2*b*c*Cos[y] &&x >= 0 &&
x <= Pi && y >= 0 && y <= Pi && {a, b, c, d} > 0 && a + b + c > d &&
b + c + d > a && c + d + a > b && d + a + b > c}, {x, y}]

, where the law of cosines the triangle inequality are used in the constraints , returns the input.

However, putting a = 1; b = 1; c = 1; d = 2; and a = 2; b = Pi/3; c = 3/2; d = 2; before the above code, one obtains correspondently {3/2, {x -> \[Pi]/2, y -> \[Pi]/2}} and {1/2 (4 Sin[2 ArcTan[Sqrt[(25 + 4 \[Pi])/(39 + 4 \[Pi])]]] + 1/2 \[Pi] Sin[2 ArcTan[Sqrt[(39 + 4 \[Pi])/(25 + 4 \[Pi])]]]), {x -> 2 ArcTan[Sqrt[(25 + 4 \[Pi])/(39 + 4 \[Pi])]], y -> 2 ArcTan[Sqrt[(39 + 4 \[Pi])/(25 + 4 \[Pi])]]}}.

Unfortunately, this is far away from Brahmagupta's formula and the general case and cyclic quadrilaterals.

How to derive the general case in Mathematica?

Addition. The above constraints are not complete and do not guarantee that a polygon is formed. Here is my next attempt. a = 1; b = 1; c = 1; d = 2;
NMaximize[{(a*d*Sin[x] + b*c*Sin[y])/2, a^2 + b^2 -
2*a*d*Cos[x] == b^2 + c^2 - 2*b*c*Cos[y] && a^2 + b^-2*a*b*Cos[z] ==
c^2 + d^2 - 2*c*d*Cos[2*Pi - x - y - z] &&
x >= 0 && x <= Pi && y >= 0 && y <= Pi && z >= 0 && z <= Pi &&
x + y + z <= 2*Pi && {a, b, c, d} > 0 && a + b + c > d &&
b + c + d > a && c + d + a > b && d + a + b > c}, {x, y, z}]

{2.75082, {x -> 1.41306, y -> 1.72849, z -> 2.15335}}

• I don't have access to that article. Nov 7 '21 at 7:13
• The result {3/2, {x -> \[Pi]/2, y -> \[Pi]/2}} for a = 1; b = 1; c = 1; d = 2; is not correct. Nov 7 '21 at 8:12