R0 prefactors inconsistent with those from ResLib

[js54434]

If we run the following scripts, we get different values for R0.

Here's the ResLib code which I ran in GNU Octave for a single value at (2,-0.18,0), 4 meV:

setup=MakeExp; 
setup.sample.a=4;  
setup.sample.b=4;
setup.sample.c=4;
setup.hcol=[50 80 50 120];  
setup.vcol = [300 300 300 300];

setup.sample.mosaic = 60.;
setup.sample.vmosaic = 1.;

setup.moncor = 1;
setup.method=0;
setup.mono.mosaic = 35.;
setup.ana.mosaic = 35.;

setup.orient1=[1 0 0];  
setup.orient2=[0 1 0];  

W1 = 4.;
H1 = 2.;
K1 = -0.18;
L1 = 0;

[R0,RMS] = ResMatS(H1,K1,L1,W1,setup)

Output:

R0 =  200.43
RMS =

   2.1775e+03  -3.4158e+02  -2.9914e+01   0.0000e+00
  -3.4158e+02   1.7115e+03   1.4749e+02   0.0000e+00
  -2.9914e+01   1.4749e+02   1.4298e+01   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00   9.0076e+01

Here's the neutronpy code:

import numpy as np
from neutronpy.resolution import Sample
from neutronpy.resolution import Instrument
efixed = 14.7
hcol = [50, 80, 50, 120]
vcol = [300, 300, 300, 300]
ana = 'PG(002)'
mono = 'PG(002)'

u = [1, 0, 0.]
v = [0, 1, 0.]
sample = Sample(4., 4., 4., 
                90, 90, 90, 
                mosaic=60.,
                u=u, v=v)
W = [4.]
H = [2.]
K = [-0.18]
L = [0.]

q = [H, K, L, W]

sample.vmosaic = 1.

EXP = Instrument(efixed=efixed, sample=sample, hcol=hcol, vcol=vcol, 
                ana=ana, mono=mono,
                moncor = 1,
                mono_mosaic=35., ana_mosaic=35.)
EXP.mono.dir = 1
EXP.sample.dir = -1
EXP.ana.dir = 1
EXP.calc_resolution(q)
RMS = EXP.RMS
R0 = EXP.R0

print(R0)
print(RMS)    

Output:

156.12974869026522
[[ 2177.45544664  -341.58347367     0.           -29.91424442]
 [ -341.58347367  1711.47074007     0.           147.48766283]
 [    0.             0.            90.07679538     0.        ]
 [  -29.91424442   147.48766283     0.            14.29832238]]

Also, their energy dependences seem to be different. Here are some plots I made of R0 vs energy, still at Q=(2,-0.18,0). The top plots are from neutronpy, and the bottom are from ResLib.

[pseudocubic]

I think there are two separate issues here.

First, in ResLib there is a factor of 1/sqrt(8*log(2)) intended to convert FWHM in arc min to statistical gaussian RMS width (standard deviation) for all of the mosaics in the variable CONVERT1, but in my translation I had removed this factor from CONVERT1 in favor of multiplying by it later in the calculation. Obviously, I got this right for the resolution matrix RM, but after fixing the second issue I discuss below it appears I still missed a factor somewhere for R0. In light of this it is probably better to just leave this factor in the original variable for now.

The second issue was in the calculation for sample mosaics, where R0 in ResLib was given by

R0_=R0_/sqrt((1+(q*etas)^2*RM_(4,4))*(1+(q*etasv)^2*RM_(2,2)));

When I translated this code it looks like I changed it such that RM(4,4) multiplied with the vertical mosaic etasv and RM(2,2) with the horizontal etas, essentially switching them. Checking the original paper (J. Appl. Phys. 42, 4736, (1971)) it appears I was wrong to make this adjustment since they define the new R0 as

where they assign the matrix indices 1: Qx, 2: Qy, 3: Qz, 4: W, as I do in my code, but different from ResLib which switches Qz and W.

I am working on pushing a fix for this issue today.