search

gnudatalanguage / gdl

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

Problem with "/=" operator in 1-element arrays #1777

Closed jaymurthy closed 5 months ago

jaymurthy commented 5 months ago

I'm getting the wrong answer from adxy. I can't attach .sav files so I attach the fits file with the header which I read with mrdfits. fNCOB15.fits.gz

adxy,hdr,80.4112,42.1945,x,y print,x,y 180.78857 319.06283 where x and y should be 60,60. I've tracked it to a line in wcssph2xy. GDL> print,radeg 57.295780 GDL> phi /= radeg GDL> print,phi 88.402565 GDL> phi /= radeg GDL> print,phi 88.402565 Oddly it works when I type it in a new instance of gdl and it works in FL GDL> phi=88.402565 GDL> radeg=57.295780 GDL> phi /= radeg GDL> print,phi 1.54292

I'm using GDL - GNU Data Language, Version v1.0.2-534-g20545ef2 on Sonoma 14.3.1 with the M2 chip.

jaymurthy commented 5 months ago

SIMPLE = T / Written by IDL: Wed Jan 24 15:29:30 2024 BITPIX = -32 / Number of bits per data pixel NAXIS = 2 / Number of data axes NAXIS1 = 121 / NAXIS2 = 121 / EXTEND = T / FITS data may contain extensions DATE = '2024-01-24' / Creation UTC (CCCC-MM-DD) date of FITS header COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy COMMENT and Astrophysics', volume 376, page 359; bibcode 2001A&A...376..359H COMMENT Diffuse background CRVAL1 = 80.4112 / CRVAL2 = 42.1945 / CDELT1 = -0.100000 / CDELT2 = 0.100000 / CTYPE1 = 'GLON-TAN' / CTYPE2 = 'GLAT-TAN' / CRPIX1 = 60.0000 / CRPIX2 = 60.0000 / HISTORY Written on Wed Jan 24 15:29:30 2024 BUNIT = 'photon units' / END

jaymurthy commented 5 months ago

iTerm2 - gdl 2024-03-10 at 7 56 32 PM

alaingdl commented 5 months ago

OK, I got the same numbers with GDL hum, seems to be not so simple because on IDL I have : 

IDL> print, x,y
       58.999991       59.000000

Which version of the idlastro lib do you use ?

jaymurthy commented 5 months ago

I'm not entirely sure. I downloaded them from idlastro just before Wayne retired. I'll do a pull from github now. My adxy and wcss2phy are both from 2016. The IDL numbers are correct. I don't have IDL but I get the same answer with FL.

jaymurthy commented 5 months ago

OK. Just synced to the latest version with same result:

GDL> im=mrdfits('data/fNCOB15.fits.gz',0,hdr) MRDFITS: Image array (121,121) Type=Real*4 GDL> adxy,hdr,80.4112,42.1945,x,y & print,x,y 180.78857 319.06283 GDL>

$ fl

Fawlty Language 0.79.53 (macos amd64 m64) Copyright (c) 2000-2022, Lajos Foldy

                     https://www.flxpert.hu/fl/

THIS SOFTWARE COMES WITH ABSOLUTELY NO WARRANTY! USE AT YOUR OWN RISK!

Type 'help, /lib' for system routines info.

FL> im=mrdfits('data/fNCOB15.fits.gz',0,hdr) MRDFITS: Image array (121,121) Type=Real*4 FL> adxy,hdr,80.4112,42.1945,x,y & print,x,y 58.999991 59.000000

alaingdl commented 5 months ago

OK, FL & IDL give the same values for me, with the same idlastron than the one I use for GDL ...

jaymurthy commented 5 months ago

Yeah. It's because of these two lines in wcssph2xy.pro phi /= radeg theta /= radeg lines 345 and 346.

jaymurthy commented 5 months ago

OK. Here is the problem GDL> phi=[double(10)] GDL> print,phi/2 10.000000

alaingdl commented 5 months ago

Thanks, in fact I was converging to the same conclusions : 

GDL> a=[1] & a/=!pi & print, a
     0.318310
GDL> a=[1.] & a/=!pi & print, a
      1.00000
GDL> a=[10.] & a/=2 & print, a
      10.0000
GDL> a=[10] & a/=2 & print, a
      10
GDL> a=[10] & a/=2. & print, a
      5.00000
GDL> a=[10.] & a/=2. & print, a
      10.0000

Hum, /= operator is ... hum ... buggy !

alaingdl commented 5 months ago

Ok, after a fast check, only the /= is wrong. +=, -= & *= are OK. general relief :disappointed:

test_operators.pro.txt

alaingdl commented 5 months ago

And I check, it is OK when A when A is an array of size > 1

GillesDuvert commented 5 months ago

... So it is in the specialization for size 1 added in revision 687f1f7 of basic_op_div.cpp

GillesDuvert commented 5 months ago

i.e. in PR #1683

GillesDuvert commented 5 months ago

Solved. Thanks for the report and also the detective job! incidentally, this was only for /=, a 1 element variable, itself being a 'short int'.

jaymurthy commented 5 months ago

Great job!

slayoo commented 5 months ago

here's a test case draft: #1779

GillesDuvert commented 5 months ago

Thanks Sylwester! However I think a more general test for all operators (like in test_all_basic_functions.pro) is necessary. There are often more than 4 specialized code variants for each basic operator, and for each variant, internal switchs between 3 different code, depending on the operand type and size!

  1. scalar variable
  2. arrays < tpool_min_elts
  3. arrays > tpool_min_elts

I'll think about it.

GillesDuvert commented 5 months ago

Just produced such a general test. Using it, I find

GillesDuvert commented 5 months ago

actually, the OR and XOR code is wrong only in the subcase "scalar operand XOR something", i.e., one over 5 different cases.

GillesDuvert commented 5 months ago

Please refer to Discussion #1782 , the case is not clear.