- Initial commit of a UB calculation setup for four circle

diffractometers

cell.c

@@ -0,0 +1,150 @@
+/**
+ * this is a little library for performing crystallographic cell transformations
+ * for SICS. Some of the actual code was lifted from the Risoe program tascom.
+ *
+ * copyright: see file COPYRIGHT
+ *
+ * Mark Koennecke, March 2005
+ */
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include "trigd.h"
+#include "cell.h"
+
+/* define constants */
+#ifndef PI
+#define PI              (3.1415926536)  /* pi */
+#endif
+#define TWOPI           (2*PI)          /* 2*pi */
+
+/*******************************************************************************
+* Transform direct lattice to reciprocal lattice.
+*******************************************************************************/
+int directToReciprocalLattice(lattice direct, plattice reciprocal)
+{
+  double alfa, beta, gamma;
+  double cos_alfa, cos_beta, cos_gamma;
+  double sin_alfa, sin_beta, sin_gamma;
+  double ad, bd, cd;
+  double arg, vol;
+
+  alfa = direct.alpha;
+  beta = direct.beta;
+  gamma = direct.gamma;
+
+  cos_alfa = Cosd (alfa);
+  cos_beta = Cosd (beta);
+  cos_gamma = Cosd (gamma);
+
+  sin_alfa = Sind (alfa);
+  sin_beta = Sind (beta);
+  sin_gamma = Sind (gamma);
+
+  reciprocal->alpha = Acosd ((cos_beta*cos_gamma - cos_alfa)/sin_beta/sin_gamma);
+  reciprocal->beta =Acosd ((cos_alfa*cos_gamma - cos_beta)/sin_alfa/sin_gamma);
+  reciprocal->gamma = Acosd ((cos_alfa*cos_beta - cos_gamma)/sin_alfa/sin_beta);
+
+  ad = direct.a;
+  bd = direct.b;
+  cd = direct.c;
+
+  arg = 1 + 2*cos_alfa*cos_beta*cos_gamma - cos_alfa*cos_alfa -
+					    cos_beta*cos_beta -
+					    cos_gamma*cos_gamma;
+  if (arg < 0.0)
+    {
+      return REC_NO_VOLUME;
+    }
+
+  vol = ad*bd*cd*sqrt (arg);
+  reciprocal->a = bd*cd*sin_alfa/vol;
+  reciprocal->b = ad*cd*sin_beta/vol;
+  reciprocal->c = bd*ad*sin_gamma/vol;
+
+  return (0);
+}
+/*******************************************************************************
+* Transform reciprocal lattice to direct lattice.
+*******************************************************************************/
+int reciprocalToDirectLattice(lattice reciprocal, plattice direct)
+{
+  double alfa, beta, gamma;
+  double cos_alfa, cos_beta, cos_gamma;
+  double sin_alfa, sin_beta, sin_gamma;
+  double ar, br, cr;
+  double arg, vol;
+
+  alfa = reciprocal.alpha;
+  beta = reciprocal.beta;
+  gamma = reciprocal.gamma;
+
+  cos_alfa = Cosd (alfa);
+  cos_beta = Cosd (beta);
+  cos_gamma = Cosd (gamma);
+
+  sin_alfa = Sind (alfa);
+  sin_beta = Sind (beta);
+  sin_gamma = Sind (gamma);
+
+  direct->alpha = Acosd ((cos_beta*cos_gamma - cos_alfa)/sin_beta/sin_gamma);
+  direct->beta = Acosd ((cos_alfa*cos_gamma - cos_beta)/sin_alfa/sin_gamma);
+  direct->gamma =  Acosd ((cos_alfa*cos_beta - cos_gamma)/sin_alfa/sin_beta);
+
+  ar = reciprocal.a;
+  br = reciprocal.b;
+  cr = reciprocal.c;
+
+  arg = 1 + 2*cos_alfa*cos_beta*cos_gamma - cos_alfa*cos_alfa -
+					    cos_beta*cos_beta -
+					    cos_gamma*cos_gamma;
+  if (arg < 0.0)
+    {
+      return REC_NO_VOLUME;
+    }
+
+  vol = ar*br*cr*sqrt (arg);
+  direct->a = br*cr*sin_alfa/vol;
+  direct->b = ar*cr*sin_beta/vol;
+  direct->c = br*ar*sin_gamma/vol;
+
+  return (0);
+}
+/***************************************************************************************
+ * Build a B matrix
+ ***************************************************************************************/
+int calculateBMatrix(lattice direct, MATRIX B) {
+  lattice reciprocal;
+  int status;
+
+  assert(MatRow(B) == 3);
+  assert(MatCol(B) == 3);
+
+  status = directToReciprocalLattice(direct,&reciprocal);
+  if(status < 0) {
+    return status;
+  }
+  
+  mat_fill(B,ZERO_MATRIX);
+  
+  /*
+    top row
+  */
+  B[0][0] = reciprocal.a;
+  B[0][1] = reciprocal.b*Cosd(reciprocal.gamma);
+  B[0][2] = reciprocal.c*Cosd(reciprocal.beta);
+
+  /*
+    middle row
+  */
+  B[1][1] = reciprocal.b*Sind(reciprocal.gamma);
+  B[1][2] = -reciprocal.c*Sind(reciprocal.beta)*Cosd(direct.alpha);
+  
+  /*
+    bottom row
+  */
+  B[2][2] = 1./direct.c;
+  
+  return 1;
+}
+

cell.h

@@ -0,0 +1,46 @@
+/**
+ * this is a little library for performing crystallographic cell transformations
+ * for SICS. Some of the actual code was lifted from the Risoe program tascom.
+ *
+ * copyright: see file COPYRIGHT
+ *
+ * Mark Koennecke, March 2005
+ */
+#ifndef SICSCELL
+#define SICSCELL
+#include "matrix/matrix.h"
+
+/**
+ * error codes
+ */
+#define REC_NO_VOLUME -100
+
+/**
+ * lattice parameters: either reciprocal or direct
+ */
+  typedef struct {
+    double a,b,c;
+    double alpha, beta, gamma;
+  }lattice, *plattice;
+/**
+ * conversion from a direct lattice to the recipcrocal one.
+ * @param direct The input direct cell parameters.
+ * @param reciprocal The output reciprocal cell constants
+ * @return 0 on success, > 0 else
+ */ 
+  int directToReciprocalLattice(lattice direct, plattice reciprocal);
+/**
+ * conversion from a reciprocal  lattice to the directone.
+ * @param reciprocal  The input reciprocal  cell parameters.
+ * @param direct The output direct  cell constants
+ * @return 0 on success, > 0 else
+ */ 
+  int reciprocalToDirectLattice(lattice reciprocal, plattice direct);
+/**
+   * calculate a crystallographic B matrix from the cell constants
+   * @param direct The direct cell lattice to calculate B from
+   * @param B will be filled with the B matrix. MUST be 3x3
+   * @return 1 on success, an negative error code else
+   */
+int calculateBMatrix(lattice direct, MATRIX B);
+#endif

exeman.c

@@ -721,8 +721,9 @@ int ExeManagerWrapper(SConnection *pCon, SicsInterp *pSics, void *pData,
   self = (pExeMan)pData;
   assert(self != NULL);
 
-  strncpy(pBufferName,argv[1],255);
+
   if(argc > 1){
+    strncpy(pBufferName,argv[1],255);
     strtolower(argv[1]);
     status = handleBatchPath(self,pCon,argc,argv);
     if(status >= 0){

macro.c

@@ -958,8 +958,8 @@ static int ProtectedExec(ClientData clientData, Tcl_Interp *interp,
     Arg2Text(argc-1,&argv[1],pCommand,1023);
     strtolower(argv[0]);
     if(strcmp(argv[0],"fulltransact") == 0){
-      snprintf(pStart,1024,"TRANSACTIONSTART %s",pCommand);
+      snprintf(pStart,1023,"TRANSACTIONSTART %s",pCommand);
-      SCWrite(pCon,pStart,eError);
+      SCWrite(pCon,pStart,eError); 
     }
     iRet = InterpExecute(pSics,pCon,pCommand);
     SCWrite(pCon,"TRANSACTIONFINISHED",eError);

make_gen

@@ -15,13 +15,13 @@ SOBJ =  network.o ifile.o conman.o SCinter.o splitter.o passwd.o \
        script.o o2t.o  alias.o napi.o nxdata.o stringdict.o sdynar.o\
        histmem.o histdriv.o histsim.o  interface.o callback.o \
         event.o emon.o evcontroller.o evdriver.o simev.o perfmon.o \
-       danu.o  nxdict.o varlog.o stptok.o nread.o \
+       danu.o  nxdict.o varlog.o stptok.o nread.o trigd.o cell.o\
        scan.o fitcenter.o telnet.o token.o  wwildcard.o hklmot.o\
        tclev.o hkl.o integrate.o optimise.o dynstring.o nxutil.o \
        mesure.o uubuffer.o  commandlog.o udpquieck.o fourtable.o\
-       rmtrail.o  help.o nxupdate.o confvirtualmot.o \
+       rmtrail.o  help.o nxupdate.o confvirtualmot.o vector.o\
        simchop.o choco.o chadapter.o  trim.o scaldate.o \
-       hklscan.o xytable.o   exebuf.o exeman.o\
+       hklscan.o xytable.o   exebuf.o exeman.o ubfour.o ubcalc.o\
        circular.o  maximize.o sicscron.o  scanvar.o \
        d_sign.o d_mod.o tcldrivable.o stdscan.o diffscan.o\
        synchronize.o definealias.o  oscillate.o \

ofac.c

@@ -109,6 +109,7 @@
 #include "oscillate.h"
 #include "diffscan.h"
 #include "hklmot.h"
+#include "ubcalc.h"
 /*----------------------- Server options creation -------------------------*/
   static int IFServerOption(SConnection *pCon, SicsInterp *pSics, void *pData,
                      int argc, char *argv[])
@@ -294,6 +295,8 @@
 	       MakeDiffScan,NULL,NULL);
     AddCommand(pInter,"MakeHKLMot",
 	       HKLMotInstall,NULL,NULL);
+    AddCommand(pInter,"MakeUBCalc",
+	       MakeUBCalc,NULL,NULL);
 
 
     /*
@@ -356,6 +359,7 @@
     RemoveCommand(pSics,"MakeOscillator");
     RemoveCommand(pSics,"MakeDiffScan");
     RemoveCommand(pSics,"MakeHKLMot");
+    RemoveCommand(pSics,"MakeUBCalc");
 
     /*
       remove site specific installation commands
@@ -365,10 +369,7 @@
       site->RemoveSiteCommands(pSics);
     }
   }
-  
-   
 /*--------------------------------------------------------------------------*/
-
   int InitObjectCommands(pServer pServ, char *file)
   {
      SConnection *pCon = NULL;

scan.c

@@ -42,37 +42,6 @@
 #include "lld.h"
 #include "stdscan.h"
 
-/*-----------------------------------------------------------------------*/
-static char *fixExtension(char *filename)
-{
-  if(strstr(filename,".hdf") != NULL)
-  {
-    changeExtension(filename,"dat");
-  }
-  return filename;
-}
-/*------------------------------------------------------------------------*/
-  char *ScanMakeFileName(SicsInterp *pSics, SConnection *pCon)
-  {
-      pSicsVariable pPath = NULL, pPref = NULL, pEnd = NULL;
-      char *pRes = NULL;
-      int iLen, iNum, iYear;
-      char pNumText[10];
-      CommandList *pCom = NULL;
-      
-      /*
-       make a simulated filename if in simulation mode
-      */
-      if(pServ->simMode)
-	return strdup("sim001001901.sim");
-
-      pRes = makeFilename(pSics,pCon);
-      if(pRes == NULL) 
-      {
-	pRes = strdup("emergency.scn");
-      }
-      return fixExtension(pRes);
-  }
 /*---------------------------------------------------------------------------*/
   static void DeleteCountEntry(void *pData)
   {
@@ -688,21 +657,7 @@ CountEntry CollectCounterData(pScanData self)
          return 0;
        }
        
-       /* allocate a new data file */
+
-       pPtr = ScanMakeFileName(self->pSics,self->pCon);
-       if(!pPtr)
-       {
-          SCWrite(self->pCon,
-                  "ERROR: cannot allocate new data filename, Scan aborted",
-                  eError);
-          self->pCon = NULL;
-          self->pSics = NULL;
-          return 0;
-       }
-       sprintf(pBueffel,"Writing data file: %s ...",pPtr);
-       SCWrite(self->pCon,pBueffel,eWarning);
-       strcpy(self->pFile,pPtr);
-       free(pPtr);
        iRet = self->WriteHeader(self);
        if(!iRet)
        {

stdscan.c

@@ -30,6 +30,37 @@
 #include "site.h"
 #include "lld.h"
 
+/*-----------------------------------------------------------------------*/
+static char *fixExtension(char *filename)
+{
+  if(strstr(filename,".hdf") != NULL)
+  {
+    changeExtension(filename,"dat");
+  }
+  return filename;
+}
+/*------------------------------------------------------------------------*/
+  char *ScanMakeFileName(SicsInterp *pSics, SConnection *pCon)
+  {
+      pSicsVariable pPath = NULL, pPref = NULL, pEnd = NULL;
+      char *pRes = NULL;
+      int iLen, iNum, iYear;
+      char pNumText[10];
+      CommandList *pCom = NULL;
+      
+      /*
+       make a simulated filename if in simulation mode
+      */
+      if(pServ->simMode)
+	return strdup("sim001001901.sim");
+
+      pRes = makeFilename(pSics,pCon);
+      if(pRes == NULL) 
+      {
+	pRes = strdup("emergency.scn");
+      }
+      return fixExtension(pRes);
+  }
 /*---------------------------------------------------------------------------*/
   int WriteHeader(pScanData self)
   {
@@ -343,6 +374,7 @@
      float fVal;
      char pBueffel[512];
      char pMessage[1024];
+     char *pPtr = NULL;
      
      assert(self);
      assert(self->iNP > 0);
@@ -390,6 +422,22 @@
      SetCounterMode((pCounter)self->pCounterData,self->iMode);
      SetCounterPreset((pCounter)self->pCounterData, self->fPreset);
      self->iCounts = 0;
+
+     /* allocate a new data file */
+     pPtr = ScanMakeFileName(self->pSics,self->pCon);
+     if(!pPtr)
+     {
+	 SCWrite(self->pCon,
+                  "ERROR: cannot allocate new data filename, Scan aborted",
+                  eError);
+	 self->pCon = NULL;
+	 self->pSics = NULL;
+	 return 0;
+     }
+     snprintf(pBueffel,511,"Writing data file: %s ...",pPtr);
+     SCWrite(self->pCon,pBueffel,eWarning);
+     strcpy(self->pFile,pPtr);
+     free(pPtr);
      
      return 1;
   }  
@@ -402,6 +450,7 @@
      float fVal;
      char pBueffel[512];
      char pMessage[1024];
+     char *pPtr  = NULL;
      
      assert(self);
      assert(self->iNP > 0);
@@ -431,6 +480,22 @@
      SetCounterMode((pCounter)self->pCounterData,self->iMode);
      SetCounterPreset((pCounter)self->pCounterData, self->fPreset);
      self->iCounts = 0;
+
+     /* allocate a new data file */
+     pPtr = ScanMakeFileName(self->pSics,self->pCon);
+     if(!pPtr)
+     {
+	 SCWrite(self->pCon,
+                  "ERROR: cannot allocate new data filename, Scan aborted",
+                  eError);
+	 self->pCon = NULL;
+	 self->pSics = NULL;
+	 return 0;
+     }
+     snprintf(pBueffel,511,"Writing data file: %s ...",pPtr);
+     SCWrite(self->pCon,pBueffel,eWarning);
+     strcpy(self->pFile,pPtr);
+     free(pPtr);
      
      return 1;
   }  

stdscan.h

@@ -11,6 +11,10 @@
 #ifndef SICSSTDSCAN
 #define SICSSTDSCAN
 
+  /**
+   * make a scan file name
+   */
+  char *ScanMakeFileName(SicsInterp *pSics, SConnection *pCon);
   /**
     * write the header of the scan file
     */

trigd.c

@@ -0,0 +1,92 @@
+/**
+ * This is a library of trigonmetric functions acting upon proper
+ * angles and not radians. Lifted from the Risoe tascom code
+ *
+ * March 2005
+ */
+#include <math.h>
+#include "trigd.h"
+
+/* define constants */
+#ifndef PI
+#define PI              (3.1415926536)  /* pi */
+#endif
+#define DEGREE_RAD      (PI/180.0)      /* Radians per degree */
+/*******************************************************************************/
+double Sign(double d)
+{
+  if(d < .0){
+    return -1;
+  } else {
+    return 1;
+  }
+}
+/*******************************************************************************
+* Sinus of angle in degrees.
+*******************************************************************************/
+extern double Sind (double x)
+{
+  return (sin (x*DEGREE_RAD));
+}
+/*******************************************************************************
+* Tangens of angle in degrees.
+*******************************************************************************/
+extern double Tand(double x)
+{
+  return (tan(x*DEGREE_RAD));
+}
+/*******************************************************************************
+* Cosinus of angle in degrees.
+*******************************************************************************/
+extern double Cosd (double x)
+{
+  return (cos (x*DEGREE_RAD));
+}
+/*******************************************************************************
+* Atan of angle in degrees.
+*******************************************************************************/
+extern double Atand (double x)
+{
+  double data;
+
+  data = (atan(x)/DEGREE_RAD);
+  return (data);
+}
+/*******************************************************************************
+* Atan2 of angle in degrees.
+*******************************************************************************/
+extern double Atand2 (double x)
+{
+  double data;
+
+  data = (atan(x)/DEGREE_RAD);
+  return (data);
+}
+/*******************************************************************************
+* Acos of angle in degrees.
+*******************************************************************************/
+extern double Acosd (double x)
+{
+  double data;
+
+  data = acos(x)/DEGREE_RAD;
+  return (data);
+}
+/*******************************************************************************
+* Asin of angle in degrees.
+*******************************************************************************/
+extern double Asind (double x)
+{
+  double data;
+
+  data = x*x;
+  if (data == 1.0)
+    return (180.00 - Sign (x)*90.00);
+  else if (data > 1)
+  {
+    return (0);
+  }
+  else
+    return (Atand (x/sqrt (1 - data)));
+
+}

trigd.h

@@ -0,0 +1,17 @@
+/**
+ * This is a library of trigonmetric functions acting upon proper
+ * angles and not radians. Lifted from the Risoe tascom code
+ *
+ * March 2005
+ */
+#ifndef SICSTRIGD
+#define SICSTRIGD
+  double Sign(double d);
+  double Sind (double x);
+  double Tand(double x);
+  double Cosd (double x);
+  double Atand (double x);
+  double Atand2 (double x);
+  double Acosd (double x);
+  double Asind (double x);
+#endif

ubcalc.c

@@ -0,0 +1,389 @@
+/*----------------------------------------------------------------------
+  UB calculation routines for four circle diffraction.
+  This is the interpreter interface to functionality implemented
+  in fourlib.c
+
+  copyright: see file COPYRIGHT
+
+  Mark Koennecke, March 2005
+-----------------------------------------------------------------------*/
+#include <stdlib.h>
+#include <assert.h>
+#include <errno.h>
+#include "tcl.h"
+#include "fortify.h"
+#include "sics.h"
+#include "ubfour.h"
+#include "ubcalc.h"
+#include "motor.h"
+#include "hkl.h"
+/*---------------------------------------------------------------------*/
+typedef struct {
+  pObjectDescriptor pDes;
+  lattice direct;
+  reflection r1, r2;
+  MATRIX UB;
+} UBCALC, *pUBCALC;
+/*----------------------------------------------------------------------*/
+static void killUBCALC(void *pData){
+  pUBCALC self = (pUBCALC)pData;
+  if(self == NULL){
+    return;
+  }
+  if(self->pDes != NULL){
+    DeleteDescriptor(self->pDes);
+  }
+  if(self->UB != NULL){
+    mat_free(self->UB);
+  }
+  free(self);
+}
+/*--------------------------------------------------------------------*/
+static int SaveUBCalc(void *data, char *name, FILE *fd){
+  pUBCALC self = (pUBCALC)data;
+  if(self == NULL){
+    return 0;
+  }
+  fprintf(fd,"%s cell %f %f %f %f %f %f\n", name, self->direct.a,
+	  self->direct.b, self->direct.c, self->direct.alpha,
+	  self->direct.beta, self->direct.gamma);
+  fprintf(fd,"%s ref1 %f %f %f  %f %f %f %f\n",name,
+	  self->r1.h, self->r1.k, self->r1.l,
+	  self->r1.s2t, self->r1.om, self->r1.chi, self->r1.phi);
+  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);
+  return 1;
+}
+/*---------------------------------------------------------------------*/
+static pUBCALC makeUBCALC(){
+  pUBCALC pNew = NULL;
+
+  pNew = (pUBCALC)malloc(sizeof(UBCALC));
+  if(pNew == NULL){
+    return NULL;
+  }
+  memset(pNew,0,sizeof(UBCALC));
+  pNew->pDes = CreateDescriptor("UBcalc");
+  pNew->UB = mat_creat(3,3,UNIT_MATRIX);
+  if(pNew->pDes == NULL || pNew->UB == NULL){
+    return NULL;
+  }
+  pNew->pDes->SaveStatus = SaveUBCalc;
+  return pNew;
+}
+/*----------------------------------------------------------------------*/
+int MakeUBCalc(SConnection *pCon, SicsInterp *pSics, void *pData,		
+	       int argc, char *argv[]){
+  pUBCALC pNew = NULL;
+  int status;
+
+  pNew = makeUBCALC();
+  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);
+  } else {
+    status = AddCommand(pSics,"ubcalc",UBCalcWrapper,killUBCALC,pNew);
+  }
+  if(status != 1){
+    SCWrite(pCon,"ERROR: failed to create duplicate UBCALC module",eError);
+  }
+  return status;
+}
+/*---------------------------------------------------------------------*/
+static int readCell(SConnection *pCon, pUBCALC self, int argc, char *argv[]){
+  int status;
+  Tcl_Interp *pTcl = InterpGetTcl(pServ->pSics);
+  char pBueffel[256];
+
+  if(argc < 8){
+    SCWrite(pCon,"ERROR: insufficient number of arguments to ubcalc cell",
+	    eError);
+    return 0;
+  }
+  status = Tcl_GetDouble(pTcl,argv[2],&self->direct.a);
+  if(status != TCL_OK){
+    snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[2]);
+    SCWrite(pCon,pBueffel,eError);
+    return 0;
+  }
+  status = Tcl_GetDouble(pTcl,argv[3],&self->direct.b);
+  if(status != TCL_OK){
+    snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[3]);
+    SCWrite(pCon,pBueffel,eError);
+    return 0;
+  }
+  status = Tcl_GetDouble(pTcl,argv[4],&self->direct.c);
+  if(status != TCL_OK){
+    snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[4]);
+    SCWrite(pCon,pBueffel,eError);
+    return 0;
+  }
+  status = Tcl_GetDouble(pTcl,argv[5],&self->direct.alpha);
+  if(status != TCL_OK){
+    snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[5]);
+    SCWrite(pCon,pBueffel,eError);
+    return 0;
+  }
+  status = Tcl_GetDouble(pTcl,argv[6],&self->direct.beta);
+  if(status != TCL_OK){
+    snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[6]);
+    SCWrite(pCon,pBueffel,eError);
+    return 0;
+  }
+  status = Tcl_GetDouble(pTcl,argv[7],&self->direct.gamma);
+  if(status != TCL_OK){
+    snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[7]);
+    SCWrite(pCon,pBueffel,eError);
+    return 0;
+  }
+  SCparChange(pCon);
+  SCSendOK(pCon);
+  return 1;
+}
+/*--------------------------------------------------------------------*/
+static int readRefMotors(SConnection *pCon, SicsInterp *pSics, 
+			 reflection *r){
+  pMotor pMot = NULL;
+  float val;
+
+  pMot = FindMotor(pSics,"stt");
+  if(pMot == NULL){
+    SCWrite(pCon,"ERROR: cannot find stt motor",eError);
+    return 0;
+  }
+  MotorGetSoftPosition(pMot,pCon,&val);
+  r->s2t = val;
+
+  pMot = FindMotor(pSics,"om");
+  if(pMot == NULL){
+    SCWrite(pCon,"ERROR: cannot find om motor",eError);
+    return 0;
+  }
+  MotorGetSoftPosition(pMot,pCon,&val);
+  r->om = val;
+
+  pMot = FindMotor(pSics,"chi");
+  if(pMot == NULL){
+    SCWrite(pCon,"ERROR: cannot find chi motor",eError);
+    return 0;
+  }
+  MotorGetSoftPosition(pMot,pCon,&val);
+  r->chi = val;
+
+  pMot = FindMotor(pSics,"phi");
+  if(pMot == NULL){
+    SCWrite(pCon,"ERROR: cannot find phi motor",eError);
+    return 0;
+  }
+  MotorGetSoftPosition(pMot,pCon,&val);
+  r->phi = val;
+
+  SCparChange(pCon);
+  SCSendOK(pCon);
+  return 1;
+}
+/*---------------------------------------------------------------------*/
+static int readReflection(SConnection *pCon, SicsInterp *pSics, reflection *r,
+			  int argc, char *argv[]){
+  char pBueffel[255];
+  int status;
+  Tcl_Interp *pTcl = InterpGetTcl(pSics);
+
+  if(argc < 5){
+    SCWrite(pCon,"ERROR: not enough arguments to ubcalc ref1,2",eError);
+    return 0;
+  }
+  status = Tcl_GetDouble(pTcl,argv[2],&r->h);
+  if(status != TCL_OK){
+    snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[2]);
+    SCWrite(pCon,pBueffel,eError);
+    return 0;
+  }
+  status = Tcl_GetDouble(pTcl,argv[3],&r->k);
+  if(status != TCL_OK){
+    snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[3]);
+    SCWrite(pCon,pBueffel,eError);
+    return 0;
+  }
+  status = Tcl_GetDouble(pTcl,argv[4],&r->l);
+  if(status != TCL_OK){
+    snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[4]);
+    SCWrite(pCon,pBueffel,eError);
+    return 0;
+  }
+
+  if(argc >= 9){
+    status = Tcl_GetDouble(pTcl,argv[5],&r->s2t);
+    if(status != TCL_OK){
+      snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[5]);
+      SCWrite(pCon,pBueffel,eError);
+      return 0;
+    }
+    status = Tcl_GetDouble(pTcl,argv[6],&r->om);
+    if(status != TCL_OK){
+      snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[6]);
+      SCWrite(pCon,pBueffel,eError);
+      return 0;
+    }
+    status = Tcl_GetDouble(pTcl,argv[7],&r->chi);
+    if(status != TCL_OK){
+      snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[7]);
+      SCWrite(pCon,pBueffel,eError);
+      return 0;
+    }
+    status = Tcl_GetDouble(pTcl,argv[8],&r->phi);
+    if(status != TCL_OK){
+      snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[8]);
+      SCWrite(pCon,pBueffel,eError);
+      return 0;
+    }
+    SCparChange(pCon);
+    SCSendOK(pCon);
+    return 1;
+  } else {
+    return readRefMotors(pCon,pSics,r);
+  }
+} 
+/*---------------------------------------------------------------------*/
+static void listUB(SConnection *pCon, MATRIX UB){
+  Tcl_DString list;
+  char pBueffel[255];
+  int i;
+
+
+  Tcl_DStringInit(&list);
+  if(UB == NULL){
+    Tcl_DStringAppend(&list,"NO UB",-1);
+  } else {
+    Tcl_DStringAppend(&list,"UB = ", -1);
+    for(i = 0; i < 3; i++){
+      snprintf(pBueffel,255,"%f %f %f\n", UB[i][0],
+	       UB[i][1],UB[i][2]);
+      Tcl_DStringAppend(&list,pBueffel,-1);
+    }
+  }
+  SCWrite(pCon,Tcl_DStringValue(&list),eValue);
+  Tcl_DStringFree(&list);
+}
+/*---------------------------------------------------------------------*/
+static int calcUB(pUBCALC self, SConnection *pCon){
+  MATRIX newUB = NULL;
+  int err = 1;
+
+  newUB = calcUBFromCellAndReflections(self->direct, self->r1, self->r2, &err);
+  if(newUB == NULL){
+    switch(err){
+    case UBNOMEMORY:
+      SCWrite(pCon,"ERROR: out of memory while calculating UB",eError);
+      return 0;
+      break;
+    case REC_NO_VOLUME:
+      SCWrite(pCon,"ERROR: bad cell constants",eError);
+      return 0;
+      break;
+    default:
+      SCWrite(pCon,"ERROR: unknown error on UB matrix calculation",eError);
+      return 0;
+    }
+  } else {
+    if(self->UB != NULL){
+      mat_free(self->UB);
+    }
+    self->UB = newUB;
+    SCSendOK(pCon);
+    return 1;
+  }
+}
+/*---------------------------------------------------------------------*/
+static void listCell(SConnection *pCon, char *name, lattice direct){
+  char pBueffel[255];
+
+  snprintf(pBueffel,255,"%s.cell = %f %f %f %f %f %f",
+	   name,direct.a, direct.b,direct.c,
+	   direct.alpha,direct.beta,direct.gamma);
+  SCWrite(pCon,pBueffel,eValue);
+  
+}
+/*---------------------------------------------------------------------*/
+static void listReflection(SConnection *pCon, char *name, 
+			   char *refName, reflection r){
+  char pBueffel[255];
+
+  snprintf(pBueffel,255,"%s.%s = %f %f %f   %f %f %f %f",
+	   name,refName, r.h, r.k, r.l, r.s2t, r.om, r.chi, r.phi);
+  SCWrite(pCon,pBueffel,eValue);
+}
+/*---------------------------------------------------------------------*/
+static int sendUBToHKL(SConnection *pCon, SicsInterp *pSics,
+			char *name, MATRIX UB){
+  pHKL hkl = NULL;
+  char pBueffel[256];
+  float ub[9];
+  int i;
+
+  hkl = FindCommandData(pSics,name,"4-Circle-Calculus");
+  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];
+    ub[i+6] = UB[2][i];
+  }
+  SetUB(hkl,ub);
+  SCSendOK(pCon);
+  return 1;
+} 
+/*----------------------------------------------------------------------*/
+int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData,		
+		  int argc, char *argv[]){
+  pUBCALC self = (pUBCALC)pData;
+  
+  assert(self);
+  if(argc < 2){
+    SCWrite(pCon,"Insuffcient number of arguments to ubcalc",eError);
+    return 0;
+  }
+  strtolower(argv[1]);
+
+  if(strcmp(argv[1],"cell") == 0){
+    return readCell(pCon, self, argc, argv);
+  } else if(strcmp(argv[1],"ref1") == 0){
+    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],"listub") == 0){
+    listUB(pCon,self->UB);
+    return 1;
+  } else if(strcmp(argv[1],"calcub") == 0){
+    return calcUB(self,pCon);
+  } else if(strcmp(argv[1],"listcell") == 0){
+    listCell(pCon,argv[0],self->direct);
+    return 1;
+  } else if(strcmp(argv[1],"listref1") == 0){
+    listReflection(pCon,argv[0],"ref1",self->r1);
+    return 1;
+  } else if(strcmp(argv[1],"listref2") == 0){
+    listReflection(pCon,argv[0],"ref2",self->r2);
+    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);
+    listUB(pCon,self->UB);
+    return 1;
+  } else if(strcmp(argv[1],"sendto") == 0){
+    if(argc < 3){
+      SCWrite(pCon,"ERROR: Nothing to send UB too",eError);
+      return 0;
+    }
+    return sendUBToHKL(pCon,pSics,argv[2],self->UB);
+  }
+  return 1;
+}

ubcalc.h

@@ -0,0 +1,19 @@
+
+/*----------------------------------------------------------------------
+  UB caclualtion routines for four circle diffraction.
+  This is the interpreter interface to functionality implemented
+  in fourlib.c
+
+  copyright: see file COPYRIGHT
+
+  Mark Koennecke, March 2005
+-----------------------------------------------------------------------*/
+#ifndef SICSUBCALC
+#define SICSUBCALC
+
+int MakeUBCalc(SConnection *pCon,SicsInterp *pSics, void *pData,                
+                int argc, char *argv[]);
+int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData,            
+                int argc, char *argv[]);
+
+#endif   

ubcalc.w

@@ -0,0 +1,28 @@
+\subsection{UB Matrix Calculation for Four Circle Diffractometers}
+This module helps in the calculation of UB matrices for four
+ circle diffraction. This is only an interpreter interface, the
+ actual functionality is in the ubfour.h, .c files. These are documented
+ in ubfour.h in order to be able to give this away separatly as a
+ library.
+
+@o ubcalc.h @{
+/*----------------------------------------------------------------------
+  UB caclualtion routines for four circle diffraction.
+  This is the interpreter interface to functionality implemented
+  in fourlib.c
+
+  copyright: see file COPYRIGHT
+
+  Mark Koennecke, March 2005
+-----------------------------------------------------------------------*/
+#ifndef SICSUBCALC
+#define SICSUBCALC
+
+int MakeUBCalc(SConnection *pCon,SicsInterp *pSics, void *pData,		
+		int argc, char *argv[]);
+int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData,		
+		int argc, char *argv[]);
+
+#endif   
+@} 
+ 

ubfour.c

@@ -0,0 +1,136 @@
+/**
+ * This is a library for calculating UB matrices for four circle diffraction.
+ * The algorithm and setting definitions is from:
+ *
+ *            Busing & Levy, Acta Cryst. (1967), 22, 457ff 
+ *
+ * Implemented:
+ * - UB from cell cell constants and two reflections.
+ *
+ * Mark Koennecke, March 2005
+ */
+#include "ubfour.h"
+#include <assert.h>
+#include "vector.h"
+#include "trigd.h"
+/*--------------------------------------------------------------------------------------*/
+static MATRIX calcUVectorFromAngles(reflection r){
+  MATRIX u;
+
+  u = makeVector();
+  if(u == NULL){
+    return NULL;
+  }
+  vectorSet(u,0,Cosd(r.om)*Cosd(r.chi)*Cosd(r.phi) - Sind(r.om)*Sind(r.phi));
+  vectorSet(u,1,Cosd(r.om)*Cosd(r.chi)*Sind(r.phi) + Sind(r.om)*Cosd(r.phi));
+  vectorSet(u,2,Cosd(r.om)*Sind(r.chi));
+  
+  return u;	    
+}
+/*--------------------------------------------------------------------------------------*/
+static MATRIX reflectionToHC(reflection r, MATRIX B){
+  MATRIX h = NULL, hc = NULL;
+
+  h = makeVector();
+  if(h == NULL){
+    return NULL;
+  }
+  vectorSet(h,0,r.h);
+  vectorSet(h,1,r.k);
+  vectorSet(h,2,r.l);
+
+  hc = mat_mul(B,h);
+  killVector(h);
+  return hc;
+}
+/*---------------------------------------------------------------------------------------*/
+MATRIX calcUBFromCellAndReflections(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;
+
+  *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 = calcUVectorFromAngles(r1);
+  u2 = calcUVectorFromAngles(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;
+}

ubfour.h

@@ -0,0 +1,40 @@
+/**
+ * This is a library for calculating UB matrices for four circle diffraction.
+ * 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
+ * cell constants and two reflections.
+ *
+ * Mark Koennecke, march 2005
+ */
+#ifndef SICSUBFOUR
+#define SICSUBFOUR
+
+#include <stdio.h>
+#include "matrix/matrix.h"
+#include "cell.h"
+
+/**
+ * error codes: also see cell.h
+ */
+#define UBNOMEMORY -200
+
+/**
+ * a reflection data structure holding the indices h,k,l and
+ * the magic four circle angles two_theta, om, chi and phi.
+ */
+typedef struct{
+  double h,k,l;
+  double s2t,om,chi,phi;
+}reflection;
+/**
+ * calculate a UB matrix from cell constants and two reflections
+ * @param direct The direct cell constants
+ * @param r1 The first reflection.
+ * @param r2 The second reflection.
+ * @param errCode an error indicator if things go wrong.
+ * @return The resulting UB matrix or NULL on errors
+ */
+MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1, reflection r2, int *errCode);
+#endif 

vector.c

@@ -0,0 +1,145 @@
+/**
+ * This is a little collection of vector functions for use with the UB
+ * calculation routines. Therefore vectors are mapped onto the matrix
+ * package. Strictly speaking all routines only work properly on 3D
+ * vectors.
+ *
+ * copyright: see file COPYRIGHT
+ *
+ * Mark Koennecke, March 2005
+ */
+#include <stdlib.h>
+#include <assert.h>
+#include <math.h>
+#include "vector.h"
+/*----------------------------------------------------------------------*/
+MATRIX makeVector(){
+  return mat_creat(3,1,ZERO_MATRIX);
+}
+/*---------------------------------------------------------------------*/
+MATRIX makeVectorInit(double val[3]){
+  int i;
+  MATRIX result = NULL;
+
+  result = makeVector();
+  if(result == NULL){
+    return result;
+  }
+  for(i = 0; i < 3; i++){
+    result[i][0] = val[i];
+  } 
+  return result;
+}
+/*---------------------------------------------------------------------*/
+void vectorToArray(MATRIX v, double val[3]){
+  int i;
+
+  assert(MatCol(v) == 3);
+  for(i = 0; i < 3; i++){
+    val[i] = v[i][0];
+  }
+}
+/*----------------------------------------------------------------------*/
+void killVector(MATRIX v){
+  mat_free(v);
+}
+/*----------------------------------------------------------------------*/
+void vectorSet(MATRIX v, int idx, double value){
+  assert(idx >= 0 && idx < 3);
+
+  v[idx][0] = value;
+}
+/*----------------------------------------------------------------------*/
+double vectorGet(MATRIX v, int idx){
+  assert(idx >= 0 && idx < 3);
+
+  return v[idx][0];
+}
+/*----------------------------------------------------------------------*/
+double vectorLength(MATRIX v){
+  assert(MatRow(v) == 3);
+
+  return sqrt(v[0][0]*v[0][0] + v[1][0]*v[1][0] + v[2][0]*v[2][0]);
+}
+/*---------------------------------------------------------------------*/
+void normalizeVector(MATRIX v){
+  int i;
+  double norm;
+
+  norm = vectorLength(v);
+  if(norm > .001) {
+    for(i = 0; i < 3; i++){
+      v[i][0] /= norm;
+    }
+  } else {
+    for(i = 0; i < 3; i++){
+      v[i][0] = .0;
+    }
+  }
+}
+/*----------------------------------------------------------------------*/
+double vectorDotProduct(MATRIX v1, MATRIX v2){
+  double sum;
+  int i;
+
+  assert(MatRow(v1) == MatRow(v2));
+
+  sum = .0;
+  for(i = 0; i < MatRow(v1); i++){
+    sum += v1[i][0]*v2[i][0];
+  }
+  return sum;
+}
+/*----------------------------------------------------------------------*/
+MATRIX vectorCrossProduct(MATRIX v1, MATRIX v2){
+  MATRIX result;
+
+  assert(MatRow(v1) == 3);
+  assert(MatRow(v2) == 3);
+
+  result = makeVector();
+  if(result == NULL){
+    return NULL;
+  }
+  vectorSet(result,0,vectorGet(v1,1)*vectorGet(v2,2) - vectorGet(v1,2)*vectorGet(v2,1));
+  vectorSet(result,1,vectorGet(v1,2)*vectorGet(v2,0) - vectorGet(v1,0)*vectorGet(v2,2));
+  vectorSet(result,2,vectorGet(v1,0)*vectorGet(v2,1) - vectorGet(v1,1)*vectorGet(v2,0));
+  return result;
+}
+/*-------------------------------------------------------------------------*/
+MATRIX matFromTwoVectors(MATRIX v1, MATRIX v2){
+  MATRIX a1, a2, a3, result;
+  int i;
+
+  a1 = mat_copy(v1);
+  if(a1 == NULL){
+    return NULL;
+  }
+  normalizeVector(a1);
+
+  a3 = vectorCrossProduct(v1,v2);
+  if(a3 == NULL){
+    return NULL;
+  }
+  normalizeVector(a3);
+
+  a2 = vectorCrossProduct(a1,a3);
+  if(a2 == NULL){
+    return NULL;
+  }
+
+  result = mat_creat(3,3,ZERO_MATRIX);
+  if(result == NULL){
+    return NULL;
+  }
+
+  for(i = 0; i < 3; i++){
+    result[i][0] = vectorGet(a1,i);
+    result[i][1] = vectorGet(a2,i);
+    result[i][2] = vectorGet(a3,i);
+  }
+  killVector(a1);
+  killVector(a2);
+  killVector(a3);
+  return result;
+}

vector.h

@@ -0,0 +1,93 @@
+/**
+ * This is a little collection of vector functions for use with the UB
+ * calculation routines. Therefore vectors are mapped onto the matrix
+ * package. Strictly speaking all routine only work properly on 3D
+ * vectors.
+ *
+ * copyright: see file COPYRIGHT
+ *
+ * Mark Koennecke, March 2005
+ */
+#ifndef SICSVECTOR
+#define SICSVECTOR
+#include <stdio.h>
+#include "matrix/matrix.h"
+
+/**
+ * make a 3D vector
+ * @return a 3D vector expressed as a MATRIX. The caller is resonsible
+ * for removing any memory associated with this vector.
+ */
+MATRIX makeVector();
+/**
+ * make a vector and initialize with values given in val
+ * @param val Array of values for the new vector
+ * @return a new vector with the values in val
+ */
+MATRIX makeVectorInit(double val[3]);
+/**
+ * copy the value in v to the array val
+ * @param v The vector to copy
+ * @param val The array to which to copy the values
+ */
+void vectorToArray(MATRIX v, double val[3]);
+/**
+ * delete the vector v
+ * @param v The vector to delete.
+ */
+void killVector(MATRIX v);
+/**
+ * set a component of a vector
+ * @param v The vector to which to apply the new value
+ * @param idx The index of the vector component to set
+ * @param val The new value for the component.
+ */
+void vectorSet(MATRIX v, int idx, double value);
+/**
+ * get the value of a vector component
+ * @param v The vector to query
+ * @param idx The index of the component
+ * @return The value of the vectors component.
+ */
+double vectorGet(MATRIX v, int idx);
+/**
+ * calculate the length of the vector v
+ * @param v The vector to calculate the length for
+ * @return The length of the vector
+ */
+double vectorLength(MATRIX v);
+/**
+ * normalize the vector by making it unit length
+ * @param v The vector to normalize.
+ */
+void normalizeVector(MATRIX v);
+/**
+ * calculate the dot or scalar product of two vectors
+ * @param v1 First vector
+ * @param v2 Second vector
+ * @return The dot product v1*v2
+ */
+double vectorDotProduct(MATRIX v1, MATRIX v2);
+/**
+ * calculate the cross product of vectors v1 and v2. This is only
+ * correct when both vectors are expressed in terms of a 
+ * cartesian coordinate system.
+ * @param v1 The first vector
+ * @param v2 The second vector
+ * @return The cross product of the vectors v1, v2
+ */
+MATRIX vectorCrossProduct(MATRIX v1, MATRIX v2);
+/**
+ * this is a special function used in UB matrix calculations. 
+ * It first calculates a coordinate system from the two vectors
+ * where:
+ * a1  = v1/|v1|
+ * a2 = v1*v2/|v1*v2|
+ * a3 = a1*a2
+ * and then constructs a matrix with the ai as columns. 
+ * @param v1 The first vector
+ * @param v2 The second vector
+ * @return A matrix as descibed above.
+ */
+MATRIX matFromTwoVectors(MATRIX v1, MATRIX v2);
+#endif