The math is correct, or is in the right ballpark. The tool you used to calculate the FOV, does not do just the geometrical calculation, but also takes into account the HAM (hidden area mask) which may limit a diagonal FOV.
If you have a look at the HMDGDB database at, for example, Pimax 5k Plus page (Pimax 5K Plus (Normal FOV, Parallel Projection ON, 90Hz) | HMD Geometry Database), you can notice in the FOV visualization render that diagonal FOV is limited by the HAM in the peripheral corners and does not even reach the theoretical value the perfect frustum would have.
Just a note: The different FOVs (angles) are marked by triangles of different colors: yellow - vertical, red - overlap, magenta - peripheral, and cyan - diagonal.
Now coming back to the calculation for simple FOV of the shape of perfect frustum - i.e. one which does not have a HAM.
If I take your example, with vertical and horizontal FOVs being the same (i.e. 103°) and the frustum being symmetric, i.e. having the “left”, “right”, “top” and “bottom” FOVs being the same (i.e. 51.5°), calculating the diagonal FOV can be done relatively easily this way (I am using Python console to demonstrate the arithmetic operations):
IPython 8.4.0 -- An enhanced Interactive Python. Type '?' for help.
In [1]: from math import *
In [2]: hor_fov = radians(103)
In [3]: ver_fov = radians(103)
In [4]: half_tan_hor = tan(hor_fov/2)
In [5]: half_tan_ver = tan(ver_fov/2)
In [6]: half_tan_diag = sqrt(half_tan_hor**2 + half_tan_ver**2)
In [7]: half_diag = atan(half_tan_diag)
In [8]: diag_fov = 2 * half_diag
In [9]: degrees(diag_fov)
Out[9]: 121.28813765330132
Just to check the intermediate values:
In [10]: hor_fov, ver_fov, half_tan_hor, half_tan_ver, half_tan_diag, half_diag, diag_fov
Out[10]:
(1.7976891295541595,
1.7976891295541595,
1.2571722989189547,
1.2571722989189547,
1.7779101153709482,
1.0584386728311084,
2.1168773456622167)
You may notice that I am also using Pythagoras theorem in my calculation, but I am using it on tangent values (to get another tangent), and I can use it only because the frustum is symmetric and the diagonal “crosses” the frustum axis.
For non symmetric frustum the calculation is done purely in 3D Euclidean algebra by using vectors and operations on them.