Determine and plot major and minor axes of ellipse

I have an element (El) constructed from a list of nodal coordinates (Nodes).

Note: The list of nodes can be arbitrarily long.

I then create a minimal area ellipse that encloses all of the node using a BoundingRegion
and determine the lengths of the major and minor radii of the ellipse. See example below.

I would like to extend this to include the calculation of the angles of the major and minor axes relative to the horizontal axis, as well as to plot these axes, preferably as dashed lines.

My code thus far is shown below.

(* Compute *)
Nodes = {{0, 0}, {48, 44}, {48, 60}, {0, 44}};
El = Polygon[Nodes];
ellipsoidBR = BoundingRegion[Nodes, "FastEllipsoid"];
ellipsoidB = RegionBoundary[ellipsoidBR];
area = RegionMeasure[ellipsoidBR];
Xc = RegionCentroid[ellipsoidBR];
a = RegionDistance[RegionBoundary[ellipsoidBR], Xc];
b = area/(Pi a );
(* Plot *)
Theta1 = 190;
Theta2 = 65;
Show[Graphics[{Opacity[0.20], ellipsoidBR}], 
 Graphics[{EdgeForm[Thick], Opacity[0.25], El}], 
 Graphics[{Thick, Circle[Xc, a]}], Graphics[{Thick, Circle[Xc, b]}], 
 Graphics[{PointSize[Large], Point[Xc]}], 
 Graphics[{PointSize[Large], Point[Nodes]}], 
 Graphics[{Thick, 
   Arrow[{Xc, {Xc[[1]] + a Cos[Theta1 Degree], 
      Xc[[2]] + a Sin[Theta1 Degree]}}]}], 
 Graphics[{Thick, 
   Arrow[{Xc, {Xc[[1]] + b Cos[Theta2 Degree], 
      Xc[[2]] + b Sin[Theta2 Degree]}}]}]]

Any help would be appreciated!

Nodes = {{0, 0}, {48, 44}, {48, 60}, {0, 44}};
ellipsoidBR = BoundingRegion[Nodes, "FastEllipsoid"];
center = ellipsoidBR[[1]];
{vals, vecs} = Eigensystem[ellipsoidBR[[2]]];
{a, b} = Sqrt[vals];
major = {center - a vecs[[1]], center + a vecs[[1]]};
minor = {center - b vecs[[2]], center + b vecs[[2]]};

Graphics[{{Opacity[0.20], ellipsoidBR},
          {Arrowheads[{-0.05, 0.05}], Arrow[major], Arrow[minor]}, 
          {Directive[AbsolutePointSize[6], Red], Point[Join[major, minor]]},
          {Directive[AbsolutePointSize[4], Purple], Point[Nodes]}}]

Nodes = {{0, 0}, {48, 44}, {48, 60}, {0, 44}};
ellipsoidBR = BoundingRegion[Nodes, "FastEllipsoid"];
RegionMember[ellipsoidBR, {x, y}]

Nodes = {{0, 0}, {48, 44}, {48, 60}, {0, 44}};
ellipsoidBR = BoundingRegion[Nodes, "FastEllipsoid"];
ellipsoidB = RegionBoundary[ellipsoidBR];
result = NMaximize[{EuclideanDistance[{x1, y1}, {x2, y2}], {x1, 
      y1} ∈ ellipsoidBR, {x2, y2} ∈ 
     ellipsoidBR}, {x1, y1, x2, y2}];
{pt1, pt2} = {{x1, y1}, {x2, y2}} /. result[[2]];
center = Mean[{pt1, pt2}];
result2 = 
  NMaximize[{t, 
    center + t*Cross[pt2 - center] ∈ ellipsoidB}, {t}];
pt3 = center + t*Cross[pt2 - center] /. result2[[2]];
pt4 = center - t*Cross[pt2 - center] /. result2[[2]];
Show[Graphics[{Opacity[0.20], ellipsoidBR}], 
 Graphics[{Arrow[{center, pt1}], Arrow[{center, pt2}], 
   Arrow[{center, pt3}], Arrow[{center, pt4}], Red, 
   Point[{pt1, pt2, pt3, pt4}]}]]