irfu / irfu-matlab

Matlab routines to work with space data, particularly with MMS and Cluster/CAA data. Also some general plasma routines.
56 stars 46 forks source link

BUG report: irf_plot fails with certain inputs #56

Closed RibomBalt closed 3 years ago

RibomBalt commented 3 years ago

Step 1: Latest code?

Step 2: Describe your environment

Step 3: Describe the problem

With certain input, irf_plot fails to plot a figure and cast an error like this: (with original output translated to English)

错误使用 matlab.graphics.axis.Axes/set
// Error using matlab.graphics.axis.Axes/set
值必须是数值类型的 1x2 向量,其中第二个元素大于第一个元素,并且可能为 Inf
// Value should be numeric 1x2 vectors, with second element larger than the first element, maybe Inf
出错 irf_zoom>zoom_y_auto (line 278)
// Error irf_zoom>zoom_y_auto (line 278)
        set(h,'ylim',[ymin ymax]);

出错 irf_zoom (line 177)
// Error irf_zoom (line 177)
                    zoom_y_auto(h)

出错 irf_plot (line 231)
// Error irf_plot (line 231)
    if ~ishold(hca), irf_zoom(hca,'y'); end % automatic zoom only if hold is not on

Relevant code:

The code reproducing this error can be very simple. Any of the following lines could reproduce this problem.

irf_plot([1:5;(0.78:0.1:1.18)]')
irf_plot([1:5 ; 0.43:0.01:0.47]')
irf_plot([1:5 ; 0.05:0.01:0.09]')

After debugging, I guess this may be related to some floating-point arithmetic error stuff. irf_plot seems to automatically call irf_zoom to adjust the y-axis. During this process, when it reaches the following parts of irf_zoom.m:

262  dy=double(diffy)/4; % 1st approx
263  dy10power=10^(floor(log10(dy)));
264  dy1stcipher=floor(dy/dy10power);

dy1stcipher becomes zero, which is unexpected because dy and dy10power are equal in those 3 cases and this should return 1.

thomas-nilsson-irfu commented 3 years ago

Dear @RibomBalt, thank you for reporting this issue.

I have just made a commit in the devel branch of irfu-matlab which should hopefully help in fixing this issue (this commit will be merged into the master branch and included in the next release).

The problem is down to the floating-point accuracy of the parameters dy and dy10power, they should in your first example both have a value of exactly 0.1 but they are ever so slightly different due to how the values are stored binary. We have the dy being somewhat smaller than dy10power:

K>> dy-dy10power
ans =
    -2.775557561562891e-17 

This is within "eps()", https://se.mathworks.com/help/matlab/ref/eps.html, however it means that when we divide one by the other we get something like 0.999999... which when floor'ed goes to all the way to zero.

With best regards,

Thomas Nilsson