Issue in r-interoff-lib.c

[gmangiap]

Dear McStas Team,

I was using Guide_anyshape_r, written by Peter Link, the head of the optic group of MLZ/FRM II. After some attempts to simulate a guide with a customized shape and with different m-values for the different faces, I found a bug in a library used by this component (r-interoff-lib.c) that makes Guide_anyshape_r use only the m-value defined for the first face in the .OFF file. To further investigate the bug (which seems to affect only the version 3), I have simulated a rectangular straight guide with the following dimensions: w1 = 30 mm, h1 = 180 mm, w2 = 20 mm, h2 = 30 mm, l = 15.0 m and different combinations of m-values for the four different walls (some of which are physically not very plausible, but nevertheless were used for testing purposes). The behavior of this guide has been simulated with Guide_gravity, Guide_anyshape_r, and (where possible) Guide_anyshape, in order to compare the results. As the two Guide_anyshape-based guides do not take into account gravitation accurately, all the tests have been carried out ignoring the gravitational field. In all tests a common initial seed was used, and 10**7 neutrons simulated. Furthermore, in all the tests, R0 = 0.99 and Qc = 0.02174 / Å were used. As virtual source I used one with 2.0 deg divergence in both directions, emitting neutrons from 1 Å to 21 Å with three Maxwellian distributions. The following table reports the different sets of m-, alpha- and W-values used for the simulations. Actually, the sets in positions 1, 3, 9 and 10 cannot be tested with the two Guide_anyshape based guides, but have been reported anyway because there are other examples that I will discuss below.

In the simulations two wavelength monitors were installed at the entrance and the exit of the guide, as well as two position sensitive detectors. Finally, two monitors measuring H- and V- divergences were positioned at the guide exit. All the data may be easily checked in the archive I have attached for the analysis by the McStas team. In the following table I have just reported the intensities for a quick discussion of the bug

As it is possible to see, the output of Guide_anyshape_r coincides with those of Guide_anyshape and Guide_gravity ONLY when all the four walls have the same m-, alpha- and W-values (examples n. 5, 6, 7, and 8). In examples n. 2 and 4 the output of Guide_anyshape_r differs from the results of Guide_gravity and coincides with the output of the examples 6 and 8, respectively, which use different sets of parameters. For example, for the case n.2 the intensity provided by Guide_anyshape_r is identical to the output of the case 6, where all the walls have m = 3. This is due to the fact that in the .OFF file used for the description of the guide in the example 2, the first face has m = 3:

OFF
#
8 8 0
0.01500000 0.09000000 0.00000000
-0.01500000 0.09000000 0.00000000
-0.01500000 -0.09000000 0.00000000
0.01500000 -0.09000000 0.00000000
0.0100000 0.01500000 15.00000000
-0.0100000 0.01500000 15.00000000
-0.0100000 -0.01500000 15.00000000
0.0100000 -0.01500000 15.00000000
3 0 1 4 3
3 1 4 5 3
3 1 2 6 4
3 1 5 6 4
3 2 3 7 3
3 2 6 7 3
3 0 3 4 4
3 3 4 7 4

If we switch the first and third line or even the first two lines with the third and fourth ones

OFF
#
8 8 0
0.01500000 0.09000000 0.00000000
-0.01500000 0.09000000 0.00000000
-0.01500000 -0.09000000 0.00000000
0.01500000 -0.09000000 0.00000000
0.0100000 0.01500000 15.00000000
-0.0100000 0.01500000 15.00000000
-0.0100000 -0.01500000 15.00000000
0.0100000 -0.01500000 15.00000000
3 1 2 6 4
3 1 5 6 4
3 0 1 4 3
3 1 4 5 3
3 2 3 7 3
3 2 6 7 3
3 0 3 4 4
3 3 4 7 4

we get the result reported in the row “2switched” which is identical to the output of lines 4 and 8. As said, the wrong results is due to the fact that Guide_Anyshape_r uses only the m-value reported in the first line of the section describing the faces in the .OFF file. For the same reason the output of Guide_Anyshape_r for the example n. 4 is identical to that of the example n.8 (and of 2switched). As hinted, the bug affects only the version 3 of McStas. I tried to redo the examples 2 and 2stwitched with the version 2.7.1 and the results are consistent:

The bug should be located in r-interoff-lib.c, in particular in the sub r_off_intersect_all, where it seems there is no trace of the faces’ indices in the block after line 894:

#else
    r_intersection intersect4[4];
    intersect4[0].time=-FLT_MAX;
    intersect4[1].time=FLT_MAX;
    intersect4[2].time=FLT_MAX;
    intersect4[3].time=FLT_MAX;
    int t_size=off_clip_3D_mod(intersect4, A, B,
        data->vtxArray, data->vtxSize, data->faceArray, data->faceSize, data->normalArray );
    if(t_size>0){
      int i=0;
      if (intersect4[0].time == -FLT_MAX) i=1;
      data->numintersect=t_size;
      if (t0) *t0 = intersect4[i].time;
      if (n0) *n0 = intersect4[i].normal;
      if (t3) *t3 = intersect4[i+1].time;
      if (n3) *n3 = intersect4[i+1].normal;
      /* should also return t[0].index and t[i].index as polygon ID */
      return t_size;
    }
#endif

which should be probably corrected as:

#else
    r_intersection intersect4[4];
    intersect4[0].time=-FLT_MAX;
    intersect4[1].time=FLT_MAX;
    intersect4[2].time=FLT_MAX;
    intersect4[3].time=FLT_MAX;
    int t_size=r_off_clip_3D_mod(intersect4, A, B,
        data->vtxArray, data->vtxSize, data->faceArray, data->faceSize, data->normalArray );
    if(t_size>0){
      int i=0;
      if (intersect4[0].time == -FLT_MAX) i=1;
      data->numintersect=t_size;
      if (t0) *t0 = intersect4[i].time;
      if (n0) *n0 = intersect4[i].normal;
      if (t3) *t3 = intersect4[i+1].time;
      if (n3) *n3 = intersect4[i+1].normal;

      if (faceindex0) *faceindex0 = intersect4[i].index;
      if (faceindex3) *faceindex3 = intersect4[i+1].index;

      /* should also return t[0].index and t[i].index as polygon ID */
      return t_size;
    }
#endif

In any case, since Guide_Anyshape_r allows to define only integer values for m, I modified it together with its library by adding the possibility of defining custom (and floating point) values for m, alpha and W for each face. I have tested this routine (called Guide_Anyshape_g) with all the examples reported in the tables, and the returned output seems to be consistent (including those of lines 9 and 10)

Finally, to further test Guide_Anyshape_g, I tried to simulate a guide of similar size to the one used above (w1 = 30 mm, h1 = 180 mm, w2 = 21 mm, h2 = 30 mm, l = 15.0 m), with different m values on the four walls and different along z every 5 meters, namely:

As also stressed above, the sets of values are hardly plausible in reality, but for testing purposes and in order to have a test with a strong reliability, the values have been used in Guide_gravity (simulating three consecutive 5m-guides with a negligible gap between them) and Guide_Anyshape_g, whose results are reported in the archive attached, and whose intensities are

The new files for Guide_Anyshape_g are also included in the archive, in case the McStas team judges them useful for extending the possibilities for the simulation of neutron guides.

Archive_with_Tests.tar.gz

Best regards,

Gaetano Mangiapia Helmholtz-Zentrum Hereon Forschungs-Neutronenquelle Heinz Maier-Leibnitz (FRM II) - Technische Universität München

[willend]

@gmangiap thank you for the issue and related patches.

I have decided to integrate your patches in the forthcoming 3.2 and 2.7.2 releases, as part of Guide_anyshape_r and r-interoff, woven together with the patches from @KBGrammer https://github.com/McStasMcXtrace/McCode/issues/1248

[gmangiap]

Dear Peter @willend,

thanks for the message and information. When you have time, please also have a look at the component Guide_Anyshape_g, which extends the possibility to customize the mirror properties of each face defined in the .OFF file

Best,

Gaetano Mangiapia Helmholtz-Zentrum Hereon Forschungs-Neutronenquelle Heinz Maier-Leibnitz (FRM II) - Technische Universität München

[willend]

@gmangiap Guide_anyshape_g is in fact what I am implementing, but under the name of Guide_anyshape_r. :-) Going forward, I don't see why one would want to use the old version.