[BUG]: Source_gen does not generate wavelength distribution according to my input file

[Lomholy]

Describe the bug
When placing Source_gen with my input file spectrum_anew.txt, the final entry of the spectrume is loaded and used for all wavelengths, thereby creating a flat wavelength distribution instead of my input distribution.

The instrument and the spectrum can be found here:

src_gen_bug.txt

spectrum_anew.txt

OS and McStas/McXtrace version (please complete the following information):

• OS: MacOs
• McStas
• Version 3.7.4

[mads-bertelsen]

Thanks for the bug report, I can reproduce the issue and am looking into what is going on!

[Lomholy]

Thank you very much Mads!

I didn't fix the issue, but I found that by reversing the order of the wavelengths in my spectrum file, the problem disappears.
Here is the file if you want to test:

spectrum_clean.txt

Therefore I reckon that the issue is in that the Table_value function implicitly requires all inputs to be ascending in the 0 column, but doesn't say it explicitly in its doc string.

[mads-bertelsen]

I had just made the same discovery. There is some sense to have that assumption as it can speed up the algorithm. The used function does have this information in the header, but it should also be available in the component header.

/*******************************************************************************
* double Table_Value(t_Table Table, double X, long j)
*   ACTION: read column [j] of a single Table at row which 1st column is X
*   input   Table: table containing data.
*           X : data value in the first column (index 0)
*           j : index of column from which is extracted the Value (0:Columns-1)
*   return  Value = data[index for X][j] with linear interpolation
* Returns Value from the j-th column of Table corresponding to the
* X value for the 1st column (index 0)
* Tests are performed (within Table_Index) on indexes i,j to avoid errors
* NOTE: data should rather be monotonic, and evenly sampled.
*******************************************************************************/

[Lomholy]

I think I may be misunderstanding then 😅

Doesn't the text there merely say that the input file should be monotonic and evenly sampled, i.e it should be either ascending or descending?

Maybe I am misunderstanding the use of monotonic here...

[mads-bertelsen]

No it's true, it should specify monotonic increasing!

[mads-bertelsen]

read_table-lib calculates if the data is monotonic or not, and in that calculation it seems its definition of monotonic is that its monotonic and increasing, so the documentation makes sense within its own context. Nowhere does it specify that this flag means monotonic increasing. It also checks if the step size is constant. For your reversed data, it finds that it's not monotonic, but it does have constantstep, and that means it uses the constant-time lookup below under wrong assumptions.

Necessary actions:

• The documentation within read_table-lib and affected components should be updated

Possible actions, but needs full performance / accuracy check:

• Could add a check for monotonic in the Table_Value constant-time path, should require both constantstep and monotonic to be true

double Table_Value(t_Table Table, double X, long j)
{
  long   Index = -1;
  double X1=0, Y1=0, X2=0, Y2=0;
  double ret=0;

  if (X > Table.max_x) return Table_Index(Table,Table.rows-1  ,j);
  if (X < Table.min_x) return Table_Index(Table,0  ,j);

  // Use constant-time lookup when possible
  if(Table.constantstep) {
    Index = (long)floor(
              (X - Table.min_x) / (Table.max_x - Table.min_x) * (Table.rows-1));
    X1 = Table_Index(Table,Index-1,0);
    X2 = Table_Index(Table,Index  ,0);
  }
  // Use binary search on large, monotonic tables
  else if(Table.monotonic && Table.rows > 100) {
    long left = Table.min_x;
    long right = Table.max_x;

    while (!((X1 <= X) && (X < X2)) && (right - left > 1)) {
      Index = (left + right) / 2;

      X1 = Table_Index(Table, Index-1, 0);
      X2 = Table_Index(Table, Index,   0);

      if (X < X1) {
        right = Index;
      } else {
        left  = Index;
      }
    }
  }

  // Fall back to linear search, if no-one else has set X1, X2 correctly
  if (!((X1 <= X) && (X < X2))) {
    /* look for index surrounding X in the table -> Index */
    for (Index=1; Index <= Table.rows-1; Index++) {
        X1 = Table_Index(Table, Index-1,0);
        X2 = Table_Index(Table, Index  ,0);
        if ((X1 <= X) && (X < X2)) break;
      } /* end for Index */
  }

  Y1 = Table_Index(Table,Index-1, j);
  Y2 = Table_Index(Table,Index  , j);

#ifdef OPENACC
#define strcmp(a,b) str_comp(a,b)
#endif

  if (!strcmp(Table.method,"linear")) {
    ret = Table_Interp1d(X, X1,Y1, X2,Y2);
  }
  else if (!strcmp(Table.method,"nearest")) {
    ret = Table_Interp1d_nearest(X, X1,Y1, X2,Y2);
  }

#ifdef OPENACC
#ifdef strcmp
#undef strcmp
#endif
#endif

  return ret;
} /* end Table_Value */

[willend]

Guys, maybe a contribution in the form of a “good to know” page somewhere under the wiki could contain this? (Supplementing whatever is written in the specific comp. Trouble is that lots of comps in McStas use read_table-lib under the hood…)

[mads-bertelsen]

That makes sense, in the component builder category probably. It's not all components that use read_table-lib thats affected, only those that do lookups, here is the list:

Source_div_quasi
Source_gen
Pol_guide_mirror
Monochromator_curved
Pol_bender
Pol_mirror
Isotropic_Sqw
Pol_bender_tapering
Source_gen4
Guide_m
Transmission_V_polarisator
Monochromator_2foc
Guide_four_side
Mirror_Elliptic_Bispectral
Mirror_Curved_Bispectral
Virtual_mcnp_ss_Guide

[willend]

Nice with that list. On top comes the McXtrace side of things of course.