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]