gnudatalanguage / gdl

GDL - GNU Data Language
GNU General Public License v2.0
276 stars 61 forks source link

ROT accuracy issue #1909

Open dinsmoro opened 4 days ago

dinsmoro commented 4 days ago

I noticed that ROT doesn't give the right outputs:

GDL> a=findgen(3,3)
GDL> print, a
      0.00000      1.00000      2.00000
      3.00000      4.00000      5.00000
      6.00000      7.00000      8.00000
GDL> print, rot(a,180.)
      5.00000      7.00000      6.00000
      1.00000      0.00000      3.00000
      1.00000      0.00000      0.00000

I tried some other older ROT functions from NASA personnel, but they performed the same which I suspect points to the poly_2d call (per we use the c++ internal POLY_2D code to do the job comment). I can't find that call in the repo, though.

This issue has popped up before, maybe system-dependent?: https://cow.physics.wisc.edu/~craigm/idl/archive/msg03150.html

And IDL on a separate computer does give the expected values:

IDL> a=findgen(3,3)
IDL> print, a
      0.00000      1.00000      2.00000
      3.00000      4.00000      5.00000
      6.00000      7.00000      8.00000
IDL> print, rot(a,180.)
      8.00000      7.00000      6.00000
      5.00000      4.00000      3.00000
      2.00000      1.00000      7.34788e-16
alaingdl commented 3 days ago

Thank you for this important bug report. I am sad that such an important bug was not detected before :( And no test in the testsuite. Sorry.

dinsmoro commented 3 days ago

Let me know if I can help - I'll take a stab at poly_2d if y'all need the help and you let me know where it is, but I can't guarantee anything since I don't know C++! But I can try!

alaingdl commented 9 hours ago

((I don't understand why we never realized that this code, which is big a complex, never gave problem before !))

OK, I found one problem in the math_fun_jmpg.cpp code, adding round() goes in the good direction. I have another trouble within permutation of i & j index, TBC

          DLong px = (DLong) round(P[0] + P[1] * (DDouble) i + P[2] * (DDouble) i);
          DLong py = (DLong) round(Q[0] + Q[1] * (DDouble) j + Q[2] * (DDouble) j);

I also realized that we have no systematic tests related to that part of code (TBD)