Numerical problems (NaNs) with "get_theta_pol_phi_tor_from_two_points(x1_vec, x2_vec)"

[WolfgangSuttrop]

Hi Severin et al,

ECRad_Scenario.set_up_launch_from_imas() reads in line-of-sight (LOS) data from the IMAS ECE IDS in the form of two LOS end points, from which it calculates the poloidal and toroidal inclinations of the sightline:

ECRad_PyLib/src/ecrad_pylib/ECRad_Scenario.py

Lines 179 to 180 in f6322e9

	self['diagnostic']["theta_pol"][:], self['diagnostic']["phi_tor"][:] = \
	get_theta_pol_phi_tor_from_two_points(x1_vec, x2_vec)

Now for a real case with perpendicular incidence, AUG, where both LOS end points are at the same toroidal angle, "phi_tor" should be zero.
Example:

   x1_vec = np.array([-2.3048453154126185, -0.45846223813263876, 0.04097180441021919])
   x2_vec = np.array([-0.9807852804032304, -0.1950903220161284, 0.006681805010885])
   thetap, phit = get_theta_pol_phi_tor_from_two_points(x1_vec, x2_vec)

which results in phit== Nan (and not zero).
The reason is that numerical rounding errors make the np.arccos() argument in "get_theta_pol_phi_tor_from_two_points()" very slightly larger than one. I suggest to use a numerically more stable expression, such as:

phi_1 = np.arctan2(x1_vec[1], x1_vec[0])
phi_12 = np.arctan2(x1_vec[1]-x2_vec[1], x1_vec[0]-x2_vec[0])
self['diagnostic']["phi_tor"][:,ich] = np.rad2deg(phi_1 - phi_12)

Hope this helps ...

Greetings,
Wolfgang

[AreWeDreaming]

Thank you very much for finding this. On my machine I did not get NaNs but worse I got a 1.7 deg. angle instead of 0 for these two points that are on a line. The error gets smaller as you move away from precisely 0. The angle convention that I use has a different sign though so I added a -.
The change has been applied to plasma_math_tools.
Same question here as in #7, do you need a new conda release?