- Fixed NB indexing

- Fixed an issue with setbin generated time binnings

ubfour.c

@@ -41,6 +41,35 @@ static MATRIX calcUVectorFromAngles(reflection r)
   vectorSet(u, 2, Cosd(om) * Sind(r.chi));
   return u;
 }
+/*--------------------------------------------------------------------------------------
+ * The code for this was directly lifted from rafnb from ILL, routine UBH in rafnbb.f. 
+ * I have yet to find a proper description of normal beam geometry calculations.
+ *--------------------------------------------------------------------------------------*/
+static MATRIX calcNBUVectorFromAngles(reflection r)
+{
+  MATRIX u;
+  double om, gamma, nu;
+
+  u = makeVector();
+  if (u == NULL) {
+    return NULL;
+  }     
+  /*
+  VAB(1,I) = SN(1)*CS(2)*CS(3) + SN(2)*(1.-CS(1)*CS(3))
+  VAB(2,I) = SN(1)*SN(2)*CS(3) - CS(2)*(1.-CS(1)*CS(3))
+  VAB(3,I) = SN(3)
+  1 = gamma, 2 = om, 3 = nu, CS = cos, SN = sin
+*/
+  gamma = r.gamma;
+  om = r.omnb;
+  nu = r.nu;
+  vectorSet(u, 0,
+		  Sind(gamma)*Cosd(om)*Cosd(nu) + Sind(om)*(1. - Cosd(gamma)*Cosd(nu)));
+  vectorSet(u, 1,
+		  Sind(gamma)*Sind(om)*Cosd(nu) - Cosd(om)*(1. - Cosd(gamma)*Cosd(nu)));
+  vectorSet(u, 2, Sind(nu));
+  return u;
+}
 
 /*--------------------------------------------------------------------------------------*/
 static MATRIX reflectionToHC(reflection r, MATRIX B)
@@ -153,6 +182,99 @@ MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1,
 
   return UB;
 }
+/*---------------------------------------------------------------------------------------*/
+MATRIX calcNBUBFromCellAndReflections(lattice direct, reflection r1,
+                                    reflection r2, int *errCode)
+{
+  MATRIX B, HT, UT, U, UB, HTT;
+  MATRIX u1, u2, h1, h2;
+  double ud[3];
+  int status;
+  reflection r;
+
+  *errCode = 1;
+
+  /*
+     calculate the B matrix and the HT matrix
+   */
+  B = mat_creat(3, 3, ZERO_MATRIX);
+  status = calculateBMatrix(direct, B);
+  if (status < 0) {
+    *errCode = status;
+    return NULL;
+  }
+  h1 = reflectionToHC(r1, B);
+  h2 = reflectionToHC(r2, B);
+  if (h1 == NULL || h2 == NULL) {
+    *errCode = UBNOMEMORY;
+    return NULL;
+  }
+  HT = matFromTwoVectors(h1, h2);
+  if (HT == NULL) {
+    *errCode = UBNOMEMORY;
+    return NULL;
+  }
+
+  /*
+     calculate U vectors and UT matrix
+   */
+  u1 = calcNBUVectorFromAngles(r1);
+  u2 = calcNBUVectorFromAngles(r2);
+  if (u1 == NULL || u2 == NULL) {
+    *errCode = UBNOMEMORY;
+    return NULL;
+  }
+  UT = matFromTwoVectors(u1, u2);
+  if (UT == NULL) {
+    *errCode = UBNOMEMORY;
+    return NULL;
+  }
+
+  /*
+     debugging output
+     printf("B-matrix\n");
+     mat_dump(B);
+     printf("HT-matrix\n");
+     mat_dump(HT);
+     printf("UT-matrix\n");
+     mat_dump(UT);
+   */
+
+  /*
+     UT = U * HT
+   */
+  HTT = mat_tran(HT);
+  if (HTT == NULL) {
+    *errCode = UBNOMEMORY;
+    return NULL;
+  }
+  U = mat_mul(UT, HTT);
+  if (U == NULL) {
+    *errCode = UBNOMEMORY;
+    return NULL;
+  }
+  UB = mat_mul(U, B);
+  if (UB == NULL) {
+    *errCode = UBNOMEMORY;
+  }
+
+  /*
+     clean up
+   */
+  killVector(h1);
+  killVector(h2);
+  mat_free(HT);
+  mat_free(HTT);
+
+  killVector(u1);
+  killVector(u2);
+  mat_free(UT);
+
+  mat_free(U);
+  mat_free(B);
+
+  return UB;
+}
 
 /*-----------------------------------------------------------------------------------*/
 static void storeReflection(reflection r, double two_theta_obs,
@@ -277,7 +399,7 @@ MATRIX calcUBFromThreeReflections(reflection r1, reflection r2,
                                   int *errCode)
 {
   MATRIX u1, u2, u3, HCHI, HI, HIINV, UB;
-  double det;
+  double det, z1[3];
 
   if (lambda <= .1) {
     *errCode = INVALID_LAMBDA;
@@ -292,7 +414,68 @@ MATRIX calcUBFromThreeReflections(reflection r1, reflection r2,
   uToScatteringVector(u1, r1.s2t / 2., lambda);
   uToScatteringVector(u2, r2.s2t / 2., lambda);
   uToScatteringVector(u3, r3.s2t / 2., lambda);
+  
+  HCHI = buildHCHIMatrix(u1, u2, u3);
+  HI = buildIndexMatrix(r1, r2, r3);
+  if (HCHI == NULL || HI == NULL) {
+    *errCode = UBNOMEMORY;
+    killVector(u1);
+    killVector(u2);
+    killVector(u3);
+    return NULL;
+  }
 
+  HIINV = mat_inv(HI);
+  if (HIINV == NULL) {
+    *errCode = NOTRIGHTHANDED;
+    killVector(u1);
+    killVector(u2);
+    killVector(u3);
+    mat_free(HI);
+    mat_free(HCHI);
+    return NULL;
+  }
+  UB = mat_mul(HCHI, HIINV);
+
+  det = mat_det(UB);
+  if (det < .0) {
+    mat_free(UB);
+    UB = NULL;
+    *errCode = NOTRIGHTHANDED;
+  }
+
+  killVector(u1);
+  killVector(u2);
+  killVector(u3);
+  mat_free(HI);
+  mat_free(HCHI);
+  mat_free(HIINV);
+
+  return UB;
+}
+/*-----------------------------------------------------------------------------*/
+MATRIX calcNBUBFromThreeReflections(reflection r1, reflection r2,
+                                  reflection r3, double lambda,
+                                  int *errCode)
+{
+  MATRIX u1, u2, u3, HCHI, HI, HIINV, UB;
+  double det, z1[3];
+
+  if (lambda <= .1) {
+    *errCode = INVALID_LAMBDA;
+    return NULL;
+  }
+
+  *errCode = 1;
+
+  z1FromNormalBeam(lambda,r1.gamma,r1.omnb, r1.nu,z1);
+  u1 = vectorToMatrix(z1);
+  z1FromNormalBeam(lambda,r2.gamma,r2.omnb, r2.nu,z1);
+  u2 = vectorToMatrix(z1);
+  z1FromNormalBeam(lambda,r3.gamma,r3.omnb, r3.nu,z1);
+  u3 = vectorToMatrix(z1);
+  
+  
   HCHI = buildHCHIMatrix(u1, u2, u3);
   HI = buildIndexMatrix(r1, r2, r3);
   if (HCHI == NULL || HI == NULL) {