search

AreWeDreaming / ECRad_PyLib

MIT License
1 stars 0 forks source link

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

Open WolfgangSuttrop opened 2 months ago

WolfgangSuttrop commented 2 months ago

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:

https://github.com/AreWeDreaming/ECRad_PyLib/blob/f6322e9261af3941eb48c983986c0a6e379b3c40/src/ecrad_pylib/ECRad_Scenario.py#L179-L180

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 commented 2 months ago

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?