Area of a (planar) quadrilateral using two opposite sides and four angles

I guessed the following area formula for a (planar) quadrilateral using two opposite sides and four angles:(When a=AB and c=CD,)

S = \frac{a^2}{2(\cot A + \cot B)} + \frac{c^2}{2(\cot C + \cot D)}

This formula can be applied to the quadrilateral satisfying A+B \ne \pi and C+D \ne \pi, including convex, concave, and even self-intersecting cases.

I’d like to know the applicable case , inapplicable case and reasons for the validity of those cases. Are there any known papers or more elegant ways to derive this formula?I guess there are many ways to explain those reasons. It greatly resembles S = \frac{1}{2}(ab \sin B + cd \sin D). It seems that b and d can be determined using a, c, A, B, C, D.

I just played around with your formula with this little program to check it on random quadrilaterals:

area[{a_, b_, c_, d_}] :=
  RegionMeasure[Polygon[{a, b, c, d}]];
formula[{a_, b_, c_, d_}] :=
  EuclideanDistance[a, b]^2/(
   2 (Cot[VectorAngle[b - a, d - a]] + 
      Cot[VectorAngle[c - b, a - b]])) +
   EuclideanDistance[d, c]^2/(
   2 (Cot[VectorAngle[b - c, d - c]] + Cot[VectorAngle[a - d, c - d]]));
Graphics[Polygon[rndmQuadrilateral =
   RandomInteger[{0, 10}, {4, 2}]], Frame -> True]
area[rndmQuadrilateral]
formula[rndmQuadrilateral]

This follows cleanly from the law of sines on each triangle formed by a diagonal. b and d are indeed determined (up to the constraint A+B+C+D=2π), and the formula breaks down when A+B=π or C+D=π because the diagonal length becomes indeterminate (parallel opposite sides, essentially a degenerate split).