Guide_anyshape doesn't transport neutrons correctly with gravity

[KBGrammer]

Summary

The Guide_anyshape component does not currently function with gravity. The current implementation calculates an intersection time along a straight line and then propagates there with gravity enabled, which causes particles to fall through the bottom of the guide.

I have developed a fix and it is included as attachments here. The sections below describe the key changes to each file. The changes are made in interoff-lib.h, interoff-lib.c, and Guide_anyshape.comp. Several functions need the acceleration vector passed down to them. There are also new functions for setting up the quadratic to solve, calculating the intersection points in time, and handling the list of intersection points for the entire off geometry.

I have not copied this implementation over to the r-interoff-lib.c package.

I forked the main repository and pasted in my local copies of these three files. Unfortunately, my files had a bunch of "pragma" directives that seem to be missing from the main versions, probably because I took the files from my installation of McStas 3.1. If I did something wrong along the lines there, I'm happy to fix it. I also have the diffs for these three files comparing the 3.1 stock file to the patched files.

The fork is here: https://github.com/KBGrammer/mcstas_gravity

The commits are here: Guide_anyshape.comp https://github.com/McStasMcXtrace/McCode/commit/721af629743ea493888f3a78c6525c3f2236f41b interoff-lib.h https://github.com/McStasMcXtrace/McCode/commit/2ca3dc63eaafa1f25cf8d9d3d225d223bc0b7cba interoff-lib.c https://github.com/McStasMcXtrace/McCode/commit/00b27e00d4fbce06aff569b9e597bdcf9c2fff25

Guide_anyshape.comp

1. The local gravity vector needs to be determined. Here, it's hardcoded as (0, -GRAVITY, 0) and rotated according to the component rotation. (lines 127-132).
2. The call to off_intersect() must include the local gravity in addition to position and velocity. If mcgravitation is unset, this is safely (0, 0, 0). At line 134-135.
3. The call to PROP_DT must precede the calculation of n_dot_v, since the velocity changes in the case of gravity. This is safely moved to line 148, and can be safely done before calculating n_dot_v even without gravity.

interoff-lib.h

1. polygon struct: includes the "D" parameter of the plan for the polygon.
2. off_struct struct: includes an array of the D parameters
3. off_intersect_all: Function definition needs to include the acceleration vector.
4. off_intersect: Function definition needs to include the acceleration vector.
5. New function declaration p_to_quadratic, which takes a plane and the neutron trajectory information and returns the coefficients of a quadratic in time.
6. quadraticSolve, which solves that quadratic.

interoff-lib.c

1. Lines 426-520: New function off_clip_3D_mod_grav. Replicates functionality of off_clip_3D_mod in the case of gravity. The function p_to_quadratic and quadraticSolve are called at lines 469-471. "Huge" solutions are ignored (value of 1e36) and a list of the smallest positive intersections for all planes is returned.
2. The off_init function has been changed to include the DArray list, which stores a pre-computed list of the D parameters for each plane. Memory allocation is at line 721. Calculation of the D parameter is at lines 779-781, right after the determination of the surface normals.
3. Lines 920-941: New function p_to_quadratic, which stores the coefficients of the quadratic. This could all be moved to quadraticSolve, but they're kept separate so that quadraticSolve can solve a general quadratic if required.
4. Lines 943-981: New function quadraticSolve, which takes the coefficients returned by the previous function and solves for real roots and returns them as well as 0, 1, or 2 for the number of roots found. A "huge" value of 1.0e36 is returned as the "no root" value. The solver written to avoid cancellation issues from subtraction.
5. Line 988-1003: off_intersect_all needs the acceleration vector.
6. Lines 1008-1023: if(mcgravitation) cases which calls either off_clip_3D_mod_grav or off_clip_3D_mod in the OFF_LEGACY case.
7. Lines 1064-1078: if(mcgravitation) cases which calls either off_clip_3D_mod_grav or off_clip_3D_mod in the NOT OFF_LEGACY case.
8. Lines 1111-1114: off_intersect needs the acceleration vector.
9. Line 1143: the call to off_intersect here needs the acceleration vector, it's hardcoded as (0,0,0) in this spot.

Example files

This image shows a simple instrument with a tapered Guide_anyshape.comp before the patch with gravity turned on. Neutrons fall out the bottom in a very odd way.

This image shows the same instrument after the patch has been applied with gravity turned on. It's not shown, but without moving PROP_DT above n_dot_v, a very small number of neutrons end up lost outside the guide because the reflection sometimes continues to have a -y component and they fall out the bottom.

This image shows the same instrument after the patch has been applied but with gravity turned off to show that the component still works without gravity.

[willend]

@KBGrammer thanks a lot for your thoughts, work and patch!

You are absolutely right that Guide_anyshape does not work with gravitation.

The differences between the 2.x and 3.x versions of interoff-lib are small but subtle, but I should actually think that the 3.x version would mostly work with 2.x without adjustment.

I will test your solution and cherry-pick the right bits for inclusion in both forthcoming 2.x and 3.x releases.

Best, Peter Willendrup

[huegle]

I I remember correctly, he difference between the 2.7 .comp and the 3.1 .comp was mostly just the following:

if (mcgravitation) {
   Coords locgrav;
locgrav = rot_apply(ROT_A_COMP(mccompcurname), coords_set(0,-GRAVITY,0));  	//McStas2.7
   //~ locgrav = rot_apply(_comp->_rotation_absolute, coords_set(0,-GRAVITY,0)); 	 //McStas3.X

All the best, Thomas

[KBGrammer]

@willend I'm glad to hear it. I'm perfectly happy to answer questions about decisions I made in writing it and to help testing it.

Cheers, Kyle

[willend]

Hi @KBGrammer,

I am now finally in the process of merging in your patches for release with the forthcoming McStas 2.7.2 and 3.2 releases (and by symmetry McXtrace 1.7.1 and 3.1 - interoff is shared).

Is it correctly understood that if I simply use the modified functions with g ~ (0,0,0) everything should work "as usual"?

(Across the code base there is a number of different comps needing consequence-edits - x 2 for our two branches ;-) McStas:

Monitor_nD_noacc.comp
Monitor_nD.comp
Guide_anyshape_r.comp
Lens.comp
Single_crystal.comp
Incoherent.comp
Isotropic_Sqw.comp
PowderN.comp
Single_crystal_inelastic.comp
Single_magnetic_crystal.comp
Refractor.comp
Guide_anyshape.comp
monitor_nd-lib.h
monitor_nd_noacc-lib.h

plus a few for McXtrace)

[willend]

'heck I am pushing the edits now under the assumption that ~ (0,0,0) works - and we will see in tomorrow's test results what / if there is something to fix: http://new-nightly.mcstas.org

[KBGrammer]

@willend

Apologies for the delay in responding. I checked over my script that I used to make the example plots above. The difference between the plots titled "patched, gravity on" and "patched, gravity off" was just whether -g is set and the input was the same otherwise. Without -g, it sets the acceleration to (0,0,0). The quadraticSolve function checks to see if the equation it gets is linear and returns a value accordingly, so it will work with no gravity as well.

Thank you for getting this change included in a future update. I had thought there might be some incidental cases that would need updates, but I certainly wouldn't have found them all.

[willend]

The patches are on mccode-3, will try to get the same stuff done for master later / tomorrow.