- Added brute force indexing support to ubcalc

- Added calculation of UB from 3 reflections to ubcalc
- Added calculation of lattice constants from UB to ubcalc
- Some fixes in stdscan in order to make the scripted scan work

cell.c

@@ -147,4 +147,31 @@ int calculateBMatrix(lattice direct, MATRIX B) {
   
   return 1;
 }
+/*--------------------------------------------------------------------------*/
+int cellFromUB(MATRIX UB, plattice direct){
+  MATRIX UBTRANS, GINV, G;
+
+  UBTRANS = mat_tran(UB);
+  if(UBTRANS == NULL){
+    return CELLNOMEMORY;
+  }
+  GINV = mat_mul(UBTRANS,UB);
+  if(GINV == NULL){
+    mat_free(UBTRANS);
+    return CELLNOMEMORY;
+  }
+  G = mat_inv(GINV);
+  if(G == NULL){
+    mat_free(UBTRANS);
+    mat_free(GINV);
+    return CELLNOMEMORY;
+  }
+  direct->a = sqrt(G[0][0]);
+  direct->b = sqrt(G[1][1]);
+  direct->c = sqrt(G[2][2]);
+  direct->alpha = Acosd(G[1][2]/(direct->b * direct->c));
+  direct->beta = Acosd(G[2][0]/(direct->a * direct->c));
+  direct->gamma = Acosd(G[0][1]/(direct->a * direct->c));
+  return 1;
+}
 

cell.h

@@ -14,7 +14,7 @@
  * error codes
  */
 #define REC_NO_VOLUME -100
-
+#define CELLNOMEMORY  -101
 /**
  * lattice parameters: either reciprocal or direct
  */
@@ -43,4 +43,11 @@
    * @return 1 on success, an negative error code else
    */
 int calculateBMatrix(lattice direct, MATRIX B);
+/**
+ * calculate the cell constants from a UB matrix
+ * @param UB The input UB matrix.
+ * @param direct A pointer to a structure holding the new cell constants
+ * @return 1 on success, an error c ode < 0 on failure
+ */
+int cellFromUB(MATRIX UB, plattice direct);
 #endif

fourlib.c

@@ -213,7 +213,7 @@ static void turnEquatorial(MATRIX z1, double *chi, double *phi){
 /*----------------------------------------------------------------------
 calculate d-spacing and theta from  z1
 -----------------------------------------------------------------------*/
-static int calcTheta(double lambda, MATRIX z1, double *d, double *theta){
+int calcTheta(double lambda, MATRIX z1, double *d, double *theta){
   double dstar, sintheta;
 
   dstar = sqrt(z1[0][0]*z1[0][0] + z1[1][0]*z1[1][0] + z1[2][0]*z1[2][0]);

fourlib.h

@@ -162,6 +162,12 @@ void chimat(MATRIX chim, double chi);
  * phi rotation matrix for phi into phim
 */
 void phimat(MATRIX phim, double phi);
+/**
+ * calcTheta calculates theta for a z1
+ * returns 1 on success, 0 else
+ */
+int calcTheta(double lambda, MATRIX z1, double *d, double *theta);
+ 
 #endif   
 
 

histsim.c

@@ -50,12 +50,6 @@
 #include "histsim.h"
 
 
-#define TESTVAL 128
-/*
-
-    define TESTVAL to set the simulated histogram to a fixed value for
-   testing
-*/   
   static int iSet = 0;
   static HistInt iSetVal = 0;
   static HistMode eHistMode;
@@ -86,6 +80,16 @@
       free(self->pPriv);
       self->pPriv = NewSIMCounter("HistoSim",fFail);
     }    
+    
+    /*
+      configured test value
+    */
+    status = StringDictGet(pOption,"testval",pData,131);
+    if(status)
+    {
+      iSet = 1;
+      iSetVal = atoi(pData);
+    }
 
     return 1;
   }  	
@@ -172,12 +176,7 @@
         SCWrite(pCon,"ERROR: histogram out of range",eError);
         return 0;
       }   
-#ifdef TESTVAL
+
-      for(ii = iStart; ii < iEnd; ii++)
-      {
-        lData[ii-iStart] = TESTVAL;
-      }      
-#else      
       if(iSet == 1)
       {
         for(ii = iStart; ii < iEnd; ii++)
@@ -192,7 +191,6 @@
             lData[ii-iStart] = random();
          }
       }
-#endif      
       return 1;
  }                              			
 /*------------------------------------------------------------------------*/

interface.c

@@ -42,7 +42,30 @@
 #include "fortify.h"
 #include "sics.h"
 #include "motor.h"
-
+/*=========================================================================
+  Empty driveable interface functions
+  ==========================================================================*/
+static int EmptyHalt(void *self){
+  return OKOK;
+}
+/*-----------------------------------------------------------------------*/
+static int EmptyLimits(void *self, float fVal, char *error, int errLen){
+  return 1;
+}
+/*-----------------------------------------------------------------------*/
+static long EmptyValue(void *self, SConnection *pCon, float fVal){
+  SCWrite(pCon,"WARNING: empty SetValue",eWarning);
+  return OKOK;
+}
+/*-----------------------------------------------------------------------*/
+static int EmptyStatus(void *self, SConnection *pCon){
+  return HWIdle;
+}
+/*------------------------------------------------------------------------*/
+static float EmptyGet(void *self, SConnection *pCon){
+  SCWrite(pCon,"WARNING: empty GetValue",eWarning);
+  return 555.55;
+}
 /*-------------------------------------------------------------------------*/
   pIDrivable CreateDrivableInterface(void)
   {
@@ -55,6 +78,11 @@
     }
     memset(pRes,0,sizeof(IDrivable));
     pRes->ID = DRIVEID;
+    pRes->Halt = EmptyHalt;
+    pRes->CheckLimits = EmptyLimits;
+    pRes->SetValue = EmptyValue;
+    pRes->CheckStatus = EmptyStatus;
+    pRes->GetValue = EmptyGet;
     return pRes;
   }
 /*-------------------------------------------------------------------------*/

stdscan.c

@@ -47,12 +47,15 @@ static char *fixExtension(char *filename)
       int iLen, iNum, iYear;
       char pNumText[10];
       CommandList *pCom = NULL;
+      char simName[255];
       
       /*
        make a simulated filename if in simulation mode
       */
-      if(pServ->simMode)
+      if(pServ->simMode){
-	return strdup("sim001001901.sim");
+	snprintf(simName,255,"%s/tmp/sim001001901.sim",getenv("HOME"));
+	return strdup(simName);
+      }
 
       pRes = makeFilename(pSics,pCon);
       if(pRes == NULL) 
@@ -756,6 +759,10 @@ int ScriptWriteHeader(pScanData self){
 }
 /*---------------------------------------------------------------------*/
 int ScriptPrepareScan(pScanData self){
+  /* configure counter */
+  SetCounterMode((pCounter)self->pCounterData,self->iMode);
+  SetCounterPreset((pCounter)self->pCounterData, self->fPreset);
+  self->iCounts = 0;
   return StandardScriptInvoke(self,"prepare");
 }
 /*-----------------------------------------------------------------------*/
@@ -764,7 +771,31 @@ int ScriptScanDrive(pScanData self, int iPoint){
 }
 /*------------------------------------------------------------------------*/
 int ScriptScanCount(pScanData self, int iPoint){
-  return StandardScriptInvokeWithPoint(self,"count",iPoint);
+  pDynString command = NULL;
+  int status;
+  char pNumber[50];
+
+  command = GetStandardInvocation(self,"count");
+  if(command == NULL){
+    return 0;
+  }
+  snprintf(pNumber,49,"%d",iPoint);
+  DynStringConcat(command,pNumber);
+  if(self->iMode == eTimer){
+    DynStringConcat(command," timer ");
+  } else {
+    DynStringConcat(command," monitor ");
+  }
+  snprintf(pNumber,49," %f ", self->fPreset);
+  DynStringConcat(command,pNumber);
+
+  status = InterpExecute(self->pSics, self->pCon,
+			 GetCharArray(command));
+  DeleteDynString(command);
+  if(status != 1) {
+    return 0;
+  }
+  return status;
 }
 /*-------------------------------------------------------------------------*/
 int ScriptScanCollect(pScanData self, int iPoint){

ubcalc.c

@@ -20,9 +20,13 @@
 /*---------------------------------------------------------------------*/
 typedef struct {
   pObjectDescriptor pDes;
+  pHKL hkl;
   lattice direct;
-  reflection r1, r2;
+  reflection r1, r2, r3;
   MATRIX UB;
+  double allowedDeviation;
+  int indexSearchLimit;
+  int maxSuggestions;
 } UBCALC, *pUBCALC;
 /*----------------------------------------------------------------------*/
 static void killUBCALC(void *pData){
@@ -53,10 +57,16 @@ static int SaveUBCalc(void *data, char *name, FILE *fd){
   fprintf(fd,"%s ref2 %f %f %f  %f %f %f %f\n",name,
 	  self->r2.h, self->r2.k, self->r2.l,
 	  self->r2.s2t, self->r2.om, self->r2.chi, self->r2.phi);
+  fprintf(fd,"%s ref3 %f %f %f  %f %f %f %f\n",name,
+	  self->r3.h, self->r3.k, self->r3.l,
+	  self->r3.s2t, self->r3.om, self->r3.chi, self->r3.phi);
+  fprintf(fd,"%s difftheta  %f\n", name, self->allowedDeviation);
+  fprintf(fd, "%s maxindex  %d\n", name, self->indexSearchLimit);
+  fprintf(fd ,"%s maxlist  %d\n", name, self->maxSuggestions);
   return 1;
 }
 /*---------------------------------------------------------------------*/
-static pUBCALC makeUBCALC(){
+static pUBCALC makeUBCALC(pHKL hkl){
   pUBCALC pNew = NULL;
 
   pNew = (pUBCALC)malloc(sizeof(UBCALC));
@@ -70,6 +80,10 @@ static pUBCALC makeUBCALC(){
     return NULL;
   }
   pNew->pDes->SaveStatus = SaveUBCalc;
+  pNew->hkl = hkl;
+  pNew->allowedDeviation = .3;
+  pNew->indexSearchLimit = 5;
+  pNew->maxSuggestions = 10;
   return pNew;
 }
 /*----------------------------------------------------------------------*/
@@ -77,17 +91,25 @@ int MakeUBCalc(SConnection *pCon, SicsInterp *pSics, void *pData,
 	       int argc, char *argv[]){
   pUBCALC pNew = NULL;
   int status;
+  pHKL hkl = NULL;
 
-  pNew = makeUBCALC();
+  if(argc < 3){
+    SCWrite(pCon,"ERROR: missing argument to MakeUBCalc: MakeUBCalc name hklobject",eError);
+    return 0;
+  }
+
+  hkl = FindCommandData(pSics,argv[2],"4-Circle-Calculus");
+  if(hkl == NULL){
+    SCWrite(pCon,"ERROR: HKL object not found or wrong type",eError);
+    return 0;
+  }
+
+  pNew = makeUBCALC(hkl);
   if(pNew == NULL){
     SCWrite(pCon,"ERROR: out of memory creating UBCALC",eError);
     return 0;
   }
-  if(argc > 1){
+  status = AddCommand(pSics,argv[1],UBCalcWrapper,killUBCALC,pNew);
-    status = AddCommand(pSics,argv[1],UBCalcWrapper,killUBCALC,pNew);
-  } else {
-    status = AddCommand(pSics,"ubcalc",UBCalcWrapper,killUBCALC,pNew);
-  }
   if(status != 1){
     SCWrite(pCon,"ERROR: failed to create duplicate UBCALC module",eError);
   }
@@ -319,18 +341,12 @@ static void listReflection(SConnection *pCon, char *name,
 }
 /*---------------------------------------------------------------------*/
 static int sendUBToHKL(SConnection *pCon, SicsInterp *pSics,
-			char *name, MATRIX UB){
+			pHKL hkl, MATRIX UB){
-  pHKL hkl = NULL;
-  char pBueffel[256];
   float ub[9];
   int i;
 
-  hkl = FindCommandData(pSics,name,"4-Circle-Calculus");
+  assert(hkl != NULL);
-  if(hkl == NULL){
+
-    snprintf(pBueffel,255,"ERROR: %s not found or of wrong type",name);
-    SCWrite(pCon,pBueffel,eError);
-    return 0;
-  }
   for(i = 0; i < 3; i++){
     ub[i]   = UB[0][i];
     ub[i+3] = UB[1][i];
@@ -340,6 +356,198 @@ static int sendUBToHKL(SConnection *pCon, SicsInterp *pSics,
   SCSendOK(pCon);
   return 1;
 } 
+/*---------------------------------------------------------------------*/
+static int setUBCalcParameters(pUBCALC self, SConnection *pCon, char *name, char *value){
+
+  if(strcmp(name,"difftheta") == 0){
+    if(!SCMatchRights(pCon,usUser)){
+      return 0;
+    }
+    self->allowedDeviation = atof(value);
+    SCparChange(pCon);
+    SCSendOK(pCon);
+    return 1;
+  } else if(strcmp(name,"maxindex") == 0){
+    if(!SCMatchRights(pCon,usUser)){
+      return 0;
+    }
+    self->indexSearchLimit = atoi(value);
+    SCparChange(pCon);
+    SCSendOK(pCon);
+    return 1;
+  } else if(strcmp(name,"maxlist") == 0){
+    if(!SCMatchRights(pCon,usUser)){
+      return 0;
+    }
+    self->maxSuggestions = atoi(value);
+    SCparChange(pCon);
+    SCSendOK(pCon);
+    return 1;
+  } else {
+    SCWrite(pCon,"ERROR: subcommand not recognized",eError);
+    return 0;
+  }
+}
+/*---------------------------------------------------------------------*/
+static int getUBCalcParameters(pUBCALC self, SConnection *pCon, char *name){
+  char pBueffel[255];
+  int ret = 0;
+
+  if(strcmp(name,"difftheta") == 0){
+    snprintf(pBueffel,255,"ubcalc.difftheta = %f",self->allowedDeviation);
+    ret = 1;
+  } else if(strcmp(name,"maxindex") == 0){
+    snprintf(pBueffel,255,"ubcalc.maxindex = %d", self->indexSearchLimit);
+    ret = 1;
+  } else if(strcmp(name,"maxlist") == 0){
+    snprintf(pBueffel,255,"ubcalc.maxindex = %d", self->maxSuggestions);
+    ret = 1;
+  } else {
+    snprintf(pBueffel,255,"ERROR: subcommand not known");
+  }
+  SCWrite(pCon,pBueffel,eValue);
+  return ret;
+}
+/*---------------------------------------------------------------------*/
+static void listPar(pUBCALC self, char *name, SConnection *pCon){
+  char pBueffel[255];
+
+  snprintf(pBueffel,255,"%s.difftheta = %f", name, self->allowedDeviation);
+  SCWrite(pCon,pBueffel,eValue);
+  snprintf(pBueffel,255,"%s.maxindex = %d", name, self->indexSearchLimit);
+  SCWrite(pCon,pBueffel,eValue);
+  snprintf(pBueffel,255,"%s.maxlist = %d", name, self->maxSuggestions);
+  SCWrite(pCon,pBueffel,eValue);
+
+}
+/*---------------------------------------------------------------------*/
+static int findIndex(pUBCALC self, SConnection *pCon, SicsInterp *pSics, 
+		       int argc, char *argv[]){
+  float two_theta, lambda;
+  pMotor pMot = NULL;
+  int status, numRef, i;
+  refIndex *index = NULL;
+  Tcl_DString list;
+  char pLine[255];
+
+  if(argc > 2){
+    two_theta = atof(argv[2]);
+  } else {
+    pMot = FindMotor(pSics,"stt");
+    if(pMot == NULL){
+      SCWrite(pCon,"ERROR: cannot find stt motor",eError);
+      return 0;
+    }
+    MotorGetSoftPosition(pMot,pCon,&two_theta);
+  }
+
+  GetLambda(self->hkl,&lambda);
+
+  numRef = self->maxSuggestions;
+  index = (refIndex *)malloc(numRef*sizeof(refIndex));
+  if(index == NULL){
+    SCWrite(pCon,"ERROR: out of memory allocating index list",eError);
+    return 0;
+  }
+  memset(index,0,numRef*sizeof(refIndex));
+
+  status = searchIndex(self->direct,lambda, two_theta,self->allowedDeviation,
+		       self->indexSearchLimit, index, numRef);
+
+  if(status < 0){
+    if(status == UBNOMEMORY){
+      SCWrite(pCon,"ERROR: out of memory searching indices",eError);
+      return 0;
+    } else if(status == REC_NO_VOLUME){
+      SCWrite(pCon,"ERROR: bad cell parameters",eError);
+      return 0;
+    } else {
+      SCWrite(pCon,"ERROR: unknown error code",eError);
+      return 0;
+    }
+  }
+
+  if(status < numRef) {
+    numRef = status;
+  }
+  Tcl_DStringInit(&list);
+  for(i = 0; i < numRef; i++){
+    snprintf(pLine,255," %3d %3d %3d %8.2f %8.2f %8.2f\r\n",
+	     (int)index[i].h, (int)index[i].k, (int)index[i].l,
+	     two_theta, index[i].t2calc,index[i].t2diff);
+    Tcl_DStringAppend(&list,pLine,-1);
+  }
+  if(numRef == 0){
+    Tcl_DStringAppend(&list,"No suitable reflections found",-1);
+  }
+  SCWrite(pCon,Tcl_DStringValue(&list),eValue);
+
+  Tcl_DStringFree(&list);
+  free(index);
+  return 1;
+}
+/*---------------------------------------------------------------------*/
+static int calcUB3Ref(pUBCALC self, SConnection *pCon){
+  float lambda;
+  MATRIX newUB = NULL;
+  int errCode = 1;
+  char pBueffel[256];
+
+  GetLambda(self->hkl,&lambda);
+  
+  newUB = calcUBFromThreeReflections(self->r1, self->r2, self->r3,
+				     lambda, &errCode);
+  if(newUB == NULL){
+    switch(errCode){
+    case UBNOMEMORY:
+      SCWrite(pCon,"ERROR: out of memory calculating UB",eError);
+      break;
+    case INVALID_LAMBDA:
+      SCWrite(pCon,"ERROR: bad value for wavelength",eError);
+      break;
+    case NOTRIGHTHANDED:
+      SCWrite(pCon,"ERROR: reflections are coplanar or do not form a right handed system",
+	      eError);
+      break;
+    default:
+      SCWrite(pCon,"ERROR: unknown error code from UB calculation",eError);
+      break;
+    }
+    return 0;
+  } else {
+    if(self->UB != NULL){
+      mat_free(self->UB);
+    }
+    self->UB = newUB;
+    SCSendOK(pCon);
+  }
+  return 1;
+}
+/*--------------------------------------------------------------------*/
+static int cellFromUBWrapper(pUBCALC self, SConnection *pCon){
+  int status;
+  char pBueffel[256];
+  lattice direct;
+
+  status = cellFromUB(self->UB,&direct);
+  if(status < 0){
+    switch(status){
+    case CELLNOMEMORY:
+      SCWrite(pCon,"ERROR: out of memory while calculating cell",eError);
+      break;
+    default:
+      SCWrite(pCon,"ERROR: unknown error calculating cell",eError);
+      break;
+    }
+    return 0;
+  } else {
+    snprintf(pBueffel,255," %8.4f %8.4f %8.4f %6.2f %6.2f %6.2f",
+	     direct.a, direct.b, direct.c,
+	     direct.alpha, direct.beta, direct.gamma);
+    SCWrite(pCon,pBueffel,eValue);
+  }
+  return 1;
+}
 /*----------------------------------------------------------------------*/
 int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData,		
 		  int argc, char *argv[]){
@@ -358,11 +566,17 @@ int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData,
     return readReflection(pCon,pSics, &self->r1,argc,argv);
   }  else if(strcmp(argv[1],"ref2") ==0){
     return readReflection(pCon,pSics, &self->r2,argc,argv);
+  } else if(strcmp(argv[1],"ref3") == 0){
+    return readReflection(pCon,pSics, &self->r3,argc,argv);
   } else if(strcmp(argv[1],"listub") == 0){
     listUB(pCon,self->UB);
     return 1;
   } else if(strcmp(argv[1],"calcub") == 0){
     return calcUB(self,pCon);
+  } else if(strcmp(argv[1],"calcub3ref") == 0){
+    return calcUB3Ref(self,pCon);
+  } else if(strcmp(argv[1],"cellfromub") == 0){
+    return cellFromUBWrapper(self,pCon);
   } else if(strcmp(argv[1],"listcell") == 0){
     listCell(pCon,argv[0],self->direct);
     return 1;
@@ -372,18 +586,27 @@ int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData,
   } else if(strcmp(argv[1],"listref2") == 0){
     listReflection(pCon,argv[0],"ref2",self->r2);
     return 1;
+  } else if(strcmp(argv[1],"listref3") == 0){
+    listReflection(pCon,argv[0],"ref3",self->r3);
+    return 1;
   } else if(strcmp(argv[1],"list") == 0){
     listCell(pCon,argv[0],self->direct);
     listReflection(pCon,argv[0],"ref1",self->r1);
     listReflection(pCon,argv[0],"ref2",self->r2);
+    listReflection(pCon,argv[0],"ref3",self->r3);
     listUB(pCon,self->UB);
+    listPar(self,argv[0],pCon);
     return 1;
-  } else if(strcmp(argv[1],"sendto") == 0){
+  } else if(strcmp(argv[1],"activate") == 0){
-    if(argc < 3){
+    return sendUBToHKL(pCon,pSics,self->hkl,self->UB);
-      SCWrite(pCon,"ERROR: Nothing to send UB too",eError);
+  } else if(strcmp(argv[1],"index") == 0){
-      return 0;
+    return findIndex(self,pCon,pSics,argc, argv);
+  } else {
+    if(argc > 2){
+      return setUBCalcParameters(self,pCon,argv[1],argv[2]);
+    } else {
+      return getUBCalcParameters(self,pCon,argv[1]);
     }
-    return sendUBToHKL(pCon,pSics,argv[2],self->UB);
   }
   return 1;
 }

ubfour.c

@@ -1,11 +1,12 @@
 /**
  * This is a library for calculating UB matrices for four circle diffraction.
- * The algorithm and setting definitions is from:
+ * The algorithm and settings definition is from:
  *
  *            Busing & Levy, Acta Cryst. (1967), 22, 457ff 
  *
  * Implemented:
  * - UB from cell cell constants and two reflections.
+ * - Brute force index search
  *
  * Mark Koennecke, March 2005
  */
@@ -13,6 +14,9 @@
 #include <assert.h>
 #include "vector.h"
 #include "trigd.h"
+#include "fourlib.h"
+#include "lld.h"
+#define ABS(x) (x < 0 ? -(x) : (x)) 
 /*--------------------------------------------------------------------------------------*/
 static MATRIX calcUVectorFromAngles(reflection r){
   MATRIX u;
@@ -50,6 +54,7 @@ MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1,
   MATRIX u1, u2, h1, h2;
   double ud[3];
   int status;
+  reflection r;
 
   *errCode = 1;
 
@@ -134,3 +139,234 @@ MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1,
 
   return UB;
 }
+/*-----------------------------------------------------------------------------------*/
+static void storeReflection(reflection r, double two_theta_obs, 
+			    double two_theta_calc, int list){
+  refIndex ri, test;
+  int count = 0, status, pos = 0;
+
+  ri.h = r.h;
+  ri.k = r.k;
+  ri.l = r.l;
+  ri.t2obs = two_theta_obs;
+  ri.t2calc = two_theta_calc;
+  ri.t2diff = ABS(two_theta_obs - two_theta_calc);
+
+  /*
+    locate the last entry bigger then us
+  */
+  status = LLDnodePtr2First(list);
+  while(status == 1){
+    LLDnodeDataTo(list,&test);
+    count++;
+    if(test.t2diff == ri.t2diff){
+      LLDnodeDataFrom(list,&ri);
+      return;
+    }
+    if(test.t2diff > ri.t2diff){
+      break;
+    }
+    status = LLDnodePtr2Next(list);
+  }
+  /*
+    special case: empty list
+  */
+  if(count == 0){
+    LLDnodeAppendFrom(list,&ri);
+    return;
+  }
+  /*
+    special case: append after last
+  */
+  LLDnodePtr2Last(list);
+  LLDnodeDataTo(list,&test);
+  if(ri.t2diff > test.t2diff){
+    LLDnodeAppendFrom(list,&ri);
+    return;
+  }
+
+  status = LLDnodePtr2First(list);
+  pos = 0;
+  while(status == 1){
+    LLDnodeDataTo(list,&test);
+    pos++;
+    if(pos == count){
+      LLDnodeInsertFrom(list,&ri);
+      return;
+    }
+    status = LLDnodePtr2Next(list);
+  }
+}
+/*----------------------------------------------------------------------------
+  u_transform(i) = u(i)*(2*sin(theta)/lambda)
+  -----------------------------------------------------------------------------*/ 
+static void uToScatteringVector(MATRIX u, double theta, double lambda){
+  double scale;
+  int i;
+
+  scale = (2. * Sind(theta))/lambda;
+  for(i = 0; i < 3; i++){
+    u[i][0] *= scale;
+  }
+}
+/*----------------------------------------------------------------------------*/
+static MATRIX buildHCHIMatrix(MATRIX u1, MATRIX u2, MATRIX u3){
+  int i;
+  MATRIX HCHI;
+
+  HCHI = mat_creat(3,3,ZERO_MATRIX);
+  if(HCHI == NULL){
+    return NULL;
+  }
+  for(i = 0; i < 3; i++){
+    HCHI[i][0] = u1[i][0];
+    HCHI[i][1] = u2[i][0];
+    HCHI[i][2] = u3[i][0];
+  }
+  return HCHI;
+}
+/*----------------------------------------------------------------------------*/
+static MATRIX buildIndexMatrix(reflection r1, reflection r2, reflection r3){
+  MATRIX HI;
+  int i;
+
+  HI = mat_creat(3,3,ZERO_MATRIX);
+  if(HI == NULL){
+    return NULL;
+  }
+  HI[0][0] = r1.h;
+  HI[1][0] = r1.k;
+  HI[2][0] = r1.l;
+
+  HI[0][1] = r2.h;
+  HI[1][1] = r2.k;
+  HI[2][1] = r2.l;
+
+  HI[0][2] = r3.h;
+  HI[1][2] = r3.k;
+  HI[2][2] = r3.l;
+
+  return HI;
+}
+/*-----------------------------------------------------------------------------*/
+MATRIX calcUBFromThreeReflections(reflection r1, reflection r2, reflection r3,
+				  double lambda, int *errCode){
+  MATRIX u1, u2, u3, HCHI, HI, HIINV,  UB;
+  double det;
+
+  if(lambda <= .1){
+    *errCode = INVALID_LAMBDA;
+    return NULL;
+  }
+
+  *errCode = 1;
+
+  u1 = calcUVectorFromAngles(r1);
+  u2 = calcUVectorFromAngles(r2);
+  u3 = calcUVectorFromAngles(r3);
+  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 = UBNOMEMORY;
+    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;
+}
+/*---------------------------------------------------------------------------------*/
+static int copyReflections(int list, refIndex index[], int maxIndex){
+  int count = 0, status;
+  refIndex ri;
+
+  status = LLDnodePtr2First(list);
+  while(status == 1 && count < maxIndex){
+    LLDnodeDataTo(list,&ri);
+    index[count] = ri;
+    status = LLDnodePtr2Next(list);
+    count++;
+  }
+  return count;
+}
+/*-----------------------------------------------------------------------------------
+ - matching reflections will be entered in to a list in a sorted way. This list
+   is copied into the index array.
+ - There is some waste here in allocating and deallocating the HC vector in the 
+   inner loop. I'am to lazy to resolve this.... May be I'am spared.....
+   -----------------------------------------------------------------------------------*/ 
+int searchIndex(lattice direct, double lambda, double two_theta, double max_deviation,
+		int limit, refIndex index[], int maxIndex){
+  int status, i, j, k, list;
+  MATRIX B, HC;
+  double theta, d;
+  reflection r;
+
+  B = mat_creat(3,3,UNIT_MATRIX);
+  if(B == NULL) {
+    return  UBNOMEMORY;
+  }
+
+  status = calculateBMatrix(direct,B);
+  if(status < 0) {
+    return status;
+  }
+
+  list = LLDcreate(sizeof(refIndex));
+  if(list <0) {
+    return  UBNOMEMORY;
+  }
+  
+  for(i = -limit; i < limit; i++){
+    r.h = (double)i;
+    for(j = -limit; j < limit; j++){
+      r.k = (double)j;
+      for(k = -limit; k < limit; k++){
+	r.l = (double)k;
+	HC = reflectionToHC(r,B);
+	status = calcTheta(lambda, HC, &d, &theta);
+	if(status == 1){
+	  if(ABS(two_theta - 2. * theta) <= max_deviation){
+	    storeReflection(r,two_theta, theta * 2., list);
+          }
+	}
+        killVector(HC);
+      }
+    }
+  }
+  mat_free(B);
+  status = copyReflections(list,index,maxIndex);
+  LLDdelete(list);
+  return status;
+} 

ubfour.h

@@ -3,8 +3,9 @@
  * The algorithm and setting definitions is from:
  *            Busing & Levy, Acta Cryst. (1967), 22, 457ff 
  *
- * This initial version only supports the calculation of the UB matrix from
+ * Implemented:
- * cell constants and two reflections.
+ * - UB from cell cell constants and two reflections.
+ * - Brute force index search
  *
  * Mark Koennecke, march 2005
  */
@@ -19,7 +20,8 @@
  * error codes: also see cell.h
  */
 #define UBNOMEMORY -200
-
+#define INVALID_LAMBDA -201
+#define NOTRIGHTHANDED -202
 /**
  * a reflection data structure holding the indices h,k,l and
  * the magic four circle angles two_theta, om, chi and phi.
@@ -37,4 +39,38 @@ typedef struct{
  * @return The resulting UB matrix or NULL on errors
  */
 MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1, reflection r2, int *errCode);
+/**
+ * calculate the UB matrix from three reflections. The three reflections must not be 
+ * coplanar and must reflect a right handed 
+ * @param r1 The first reflection
+ * @param r2 The second reflection
+ * @param r3 The third reflection.
+ * @param errCode a code indictang errors which happened.
+ * @return A UB matrix on success or NULL on errors. Then errcode will indicate
+ * the type of teh error.
+ */
+MATRIX calcUBFromThreeReflections(reflection r1, reflection r2, reflection r3,
+				  double lambda, int *errCode);
+/**
+ * a data structure holding an indexing suggestion
+ */
+typedef struct {
+  double h, k, l;                /* suggested index */
+  double t2obs, t2calc, t2diff;  /*  2 theta observed and calculated and the difference */
+}refIndex, *prefIndex;
+/**
+ * search for possible indexes matching the two theta value given. This is a brute force search
+ * @param direct The lattice constants of the crystal
+ * @param lambda The wavelength used
+ * @param two_theta The two theta value of the reflection 
+ * @param max_deviation maximum allowed diviation of calculated two thetas from tw-Theta given
+ * @param limit Index limit for the search
+ * @param index Preallocated array for storing index suggestions. The procedure will
+ * sort the entries in this array according to the difference in two theta
+ * @param maxIndex The number of entries allowed for index.
+ * @return The number of indexes in index or a negative error code when an error
+ * occurs.
+ */
+int searchIndex(lattice direct, double lambda, double two_theta, double max_deviation,
+		int limit, refIndex index[], int maxIndex);  
 #endif