- Removed -fwritable-string

SKIPPED:
	psi/dornier2.c
	psi/ecbdriv.c
	psi/el734hp.c
	psi/libpsi.a
	psi/make_gen
	psi/makefile_linux
	psi/pimotor.c
	psi/pipiezo.c
	psi/sinqhttp.c
	psi/tcpdornier.c
	psi/velodornier.c

SCinter.c

@@ -1006,6 +1006,7 @@ static void freeList(int listID)
     } while(0!=LLDnodePtr2First(listID));
     LLDdelete(listID);
 }
+
 /*------------------------------------------------------------------------*/
  void RemoveStartupCommands(void)
  {

anticollider.c

@@ -192,9 +192,12 @@ int StartLevel(int level, int sequenceList, int motorList, SConnection *pCon){
       pMot = FindMotEntry(motorList,seq.pMotor);
       if(pMot){
     	status = StartRegMot(pMot,pCon,seq.target);
-	if(status){
+        /*
+         * I have to ignore the problem here: if I do not increment the count
+         * all the other levels will not be drive and the anticollider
+         * gets into a mess
+         */
 	    count++;
-	}
       } else {
 	    sprintf(pBueffel,"ERROR: motor %s, requested from anticollider script",
 		   seq.pMotor);

anticollider.w

@@ -114,6 +114,9 @@ managed through the lld functions as everywhere within SICS.
 	pMotReg FindMotFromDataStructure(int iList, void *pData);
 	int     CheckAllMotors(int iList, SConnection *pCon);
 	void KillMotList(int iList);
+	void StopAllMotors(int iList);
+	void DeactivateAllMotors(int iList);
+	
 @}
 The functions:
 \begin{description}

cone.c

@@ -0,0 +1,356 @@
+/*----------------------------------------------------------------------
+  SICS cone module for cone scans. Form more details see conescan.tex
+  and cone.tex.
+  
+  COPYRIGHT: see file COPYRIGHT
+  
+  Mark Koennecke, March 2006
+------------------------------------------------------------------------*/
+#include <stdio.h>
+#include <assert.h>
+#include "cone.h"
+#include  "hkl.i"
+#include "vector.h"
+#include "fourlib.h"
+/*=================== Object Descriptor Interface ===================================================*/
+static void *ConeGetInterface(void *pData, int iID){
+  pConeData self = NULL;
+
+  self = (pConeData)pData;
+  if(self == NULL){
+    return NULL;
+  }
+  if(iID == DRIVEID){
+    return self->pDriv;
+  }
+  return NULL;
+}
+/*---------------------------------------------------------------------------------------------------*/
+static void ConeSaveStatus(void *data, char *name, FILE *fd){
+    pConeData self = (pConeData)data;
+    if(self == NULL){
+        return;
+    }
+    fprintf(fd,"%s center  %d\n", name,self->center);
+    fprintf(fd,"%s target %f %f %f\n", name, self->target.h,
+        self->target.k, self->target.l);
+}
+/*=================== Drivable Interface ============================================================*/
+static int ConeHalt(void *pData){
+  pConeData self = NULL;
+
+  self = (pConeData)pData;
+  assert(self != NULL);
+  
+  self->pHkl->pTheta->pDrivInt->Halt(self->pHkl->pTheta);
+  self->pHkl->pOmega->pDrivInt->Halt(self->pHkl->pOmega);
+  self->pHkl->pChi->pDrivInt->Halt(self->pHkl->pChi);
+  self->pHkl->pPhi->pDrivInt->Halt(self->pHkl->pPhi);
+  return 1;
+}
+/*-----------------------------------------------------------------------------------------------------*/
+static int ConeCheckLimits(void *self, float fVal, char *error, int errLen){
+  /*
+    There is no meaningful implementation here. This gets called when starting the motor.
+    At that stage not all other values may be known. If the calculation fails, this will die 
+    at status check time.
+  */
+  return 1;
+}
+/*------------------------------------------------------------------------*/
+static  MATRIX makeCsToPsiMatrix(reflection center, double lambda){
+    MATRIX psiToCs = NULL, csToPsi = NULL, t1, t2;
+    double z1[3], u;
+    
+  psiToCs = makeInstToConeVectorMatrix(center,lambda);
+  if(psiToCs == NULL){
+    return NULL;
+  }
+  csToPsi = mat_inv(psiToCs);
+  /*
+   * this is debugging code: remove together with variables
+   */
+  z1FromAngles(lambda,center.s2t,center.om,center.chi,center.phi,z1);
+  t1 = makeVectorInit(z1);
+  t2 = mat_mul(psiToCs,t1); 
+  normalizeVector(t2);
+  t1[0][0] = .0;
+  t1[1][0] = .0;
+  t1[2][0] = 1.;
+  u = angleBetween(t1,t2);   
+   
+  mat_free(psiToCs);
+  return csToPsi;
+}
+/*----------------------------------------------------------------------------------------------------
+ * I am lazy in this function: I calculate anew from all the data. This saves
+ * me a lot of trouble keeping track of parameter changes in UBCALC etc.
+ * ---------------------------------------------------------------------------*/ 
+static long ConeSetValue(void *pData, SConnection *pCon,  float fVal){
+  pConeData self = NULL;
+  float fSet[4];
+  double openingAngle, length, testAngle;
+  MATRIX csToPsi = NULL, B = NULL, newScat = NULL, cent;
+  int status;
+  reflection center;
+  char buffer[131];
+  double z1[3];
+
+  if(!SCMatchRights(pCon,usUser)){
+    return 0;
+  }
+
+  self = (pConeData)pData;
+  assert(self != NULL);
+
+  /*
+   * calculate opening angle
+   */
+  B = mat_creat(3,3,UNIT_MATRIX);
+  status =  calculateBMatrix(self->ubi->direct,B);
+  if(status < 0){
+    SCWrite(pCon,"ERROR: cell has no volume",eError);
+    return 0;
+  }
+  center = getReflection(self->ubi,self->center);
+  openingAngle = angleBetweenReflections(B,center,self->target);
+  
+  /*
+   * calculate conversion matrix from cone system to PSI system
+   */
+  csToPsi = makeCsToPsiMatrix(center,self->ubi->hkl->fLambda);
+  if(csToPsi == NULL){
+        SCWrite(pCon,"ERROR: bad parameters: failed to generate conversion matrix",
+            eError);
+        return 0;   
+  }
+  
+  /*
+   * calculate scattering vector on cone and make its length
+   * match the length of the apropriate scattering vector
+   */
+  length = scatteringVectorLength(B,self->target);
+  newScat = calcConeVector(openingAngle, fVal, length, csToPsi);
+  if(newScat == NULL){
+    SCWrite(pCon,"ERROR: fails to calculate cone vector",eError);
+    return 0;
+  }
+  
+  /**
+   * this is debugging code
+   */
+  length = vectorLength(newScat);
+  z1FromAngles(self->ubi->hkl->fLambda,center.s2t,center.om,center.chi,
+    center.phi,z1);
+  cent = makeVectorInit(z1);
+  testAngle = angleBetween(cent,newScat);
+  snprintf(buffer,131,"OpeningAngle = %f, testAngle = %f, length = %f",
+    openingAngle,testAngle,length);
+   SCWrite(pCon,buffer,eWarning); 
+  
+  
+  
+  /*
+   * try to find setting angles for this vector
+   */
+  status = findAllowedBisecting(self->pHkl->fLambda,newScat, fSet,
+    hklInRange, self->pHkl);
+  /*
+   * clean up matrices
+   */
+   mat_free(B);
+   mat_free(newScat);
+   mat_free(csToPsi);
+   if(status != 1){
+    SCWrite(pCon,"ERROR: cannot get cone vector into scattering position",
+        eError);
+        SCSetInterrupt(pCon,eAbortOperation);
+        return 0;
+   }
+   self->lastConeAngle = fVal;  
+  /*
+   * start motors
+   */
+  return startHKLMotors(self->pHkl, pCon,fSet);
+}
+/*---------------------------------------------------------------------------------------------------*/
+static int checkMotors(pConeData self, SConnection *pCon){
+  int status;
+
+  status = self->pHkl->pTheta->pDrivInt->CheckStatus(self->pHkl->pTheta, pCon);
+  if(status != HWIdle && status != OKOK){
+    return status;
+  }
+  status = self->pHkl->pOmega->pDrivInt->CheckStatus(self->pHkl->pOmega, pCon);
+  if(status != HWIdle && status != OKOK){
+    return status;
+  }
+  status = self->pHkl->pChi->pDrivInt->CheckStatus(self->pHkl->pChi, pCon);
+  if(status != HWIdle && status != OKOK){
+    return status;
+  }
+  status = self->pHkl->pPhi->pDrivInt->CheckStatus(self->pHkl->pPhi, pCon);
+  if(status != HWIdle && status != OKOK){
+    return status;
+  }
+  return HWIdle;
+}
+/*-----------------------------------------------------------------------------------------------------*/
+static int ConeCheckStatus(void *pData, SConnection *pCon){
+  pConeData self = NULL;
+  int status;
+
+  self = (pConeData)pData;
+  assert(self != NULL);
+  return checkMotors(self,pCon);
+}
+/*-----------------------------------------------------------------------------------------------------*/
+static float ConeGetValue(void *pData, SConnection *pCon){
+  pConeData self = NULL;
+  float fVal[3];
+  int status;
+
+  self = (pConeData)pData;
+  assert(self != NULL);
+
+  return self->lastConeAngle; 
+}
+/*=============================== Live and Death ====================================*/
+static pConeData MakeConeMot(pUBCALC u){
+  pConeData self = NULL;
+
+  assert(u != NULL);
+
+  self = (pConeData)malloc(sizeof(coneData));
+  if(self == NULL){
+    return NULL;
+  }
+  memset(self,0,sizeof(coneData));
+  self->pDes = CreateDescriptor("Cone");
+  self->pDriv = CreateDrivableInterface();
+  if(self->pDes == NULL || self->pDriv == NULL){
+    free(self);
+    return NULL;
+  }
+  self->pDes->GetInterface = ConeGetInterface;
+  self->pDriv->Halt = ConeHalt;
+  self->pDriv->CheckLimits = ConeCheckLimits;
+  self->pDriv->SetValue = ConeSetValue;
+  self->pDriv->CheckStatus = ConeCheckStatus;
+  self->pDriv->GetValue = ConeGetValue;
+  self->ubi = u;
+  self->pHkl = u->hkl;
+  return self;
+}
+/*----------------------------------------------------------------------------------*/
+static void KillConeMot(void *pData){
+  pConeData self = NULL;
+
+  self = (pConeData)pData;
+  if(self == NULL){
+    return;
+  }
+  if(self->pDes != NULL){
+    DeleteDescriptor(self->pDes);
+  }
+  if(self->pDriv){
+    free(self->pDriv);
+  }
+  free(self);
+}
+/*=============================== Interpreter Interface ============================*/
+int ConeAction(SConnection *pCon, SicsInterp *pSics, void *pData, int argc, char *argv[]){
+  pConeData self = NULL;
+  float value;
+  int id;
+  char pBuffer[132];
+
+  self = (pConeData)pData;
+  assert(self != NULL);
+
+   if(argc > 1) {
+        strtolower(argv[1]);
+        if(strcmp(argv[1],"center") == 0){
+            if(argc > 2){
+                if(!SCMatchRights(pCon,usUser)){
+                    return 0;
+                }
+                id = atoi(argv[2]);
+                if(id < 0 || id > 2 ){
+                    SCWrite(pCon,"ERROR: id must be between 0 - 3",eError);
+                    return 0;
+                }
+                self->center = id;
+                SCSendOK(pCon);
+                return 1;       
+            } else {
+                snprintf(pBuffer,131,"%s.center = %d", argv[0], self->center);
+                SCWrite(pCon,pBuffer,eValue);
+                return 1;
+            }
+        } else if (strcmp(argv[1],"target") == 0){
+            if(argc >= 5){
+                if(!SCMatchRights(pCon,usUser)){
+                    return 0;
+                }
+                self->target.h = atof(argv[2]);
+                self->target.k = atof(argv[3]);
+                self->target.l = atof(argv[4]);
+                SCSendOK(pCon);
+                return 1;
+            } else {
+                snprintf(pBuffer,131,"%s.target = %f %f %f", argv[0],
+                self->target.h, self->target.k, self->target.l);
+                SCWrite(pCon,pBuffer,eValue);
+                return 1; 
+            }
+        } else {
+            SCWrite(pCon,"ERROR: subcommand to cone not known",eError);
+            return 0;
+        }
+   }
+
+    /*
+     * default: print value
+     */
+  value = self->pDriv->GetValue(self,pCon);
+  if(value < -9000.){
+    snprintf(pBuffer,131,"ERROR: failed to read %s",argv[0]);
+    SCWrite(pCon,pBuffer,eError);
+    return 0;
+  }
+  snprintf(pBuffer,131,"%s = %f", argv[0], value);
+  SCWrite(pCon,pBuffer,eValue);
+  return 1;
+}
+/*------------------------------------------------------------------------------------------*/
+int MakeCone(SConnection *pCon, SicsInterp *pSics, void *pData, int argc, char *argv[]){
+  pUBCALC ubi = NULL;
+  pConeData pMot = NULL;
+  char pBuffer[131];
+  int status;
+
+  if(argc < 3){
+    SCWrite(pCon,"ERROR: insuffient number of arguments to MakeCone",eError);
+    return 0;
+  }
+  
+  ubi = FindCommandData(pSics,argv[2],"UBcalc"); 
+  if(ubi == NULL){
+    snprintf(pBuffer,131,"ERROR: %s is no UBcalc object" , argv[2]);
+    SCWrite(pCon,pBuffer,eError);
+    return 0; 
+  }
+  
+  pMot = MakeConeMot(ubi);
+  if(pMot == NULL){
+    SCWrite(pCon,"ERROR: out of memory creating cone virtual motor",eError);
+    return 0;
+  } 
+  status = AddCommand(pSics,argv[1],ConeAction,KillConeMot,pMot);
+  if(status != 1){
+    SCWrite(pCon,"ERROR: failed to create duplicate cone motor",eError);
+  }
+  return status;
+  
+}

cone.h

@@ -0,0 +1,32 @@
+
+/*----------------------------------------------------------------------
+  SICS cone module for cone scans. Form more details see conescan.tex
+  and cone.tex.
+  
+  COPYRIGHT: see file COPYRIGHT
+  
+  Mark Koennecke, March 2006
+------------------------------------------------------------------------*/
+#ifndef SICSCONE
+#define SICSCONE
+#include "sics.h"
+#include "ubcalc.h"
+
+/*-----------------------------------------------------------------------*/
+
+typedef struct {
+        pObjectDescriptor pDes;
+        pIDrivable pDriv;
+        reflection target;
+        float lastConeAngle;
+        pUBCALC ubi;
+        int center;
+        pHKL pHkl;
+} coneData, *pConeData;
+
+/*----------------------------------------------------------------------*/
+int MakeCone(SConnection *pCon, SicsInterp *pSics, void *pData,
+        int argc, char *argv[]);
+int ConeAction(SConnection *pCon, SicsInterp *pSics, void *pData,
+        int argc, char *argv[]);
+#endif

cone.w

@@ -0,0 +1,64 @@
+\subsection{Cone}
+This module is a virtual motor for driving on a cone around a give center reflection. 
+The application is in four circle diffractometry. Consider, a reflection has been 
+found and indexed. If the cell parameters are known, we now know that at a certain 
+two theta and at a certain agle around the first reflection there should be 
+another one. Serching this cone would speed up UB matrix determination. Driving on such
+a cone and therewith scanning it is implemented in this module. 
+
+The math for doing this is detailed in a separate internal paper: conescan.tex.
+For its implementation, cone works very closely with the UBCALC module. UBCALC 
+ provides the cell constants and it also provides the list of reflections from 
+ which to get the center reflection of the cone.
+
+The module requires a data structure:
+@d conedat @{
+typedef struct {
+	pObjectDescriptor pDes;
+	pIDrivable pDriv;
+	reflection target;
+	float lastConeAngle;
+    pUBCALC ubi;
+    int center;
+    pHKL pHkl;
+} coneData, *pConeData;
+@} 
+The fields are:
+\begin{description}
+\item[pDes] The standard object descriptor
+\item[pDriv] The drivable interface which implements most of this modules 
+ functionality.
+\item[target] The target reflection of the cone. This determines the length
+ of the scattering vector and the opening angle of the cone.
+\item[ubi] A pointer to the UBCALC module which holds things like lattice constants. 
+ \item[center] The reflection number in UBCALCS reflection list to use as the
+  center of the cone. 
+\end{description}
+
+The external interface to cone is mostly defined by the drivable interface. In 
+addition there are only the usal interpreter interface functions to install and
+configure cone.
+
+@o cone.h @{
+/*----------------------------------------------------------------------
+  SICS cone module for cone scans. Form more details see conescan.tex
+  and cone.tex.
+  
+  COPYRIGHT: see file COPYRIGHT
+  
+  Mark Koennecke, March 2006
+------------------------------------------------------------------------*/
+#ifndef SICSCONE
+#define SICSCONE
+#include "sics.h"
+#include "ubcalc.h"
+
+/*-----------------------------------------------------------------------*/
+@<conedat@>
+/*----------------------------------------------------------------------*/
+int MakeCone(SConnection *pCon, SicsInterp *pSics, void *pData,
+	int argc, char *argv[]);
+int ConeAction(SConnection *pCon, SicsInterp *pSics, void *pData,
+	int argc, char *argv[]);
+#endif
+@}

devexec.c

@@ -105,6 +105,7 @@
                    int iStop;
                    int iStatus;
                    int iEnd;
+                   int drivePrint;
                    long lTask;
                    pTaskMan pTask;
                    int iLock;
@@ -144,6 +145,7 @@
     pRes->pTask = pTask;
     pRes->lTask = -1;
     pRes->iLock = 0;
+    pRes->drivePrint = 0;
     pRes->pCall = CreateCallBackInterface();
     pRes->lastRun = time(NULL);
     return pRes;
@@ -187,6 +189,7 @@
      pICountable pCountInt = NULL;
      static int overwriteOwner = -1;
      char *overwriteOption;
+     float oldVal;
      
      assert(self);
      assert(pDes);
@@ -239,6 +242,13 @@
      if(pDrivInt)
      {
          iRet = pDrivInt->SetValue(pData,pCon,fNew);
+         if(iRet == OKOK && self->drivePrint == 1)
+         {
+            oldVal = pDrivInt->GetValue(pData,pCon);
+            snprintf(pBueffel,131,"Driving %s from %8.3f to %8.3f",
+                name, oldVal, fNew);
+            SCWrite(pCon,pBueffel,eWarning);
+         }
      }
      else if(pCountInt)
      {
@@ -289,6 +299,7 @@
      }
      return 0;
   }
+/*-----------------------------------------------------------------------*/  
 /*--------------------------------------------------------------------------*/
   int StartMotor(pExeList self, SicsInterp *pSics, SConnection *pCon,
                  char *name, float fVal)
@@ -988,6 +999,39 @@
        return 0;
      }
    }
+ /*---------------------------------------------------------------------*/
+ int DevexecAction(SConnection *pCon, SicsInterp *pSics, void *pData,
+    int argc, char *argv[])
+ {
+    int val;
+    char pBueffel[256];
+    
+    pExeList self = (pExeList)pData;
+    if(argc < 2)
+    {
+        SCWrite(pCon,"ERROR: not enough argumentd to devexec command",eError);
+        return 0;
+    }
+    strtolower(argv[1]);
+    if(strcmp(argv[1],"driveprint") == 0)
+    {
+        if(argc > 2 && SCMatchRights(pCon,usUser))
+        {
+            val = atoi(argv[2]);
+            self->drivePrint = val;
+            SCSendOK(pCon);
+            return 1;
+        }  else {
+            snprintf(pBueffel,255,"devexe.drivePrint = %d",
+                self->drivePrint);
+            SCWrite(pCon,pBueffel,eValue);
+            return 1;
+        }  
+    }  else {
+        SCWrite(pCon,"ERROR: unknown subcommand to devexec",eError);
+        return 0;
+    }   
+ }        
 /*-------------------------------------------------------------------------*/
    int ContinueAction(SConnection *pCon, SicsInterp *pSics, void *pData,
                                 int argc, char *argv[])

devexec.h

@@ -152,7 +152,11 @@
   /*
    continues execution
   */    
-   
+  int DevexecAction(SConnection *pCon, SicsInterp *pSics, void *pData,
+              int argc, char *argv[]);
+   /*
+    * various commands
+    */
 /*--------------------------- Locking ---------------------------------*/
  
 #line 183 "devexec.w"

fourlib.c

@@ -198,7 +198,7 @@ void pol2det(psdDescription *psd, double gamma, double nu, int *x, int *y){
   turn chi and phi in order to get Z1 into the equatorial plane
 */
 static void turnEquatorial(MATRIX z1, double *chi, double *phi){
-  if(ABS(z1[0][0]) < .0001 || ABS(z1[1][0]) < .0001){
+  if(ABS(z1[0][0]) < .0001 &&  ABS(z1[1][0]) < .00001){
     *phi = .0;
     *chi = 90.;
     if(z1[2][0] < .0){
@@ -607,6 +607,46 @@ void z1FromAllAngles(double lambda, double omega , double gamma,
   z1fromz3(z1,z3,chi,phi);
   mat_free(z3);    
 }
+/*-----------------------------------------------------------------------*/
+int findAllowedBisecting(double lambda, MATRIX z1, float fSet[4], 
+    inRange testFunc, void *userData){
+    int status, i, mask[4];
+    double stt, om, chi, phi, psi, ompsi, chipsi, phipsi;
+    float fTest[4];
+    
+    status = z1mToBisecting(lambda,z1, &stt, &om, &chi, &phi);
+    chi = circlify(chi);
+    phi = circlify(phi);
+    fSet[0] = stt;
+    fSet[1] = om;
+    fSet[2] = chi;
+    fSet[3] = phi;
+    if(status != 1) {
+        return 0;
+    }
+    
+    for(psi = .0; psi < 360.; psi += .5){
+        rotatePsi(om,chi,phi,psi,&ompsi,&chipsi,&phipsi);
+        fTest[0] = stt;
+        fTest[1] = ompsi;
+        fTest[2] = circlify(chipsi);
+        fTest[3] = circlify(phipsi);
+        status = testFunc(userData,fTest,mask);
+        if(status == 1){
+            for(i = 0; i < 4; i++){
+                fSet[i] = fTest[i];
+            }
+            return 1;
+        }
+        if(mask[0] == 0) {
+            /*
+             * useless: when two theta problem there is no solution
+             */
+            return 0;
+        }        
+    }
+    return 0;
+}        
 /*------------------- a test program ------------------------------------*/
 #ifdef TESTCODE
 int main(int argc, char *argv[]){

fourlib.h

@@ -172,6 +172,25 @@ void phimat(MATRIX phim, double phi);
  */
 int calcTheta(double lambda, MATRIX z1, double *d, double *theta);
 
+/**
+ * try very hard to calculate settings for the bisecting position
+ * within the instruments limits. Tries to tweak omega and to rotate 
+ * psi until a usable position can be found.
+ * @param lambda The wavelength
+ * @param z1 The scattering for which to find a position
+ * @param fSet The output setting angles. In case of a failure this 
+ * contains the normal setting angles for psi = 0. The order is:
+ * 2 theta, om, chi, phi
+ * @testFunc A user supplied function which test if the setting angles 
+ * are in range.
+ * @param userData A user specified pointer to some data which may be needed
+ *  in testing the range.
+ * @return 0 on failure, 1 on success.   
+ */   
+typedef int (*inRange)(void *userData, float dSet[4], int mask[4]);
+int findAllowedBisecting(double lambda, MATRIX z1, float fSet[4], 
+    inRange testFunc, void *userData);
+ 
 #endif   
 
 

hkl.c

@@ -597,7 +597,7 @@ static MATRIX calculateScatteringVector(pHKL self, float fHKL[3])
     return z1;
 }
 /*---------------------------------------------------------------------*/
-static int calculateBisecting(MATRIX z1, pHKL self, SConnection *pCon,
+static int calculateBisectingOld(MATRIX z1, pHKL self, SConnection *pCon,
 			      float fSet[4], double myPsi, int iRetry)
 {
   double stt, om, chi, phi, psi, ompsi, chipsi, phipsi;
@@ -674,6 +674,74 @@ static int calculateBisecting(MATRIX z1, pHKL self, SConnection *pCon,
 
   return 0;
 }
+/*----------------------------------------------------------------------*/
+int hklInRange(void *data, float fSet[4], int mask[4])
+{
+   pHKL self = (pHKL)data;
+   float fHard, fLimit;
+   char pError[132];
+   int i, test;
+   double dTheta;
+
+   /* check two theta  */
+   dTheta = fSet[0];
+   mask[0] = checkTheta(self, &dTheta);
+   fSet[0] = dTheta;
+
+   /* for omega check against the limits +- SCANBORDER in order to allow for
+      a omega scan
+   */
+   MotorGetPar(self->pOmega,"softlowerlim",&fLimit);
+   if((float)fSet[1] < fLimit + self->scanTolerance){
+     mask[1] = 0;
+   } else {
+     mask[1] = 1;
+     MotorGetPar(self->pOmega,"softupperlim",&fLimit);
+     if((float)fSet[1] > fLimit - self->scanTolerance){
+       mask[1] = 0;
+     } else {
+       mask[1] = 1;
+     }
+   }
+
+   /* check chi and phi*/
+   mask[2] = MotorCheckBoundary(self->pChi,fSet[2], &fHard,pError,131);
+   mask[3] = MotorCheckBoundary(self->pPhi,fSet[3], &fHard,pError,131);
+   for(i = 0, test = 0; i < 4; i++){ 
+        test += mask[i];
+   }
+   if(test != 4) {
+        return 0;
+   } else {
+        return 1;
+   }
+}
+/*---------------------------------------------------------------------*/
+static int calculateBisecting(MATRIX z1, pHKL self, SConnection *pCon,
+                  float fSet[4], double myPsi, int iRetry)
+{
+  double stt, om, chi, phi, psi, ompsi, chipsi, phipsi;
+  int i, test;
+
+  /* 
+     just the plain angle calculation
+  */
+  if(!z1mToBisecting(self->fLambda,z1,&stt,&om,&chi,&phi))
+  {
+    return 0;        
+  }
+  
+  if(iRetry == 1) {
+    rotatePsi(om,chi,phi,psi,&ompsi,&chipsi,&phipsi);
+    fSet[1] = ompsi;
+    fSet[2] = circlify(chipsi);
+    fSet[3] = circlify(phipsi);
+    return 1;
+  } else {
+    return   findAllowedBisecting(self->fLambda, z1, fSet, hklInRange,self);
+  }
+  
+}
 /*-----------------------------------------------------------------------*/
 static int calculateNormalBeam(MATRIX z1, pHKL self, SConnection *pCon,
 			      float fSet[4], double myPsi, int iRetry)
@@ -917,22 +985,12 @@ static void stopHKLMotors(pHKL self)
     self->pNu->pDrivInt->Halt(self->pNu);
   }
 }
-/*-------------------------------------------------------------------------*/
+/*------------------------------------------------------------------------*/
-  int RunHKL(pHKL self, float fHKL[3], 
+int startHKLMotors(pHKL self, SConnection *pCon, float fSet[4])
-             float fPsi, int iHamil, SConnection  *pCon)
 {
-     float fSet[4];
-     int iRet,i;
      char pBueffel[512];
      pDummy pDum;
-
+     int iRet;
-     assert(self);
-     iRet = CalculateSettings(self,fHKL,fPsi,iHamil,fSet,pCon);
-     if(!iRet)
-     {
-       SCWrite(pCon,"ERROR: NOT started",eError);
-       return 0;
-     }
 
      /* start all the motors */
      pDum = (pDummy)self->pTheta;
@@ -966,15 +1024,8 @@ static void stopHKLMotors(pHKL self)
             SCWrite(pCon,"ERROR: cannot start nu motor",eError);
             return 0;
         } 
-        else
-	{
-           for(i = 0; i < 3; i++)
-           {
-               self->fLastHKL[i] = fHKL[i];
-           }               
          return 1;
      }
-     }
 
      pDum = (pDummy)self->pChi;
      iRet = StartDevice(pServ->pExecutor, "HKL",
@@ -994,6 +1045,30 @@ static void stopHKLMotors(pHKL self)
         stopHKLMotors(self);
         return 0;
      }        
+     return 1;
+}
+/*-------------------------------------------------------------------------*/
+  int RunHKL(pHKL self, float fHKL[3], 
+             float fPsi, int iHamil, SConnection  *pCon)
+  {
+     float fSet[4];
+     int iRet,i;
+     char pBueffel[512];
+     pDummy pDum;
+
+     assert(self);
+     iRet = CalculateSettings(self,fHKL,fPsi,iHamil,fSet,pCon);
+     if(!iRet)
+     {
+       SCWrite(pCon,"ERROR: NOT started",eError);
+       return 0;
+     }
+
+     iRet = startHKLMotors(self,pCon,fSet);
+     if(iRet != 1){
+        return iRet;
+     }
+
      for(i = 0; i < 3; i++)
      {
        self->fLastHKL[i] = fHKL[i];

hkl.h

@@ -48,6 +48,7 @@
 
   int HKLAction(SConnection *pCon, SicsInterp *pSics, void *pData,
                 int argc, char *argv[]);  
-
+  int hklInRange(void *data, float fSet[4], int mask[4]);
+  int startHKLMotors(pHKL self, SConnection *pCon, float fSet[4]);
 
 #endif  

hkl.tex

@@ -125,7 +125,8 @@ $\langle$hklint {\footnotesize ?}$\rangle\equiv$
 \mbox{}\verb@@\\
 \mbox{}\verb@  int HKLAction(SConnection *pCon, SicsInterp *pSics, void *pData,@\\
 \mbox{}\verb@                int argc, char *argv[]);  @\\
-\mbox{}\verb@@\\
+\mbox{}\verb@  int hklInRange(void *data, float fSet[4], int mask[4]);@\\
+\mbox{}\verb@  int startHKLMotors(pHKL self, SConnection *pCon, float fSet[4]);@\\
 \mbox{}\verb@@$\diamond$
 \end{list}
 \vspace{-1ex}

hkl.w

@@ -108,7 +108,8 @@ module:
 
   int HKLAction(SConnection *pCon, SicsInterp *pSics, void *pData,
                 int argc, char *argv[]);  
-
+  int hklInRange(void *data, float fSet[4], int mask[4]);
+  int startHKLMotors(pHKL self, SConnection *pCon, float fSet[4]);
 @}
 All functions return 0 on failure, 1 on success if not stated otherwise.
 Most functions take a pointer to a HKL data structure as first parameter.

make_gen

@@ -30,7 +30,7 @@ SOBJ =  network.o ifile.o conman.o SCinter.o splitter.o passwd.o \
        s_rnge.o sig_die.o gpibcontroller.o $(NIOBJ) mcreader.o mccontrol.o\
        hmdata.o  nxscript.o  tclintimpl.o sicsdata.o mcstascounter.o \
        mcstashm.o initializer.o remob.o tclmotdriv.o protocol.o \
-       sinfox.o sicslist.o
+       sinfox.o sicslist.o cone.o
  
 MOTOROBJ = motor.o  simdriv.o  
 COUNTEROBJ =  countdriv.o simcter.o counter.o

makefile_slinux

@@ -16,8 +16,8 @@ include sllinux_def
 CC       = gcc
 CFLAGS = -I$(HDFROOT)/include -DHDF4 -DHDF5 -DNXXML $(NI) \
        -Ipsi/hardsup -I. \
-       -fwritable-strings -DCYGNUS -DNONINTF -g $(DFORTIFY) \
+       -Werror -DCYGNUS -DNONINTF -g $(DFORTIFY) \
-       -Wall -Wno-unused -Wno-comment -Wno-switch -Werror
+       -Wall -Wno-unused -Wno-comment -Wno-switch 
 
 BINTARGET = bin
 EXTRA=nintf.o

mcstas/dmc/DataNumber

@@ -1,3 +1,3 @@
- 119 
+ 120 
 NEVER, EVER modify or delete this file
 You'll risk eternal damnation and a reincarnation as a cockroach!|n

mcstas/dmc/vdmc.tcl

@@ -6,7 +6,7 @@
 #---------------------------------------------------------------------------
 #                        O P T I O N S
 
-set home $env(HOME)/psi/sics/mcstas/dmc
+set home $env(HOME)/src/workspace/sics/mcstas/dmc
 
 #--------------------------------- first all the server options are set
 #ServerOption RedirectFile $home/stdcdmc
@@ -139,11 +139,11 @@ commandlog auto
 commandlog intervall 5
 
 #----------- enable sycamore
-InstallProtocolHandler
+#InstallProtocolHandler
-InstallSinfox
+#InstallSinfox
-source sycFormat.tcl
+#source sycFormat.tcl
-source /usr/lib/tcllib1.6.1/stooop/stooop.tcl
+#source /usr/lib/tcllib1.6.1/stooop/stooop.tcl
-namespace import stooop::*
+#namespace import stooop::*
-source sinfo.tcl
+#source sinfo.tcl
-source sycamore.tcl
+#source sycamore.tcl
-Publish sinfo Spy
+#Publish sinfo Spy

mcstas/dmc/vdmcstatus.tcl

@@ -98,10 +98,10 @@ twothetad ignorefault 0.000000
 twothetad AccessCode 2.000000
 twothetad movecount 10.000000
 # Counter counter
-counter SetPreset 60000.000000
+counter SetPreset 60000000.000000
 counter SetMode Monitor
 banana CountMode monitor
-banana preset 60.000000
+banana preset 60000.000000
 # Motor a1
 a1 sign 1.000000
 a1 SoftZero 0.000000
@@ -215,7 +215,7 @@ comment2 UNKNOWN
 comment2 setAccess 2
 comment3 UNKNOWN
 comment3 setAccess 2
-starttime 2005-12-20 05:13:05
+starttime 2006-03-10 15:33:04
 starttime setAccess 2
 adress UNKNOWN
 adress setAccess 2

motreg.c

@@ -115,7 +115,14 @@ int StartRegMot(pMotReg self, SConnection *pCon, float fValue){
   */
 
   pDriv->SetValue = oldSet;
+  if(ret == 1){
       self->iActive = 1;
+  } else {
+      snprintf(pBueffel,131,"ERROR: failed to start motor %s",
+        self->motorName);
+      SCWrite(pCon,pBueffel,eError);
+      self->iActive = 0;
+  }
   return ret;
 }
 /*----------------------------------------------------------------------*/

motreglist.c

@@ -91,6 +91,21 @@ void StopAllMotors(int iList){
   }
 }
 /*----------------------------------------------------------------------*/
+void DeactivateAllMotors(int iList){
+  int iRet, count = 0;
+  pMotReg pMot = NULL;
+  pIDrivable pDriv;
+
+  iRet = LLDnodePtr2First(iList);
+  while(iRet != 0){
+    LLDnodeDataTo(iList,&pMot);
+    if(pMot != NULL){
+      pMot->iActive = 0;
+    }
+    iRet = LLDnodePtr2Next(iList);
+  }
+}
+/*----------------------------------------------------------------------*/
 void KillMotList(int iList){
   int iRet;
   pMotReg pMot = NULL;

motreglist.h

@@ -20,6 +20,9 @@
         int     CheckAllMotors(int iList, SConnection *pCon);
         void KillMotList(int iList);
         void StopAllMotors(int iList);
+        void DeactivateAllMotors(int iList);
+        
+
 
 
 #endif

napi4.c

@@ -21,7 +21,7 @@
              
   For further information, see <http://www.neutron.anl.gov/NeXus/>
   
-  $Id: napi4.c,v 1.7 2006/03/23 12:39:11 zolliker Exp $
+  $Id: napi4.c,v 1.8 2006/03/31 15:24:53 koennecke Exp $
 
 ----------------------------------------------------------------------------*/
 #include <stdlib.h>
@@ -1236,7 +1236,8 @@ extern	void *NXpData;
   {
     pNexusFile pFile;
     pFile = NXIassert (fid);
-     printf("HDF4 link: iTag = %ld, iRef = %ld, target=\"%s\"\n", sLink->iTag, sLink->iRef, sLink->targetPath);
+     printf("HDF4 link: iTag = %ld, iRef = %ld, target=\"%s\"\n", 
+             sLink->iTag, sLink->iRef, sLink->targetPath);
     return NX_OK;
   }
   

ofac.c

@@ -117,6 +117,7 @@
 #include "protocol.h"
 #include "sinfox.h"
 #include "sicslist.h"
+#include "cone.h"
 /*----------------------- Server options creation -------------------------*/
   static int IFServerOption(SConnection *pCon, SicsInterp *pSics, void *pData,
                      int argc, char *argv[])
@@ -252,6 +253,7 @@
      AddCommand(pInter,"Success",Success,NULL,pExe);
      AddCommand(pInter,"pause",PauseAction,NULL,pExe);
      AddCommand(pInter,"continue",ContinueAction,NULL,pExe);
+     AddCommand(pInter,"devexec",DevexecAction,NULL,pExe);
 
     /* add additional object creation commands here */
     AddCommand(pInter,"TokenInit",TokenInit,NULL,NULL);
@@ -319,6 +321,8 @@
 	       InstallProtocol,NULL,NULL);
      AddCommand(pInter,"InstallSinfox",
            InstallSinfox,NULL,NULL);
+     AddCommand(pInter,"MakeCone",
+           MakeCone,NULL,NULL);
 
     /*
       install site specific commands
@@ -386,6 +390,7 @@
     RemoveCommand(pSics,"MakemcStasReader");
     RemoveCommand(pSics,"InstallProtocolHandler");
     RemoveCommand(pSics,"InstallSinfox");
+    RemoveCommand(pSics,"MakeCone"); 
     /*
       remove site specific installation commands
     */

sicsstat.tcl

@@ -3,3 +3,5 @@ counter SetPreset 10.000000
 counter SetMode Timer
 hm CountMode timer
 hm preset 10.000000
+hm genbin 0.000000 10.000000 10000
+hm init

tasdrive.c

@@ -190,9 +190,20 @@ static int TASHalt(void *pData){
   }
   return 1;
 }
+
+/*--------------------------------------------------------------------------*/
+static void writeMotPos(SConnection *pCon, char *name, 
+    float val, float target){
+    char pBueffel[132];
+    
+    snprintf(pBueffel,131,"Driving %5s from %8.3f to %8.3f",
+        name, val, target);
+    SCWrite(pCon,pBueffel,eWarning);       
+    
+}
 /*---------------------------------------------------------------------------*/
 static int startMotors(ptasMot self, tasAngles angles, 
-		       SConnection *pCon, int driveQ){
+		       SConnection *pCon, int driveQ, int driveTilt){
   float val;
   double curve;
   int status;
@@ -209,6 +220,8 @@ static int startMotors(ptasMot self, tasAngles angles,
       return status;
     }
   }
+  writeMotPos(pCon,"a1",val, angles.monochromator_two_theta/2.);
+  
   val = self->math->motors[A2]->pDrivInt->GetValue(self->math->motors[A2],pCon);
   if(ABS(val - angles.monochromator_two_theta) > MOTPREC){
     status = self->math->motors[A2]->pDrivInt->SetValue(self->math->motors[A2],
@@ -218,6 +231,8 @@ static int startMotors(ptasMot self, tasAngles angles,
       return status;
     }
   }
+  writeMotPos(pCon,"a2",val, angles.monochromator_two_theta);
+  
   if(self->math->motors[MCV] != NULL){
     curve = maCalcVerticalCurvature(self->math->machine.monochromator,
 				    angles.monochromator_two_theta);
@@ -230,7 +245,9 @@ static int startMotors(ptasMot self, tasAngles angles,
 	   return status;
       }
     }
+    writeMotPos(pCon,"mcv",val, curve);
    }
+   
   if(self->math->motors[MCH] != NULL){
     curve = maCalcHorizontalCurvature(self->math->machine.monochromator,
 				    angles.monochromator_two_theta);
@@ -243,6 +260,7 @@ static int startMotors(ptasMot self, tasAngles angles,
 	return status;
       }
     }
+    writeMotPos(pCon,"mch",val, curve);
   }
 
   /*
@@ -259,6 +277,9 @@ static int startMotors(ptasMot self, tasAngles angles,
       		return status;
     	}
   	}
+    writeMotPos(pCon,self->math->motors[A5]->name,
+        val, angles.analyzer_two_theta/2.);
+    
   	val = self->math->motors[A6]->pDrivInt->GetValue(self->math->motors[A6],pCon);
   	if(ABS(val - angles.analyzer_two_theta) > MOTPREC){
     	status = self->math->motors[A6]->pDrivInt->SetValue(self->math->motors[A6],
@@ -268,6 +289,8 @@ static int startMotors(ptasMot self, tasAngles angles,
       		return status;
     	}
   	}
+    writeMotPos(pCon,"a6",val, angles.analyzer_two_theta);
+
   	if(self->math->motors[ACV] != NULL){
     	curve = maCalcVerticalCurvature(self->math->machine.analyzer,
 						    angles.analyzer_two_theta);
@@ -280,6 +303,7 @@ static int startMotors(ptasMot self, tasAngles angles,
 			return status;
      	 }
   	  }
+    writeMotPos(pCon,"acv",val, curve);
   	}
  	if(self->math->motors[ACH] != NULL){
 	    curve = maCalcHorizontalCurvature(self->math->machine.analyzer,
@@ -295,6 +319,7 @@ static int startMotors(ptasMot self, tasAngles angles,
  	     }
  	   }
  	 }
+    writeMotPos(pCon,"ach",val, curve);
   }
   
   if(driveQ == 0){
@@ -313,6 +338,8 @@ static int startMotors(ptasMot self, tasAngles angles,
       return status;
     }
   }
+  writeMotPos(pCon,"a3",val, angles.a3);
+  
   val = self->math->motors[A4]->pDrivInt->GetValue(self->math->motors[A4],pCon);
   if(ABS(val - angles.sample_two_theta) > MOTPREC){
     status = self->math->motors[A4]->pDrivInt->SetValue(self->math->motors[A4],
@@ -322,6 +349,9 @@ static int startMotors(ptasMot self, tasAngles angles,
       return status;
     }
   }
+  writeMotPos(pCon,"a4",val, angles.sample_two_theta);
+  
+  if(driveTilt == 1){
       val = self->math->motors[SGL]->pDrivInt->GetValue(self->math->motors[SGL],pCon);
     if(ABS(val - angles.sgl) > MOTPREC){
         status = self->math->motors[SGL]->pDrivInt->SetValue(self->math->motors[SGL],
@@ -331,6 +361,8 @@ static int startMotors(ptasMot self, tasAngles angles,
             return status;
         }
     }
+    writeMotPos(pCon,"sgl",val, angles.sgl);
+  
     val = self->math->motors[SGU]->pDrivInt->GetValue(self->math->motors[SGU],pCon);
     if(ABS(val - angles.sgu) > MOTPREC){
         status = self->math->motors[SGU]->pDrivInt->SetValue(self->math->motors[SGU],
@@ -340,12 +372,14 @@ static int startMotors(ptasMot self, tasAngles angles,
             return status;
         }
     }
+    writeMotPos(pCon,"sgu",val, angles.sgu);
+  }
   self->math->mustDrive = 0;
   return OKOK;
 }
 /*---------------------------------------------------------------------------*/
 static int checkQMotorLimits(ptasMot self, SConnection *pCon, 
-			     tasAngles angles){
+			     tasAngles angles, int driveTilt){
   int status, retVal = 1;
   char error[131];
   char pBueffel[256];
@@ -370,6 +404,7 @@ static int checkQMotorLimits(ptasMot self, SConnection *pCon,
     SCWrite(pCon,pBueffel,eError);
   }
 
+  if(driveTilt == 1){
       status = self->math->motors[A3]->pDrivInt->CheckLimits(self->math->motors[SGU],
 	   						 angles.sgu,
 		  					 error,
@@ -389,12 +424,14 @@ static int checkQMotorLimits(ptasMot self, SConnection *pCon,
             snprintf(pBueffel,256,"ERROR: limit violation an SGL: %s", error);
             SCWrite(pCon,pBueffel,eError);
       }
+  }
   return retVal;
 }
 /*-----------------------------------------------------------------------------*/
 static int calculateAndDrive(ptasMot self, SConnection *pCon){
   tasAngles angles;
-  int status, driveQ = 1;
+  int status, driveQ = 1, driveTilt = 1;
+  MATRIX scatteringPlaneNormal = NULL;
 
   if(self->math->ubValid == 0){
     SCWrite(pCon,"WARNING: UB matrix invalid",eWarning);
@@ -424,14 +461,35 @@ static int calculateAndDrive(ptasMot self, SConnection *pCon){
     driveQ = 0;
     break;
   }
-  if(!checkQMotorLimits(self,pCon,angles)){
+  /**
+   * check out of plane condition and permission
+   */
+  if(self->math->outOfPlaneAllowed != 1){
+    scatteringPlaneNormal = calcScatteringPlaneNormal(self->math->r1.qe,
+        self->math->r2.qe);
+    if(scatteringPlaneNormal == NULL){
+        SCWrite(pCon,"ERROR: unable to calculate scattering plane normal",eError);
         driveQ = 0;
+    } else {   
+        if(!isInPlane(scatteringPlaneNormal,self->math->target)){
+            SCWrite(pCon,
+            "WARNING: scattering vector is out of plane and you disallowed me to try to drive there",
+         eWarning);
+        }
+    }
+    driveTilt = 0;
+  }
+  
+  if(driveQ != 0){
+      if(!checkQMotorLimits(self,pCon,angles,driveTilt)){
+        driveQ = 0;
+    }
   }
 
   if(driveQ == 0){
     SCWrite(pCon,"WARNING: NOT driving Q-vector because of errors",eError);
   }
-  return startMotors(self,angles,pCon, driveQ);
+  return startMotors(self,angles,pCon, driveQ,driveTilt);
 }
 /*-----------------------------------------------------------------------------*/
 static int checkMotors(ptasMot self, SConnection *pCon){

tasub.c

@@ -65,6 +65,7 @@ static int tasUBSave(void *pData, char *name, FILE *fd){
 	  self->cell.b, self->cell.c, self->cell.alpha,
 	  self->cell.beta, self->cell.gamma);
   saveReflections(self,name,fd);
+  fprintf(fd,"%s outofplane %d\n", name, self->outOfPlaneAllowed);
   if(self->tasMode == KICONST){
      fprintf(fd,"%s const  ki\n",name);
   }else if(self->tasMode == ELASTIC){
@@ -136,6 +137,7 @@ static ptasUB MakeTasUB(){
   pNew->tasMode = KICONST;
   pNew->targetEn = .0;
   pNew->actualEn = .0;
+  pNew->outOfPlaneAllowed = 1;
   pNew->mustRecalculate = 1;
 
   return pNew;
@@ -1468,6 +1470,25 @@ int TasUBWrapper(SConnection *pCon,SicsInterp *pSics, void *pData,
       SCWrite(pCon,pBueffel,eValue);
       return 1;
     }
+  } else  if(strcmp(argv[1],"outofplane") == 0){
+    if(argc > 2){
+      strtolower(argv[2]);
+      if(!SCMatchRights(pCon,usUser)){
+        return 0;
+      }
+      status = Tcl_GetInt(InterpGetTcl(pSics),argv[2],&newSS);
+      if(status != TCL_OK){
+        SCWrite(pCon,"ERROR: failed to convert argument to number",eError);
+        return 0;
+      }
+      self->outOfPlaneAllowed = newSS;
+      SCSendOK(pCon);
+      return 1;
+    } else {
+      snprintf(pBueffel,131,"%s.outofplane = %d",argv[0],self->outOfPlaneAllowed);
+      SCWrite(pCon,pBueffel,eValue);
+      return 1;
+    }
   } else {
     snprintf(pBueffel,131,"ERROR: subcommand %s to %s not defined",argv[1],
 	     argv[0]);

tasub.h

@@ -23,6 +23,7 @@
                 tasQEPosition target;
                 tasQEPosition current;
                 int tasMode;
+                int outOfPlaneAllowed;
                 double targetEn, actualEn;
                 int mustRecalculate;
                 int mustDrive; 

tasub.w

@@ -25,11 +25,13 @@ A data structure:
 		tasQEPosition target;
 		tasQEPosition current;
 		int tasMode;
+		int outOfPlaneAllowed;
 		double targetEn, actualEn;
 		int mustRecalculate;
 		int mustDrive; 
 		pMotor motors[12];
-		int r1, r2;
+		tasReflection r1, r2;
+		int ubValid;
 }tasUB, *ptasUB;
 @} 
 \begin{description}
@@ -41,6 +43,7 @@ A data structure:
 \item[target] The Q energy target position
 \item[current] The actual Q energy position as calculated from angles.
 \item[tasMode] The mode: constant KI or constant KF
+\item[outOfPlaneAllowed] 0/1 determining if it is allowed to drive out of plane or not.
 \item[ptargetEn] The target energy transfer.
 \item[actualEn] The actual energy transfer as calcluated from angles.
 \item[mustRecalculate]  A flag indicatin that the current Q energy psoition has to be 

tasublib.c

@@ -20,6 +20,7 @@
 #define DEGREE_RAD      (PI/180.0)      /* Radians per degree */
 #define VERT 0
 #define HOR  1
+#define INPLANEPREC 0.01
 /*============== monochromator/analyzer stuff =========================*/
 double energyToK(double energy){
   double K;
@@ -178,6 +179,9 @@ MATRIX calcPlaneNormal(tasReflection r1, tasReflection r2){
 	planeNormal[i][0] = -1.*planeNormal[i][0];
       }
     }
+    mat_free(u1);
+    mat_free(u2);
+    normalizeVector(planeNormal);
     return planeNormal;
   } else {
     return NULL;
@@ -438,7 +442,7 @@ int calcTasQAngles(MATRIX UB, MATRIX planeNormal, int ss,  tasQEPosition qe,
   
   angles->a3 = om + theta;
   /*
-    put a3 into -180, 180 properly. We cal always turn by 180 because the
+    put a3 into -180, 180 properly. We can always turn by 180 because the
     scattering geometry is symmetric in this respect. It is like looking at
     the scattering plane from the other side
   */
@@ -674,3 +678,44 @@ double getTasPar(tasQEPosition qe, int tasVar){
     assert(0);
   }
 }
+/*-------------------------------------------------------------------------*/
+int isInPlane(MATRIX planeNormal, tasQEPosition qe){
+    MATRIX v = NULL;
+    double val;
+    
+    v = makeVector();
+    v[0][0] = qe.qh;
+    v[1][0] = qe.qk;
+    v[2][0] = qe.ql;
+    val = vectorDotProduct(planeNormal,v);
+    mat_free(v);
+    if(ABS(val) > INPLANEPREC){
+        return 0;
+    } else {
+        return 1;
+    }
+}
+/*--------------------------------------------------------------------*/
+MATRIX calcScatteringPlaneNormal(tasQEPosition qe1, tasQEPosition qe2){
+    MATRIX v1 = NULL, v2 = NULL, planeNormal;
+    
+    v1 = makeVector(); 
+    v2 = makeVector();
+    if(v1 != NULL && v2 != NULL){
+        v1[0][0] = qe1.qh;
+        v1[1][0] = qe1.qk;
+        v1[2][0] = qe1.ql;
+
+        v2[0][0] = qe2.qh;
+        v2[1][0] = qe2.qk;
+        v2[2][0] = qe2.ql;
+    
+        planeNormal = vectorCrossProduct(v1,v2);
+        normalizeVector(planeNormal);
+        mat_free(v1);
+        mat_free(v2);
+        return planeNormal;
+    } else {
+        return NULL;
+    }
+}

tasublib.h

@@ -14,7 +14,7 @@
 /*================= error codes =====================================*/
 #define ENERGYTOBIG -700
 #define BADSYNC     -701  /* mono/analyzer out of sync: 2*theta != two_theta*/
-#define UBNOMEMORY  -702
+#define UBNOMEMORY  -200
 #define TRIANGLENOTCLOSED  -703
 #define BADRMATRIX         -704
 #define BADUBORQ           -705
@@ -224,5 +224,20 @@ void setTasPar(ptasQEPosition qe, int tasMode, int tasVar, double value);
  * @return The value of the TAS variable.
  */
 double getTasPar(tasQEPosition qe, int tasVar);
+/**
+ * checks if a QE Position is in the scattering plane as defined by the 
+ * planeNormal
+ * @param planeNormal The plane normal of the scattering plane
+ * @param qe The QE position to check
+ * @return 0 when not in plane,1 when in plane
+ */
+int isInPlane(MATRIX scatteringPlaneNormal, tasQEPosition qe);  
+/**
+ * calculate the normal of the scattering plane from two reflections
+ * @param qe1 QE Position of first reflection in scattering plane
+ * @param qe2 QE position of second reflection in scattering plane
+ * @return The scattering plane normal
+ */
+MATRIX calcScatteringPlaneNormal(tasQEPosition qe1, tasQEPosition qe2);
  
 #endif

ubcalc.c

@@ -17,17 +17,6 @@
 #include "ubcalc.h"
 #include "motor.h"
 #include "hkl.h"
-/*---------------------------------------------------------------------*/
-typedef struct {
-  pObjectDescriptor pDes;
-  pHKL hkl;
-  lattice direct;
-  reflection r1, r2, r3;
-  MATRIX UB;
-  double allowedDeviation;
-  int indexSearchLimit;
-  int maxSuggestions;
-} UBCALC, *pUBCALC;
 /*----------------------------------------------------------------------*/
 static void killUBCALC(void *pData){
   pUBCALC self = (pUBCALC)pData;
@@ -610,3 +599,24 @@ int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData,
   }
   return 1;
 }
+/*------------------------------------------------------------------------*/
+reflection getReflection(void  *ubcalc, int no){
+    pUBCALC self = (pUBCALC)ubcalc;
+    
+    assert(self != NULL);
+    
+    switch(no){
+        case 0:
+            return self->r1;
+            break;
+         case 1:
+            return self->r2;
+            break;
+         case 2:
+            return self->r3;
+            break;
+          default:
+            assert(0);
+            break;
+    }
+}

ubcalc.h

@@ -1,6 +1,6 @@
 
 /*----------------------------------------------------------------------
-  UB caclualtion routines for four circle diffraction.
+  UB calculation routines for four circle diffraction.
   This is the interpreter interface to functionality implemented
   in fourlib.c
 
@@ -10,10 +10,27 @@
 -----------------------------------------------------------------------*/
 #ifndef SICSUBCALC
 #define SICSUBCALC
-
+#include "sics.h"
+#include "matrix/matrix.h"
+#include "cell.h"
+#include "hkl.h"
+#include "ubfour.h"
+/*---------------------------------------------------------------------*/
+typedef struct {
+  pObjectDescriptor pDes;
+  pHKL hkl;
+  lattice direct;
+  reflection r1, r2, r3;
+  MATRIX UB;
+  double allowedDeviation;
+  int indexSearchLimit;
+  int maxSuggestions;
+} UBCALC, *pUBCALC;
+/*-------------------------------------------------------------------*/
 int MakeUBCalc(SConnection *pCon,SicsInterp *pSics, void *pData,                
                 int argc, char *argv[]);
 int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData,            
                 int argc, char *argv[]);
                 
+reflection getReflection(void *ubcalc, int no);
 #endif   

ubcalc.w

@@ -7,7 +7,7 @@ This module helps in the calculation of UB matrices for four
 
 @o ubcalc.h @{
 /*----------------------------------------------------------------------
-  UB caclualtion routines for four circle diffraction.
+  UB calculation routines for four circle diffraction.
   This is the interpreter interface to functionality implemented
   in fourlib.c
 
@@ -17,12 +17,29 @@ This module helps in the calculation of UB matrices for four
 -----------------------------------------------------------------------*/
 #ifndef SICSUBCALC
 #define SICSUBCALC
-
+#include "sics.h"
+#include "matrix/matrix.h"
+#include "cell.h"
+#include "hkl.h"
+#include "ubfour.h"
+/*---------------------------------------------------------------------*/
+typedef struct {
+  pObjectDescriptor pDes;
+  pHKL hkl;
+  lattice direct;
+  reflection r1, r2, r3;
+  MATRIX UB;
+  double allowedDeviation;
+  int indexSearchLimit;
+  int maxSuggestions;
+} UBCALC, *pUBCALC;
+/*-------------------------------------------------------------------*/
 int MakeUBCalc(SConnection *pCon,SicsInterp *pSics, void *pData,		
 		int argc, char *argv[]);
 int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData,		
 		int argc, char *argv[]);
 		
+reflection getReflection(void *ubcalc, int no);
 #endif   
 @} 
  

ubfour.c

@@ -10,6 +10,7 @@
  *
  * Mark Koennecke, March 2005
  */
+#include <math.h>
 #include "ubfour.h"
 #include <assert.h>
 #include "vector.h"
@@ -376,3 +377,137 @@ int searchIndex(lattice direct, double lambda, double two_theta, double max_devi
   LLDdelete(list);
   return status;
 } 
+/*-------------------------------------------------------------------------*/
+double angleBetweenReflections(MATRIX B, reflection r1, reflection r2){
+    MATRIX chi1, chi2, h1, h2;
+    double angle;
+    
+    h1 = makeVector();
+    if(h1 == NULL){
+        return -9999.99;
+    }
+    h1[0][0] = r1.h;
+    h1[1][0] = r1.k;
+    h1[2][0] = r1.l;
+
+    h2 = makeVector();
+    if(h2 == NULL){
+        return -999.99;
+    }
+    h2[0][0] = r2.h;
+    h2[1][0] = r2.k;
+    h2[2][0] = r2.l;
+    
+    chi1 = mat_mul(B,h1);
+    chi2 = mat_mul(B,h2);
+    if(chi1 != NULL && chi2 != NULL){
+        angle = angleBetween(chi1,chi2);
+        killVector(chi1);
+        killVector(chi2);
+    }
+    killVector(h1);
+    killVector(h2);
+    return angle;
+}
+/*---------------------------------------------------------------------------*/
+MATRIX makeInstToConeVectorMatrix(reflection r,double lambda){
+    double z1[3], alpha, beta;
+    MATRIX mAlpha = NULL, mBeta = NULL, inst2CS = NULL;
+    
+    z1FromAngles(lambda,r.s2t,r.om,r.chi,r.phi,z1);
+    alpha =  atan2(z1[1],z1[0]);
+    beta  =  -atan2(z1[0],z1[2]);
+/*     printf("alpha = %f, beta = %f\n", alpha*57.57, beta*57.57);*/
+    
+    mAlpha = mat_creat(3,3,ZERO_MATRIX);
+    mBeta = mat_creat(3,3,ZERO_MATRIX);
+    if(mAlpha == NULL || mBeta == NULL){
+        return NULL;
+    }
+    mAlpha[0][0] = cos(alpha);
+    mAlpha[0][1] = sin(alpha);
+    mAlpha[1][0] = -sin(alpha);
+    mAlpha[1][1] = cos(alpha);
+    mAlpha[2][2] = 1.;
+    
+    mBeta[0][0] = cos(beta);
+    mBeta[0][2] = sin(beta);
+    mBeta[1][1] = 1.;
+    mBeta[2][0] = -sin(beta);
+    mBeta[2][2] =  cos(beta);
+    
+    inst2CS = mat_mul(mBeta,mAlpha);
+    mat_free(mAlpha);
+    mat_free(mBeta);
+    
+    return inst2CS;
+}
+/*--------------------------------------------------------------------------*/
+MATRIX calcConeVector(double openingAngle, double coneAngle, 
+    double length, MATRIX coneToPsi){
+    MATRIX coneRot = NULL, nullVector = NULL, coneVector = NULL;
+    MATRIX coneVectorScatter = NULL;
+    double testAngle;
+    MATRIX z;
+    
+    z = makeVector();
+    z[2][0] =  1.;
+ 
+    
+    coneRot = mat_creat(3,3,ZERO_MATRIX);
+    if(coneRot == NULL){
+        return NULL;
+    }
+    coneRot[0][0] = Cosd(coneAngle);
+    coneRot[0][1] = -Sind(coneAngle);
+    coneRot[1][0] = Sind(coneAngle);
+    coneRot[1][1] = Cosd(coneAngle);
+    coneRot[2][2] = 1.0;
+    
+    nullVector = makeVector();
+    if(nullVector == NULL){
+        return NULL;
+    }
+    nullVector[0][0] = Sind(openingAngle);
+    nullVector[1][0] = .0;
+    nullVector[2][0] = Cosd(openingAngle);
+    normalizeVector(nullVector);
+    scaleVector(nullVector,length);
+    testAngle = angleBetween(z,nullVector);
+    
+    coneVector = mat_mul(coneRot,nullVector);
+    if(coneVector == NULL){
+        return NULL;
+    }
+    testAngle = angleBetween(z,coneVector);
+    
+    coneVectorScatter = mat_mul(coneToPsi,coneVector);
+    
+    mat_free(coneRot);
+    killVector(nullVector);
+    killVector(coneVector);
+    
+    return coneVectorScatter;
+}
+/*---------------------------------------------------------------------------*/
+double scatteringVectorLength(MATRIX B, reflection r){
+    MATRIX h = NULL, psi = NULL;
+    double length;
+    
+    h = makeVector();
+    if(h == NULL){
+        return -9999.9;
+    }
+    h[0][0] = r.h;
+    h[1][0] = r.k;
+    h[2][0] = r.l;
+    
+    psi = mat_mul(B,h);
+    
+    length = vectorLength(psi);
+    killVector(h);
+    killVector(psi);
+    return  length;
+}
+    
+ 

ubfour.h

@@ -5,9 +5,14 @@
  *
  * Implemented:
  * - UB from cell cell constants and two reflections.
+ * - UB from three reflections
  * - Brute force index search
  *
- * Mark Koennecke, march 2005
+ * Mark Koennecke, March 2005
+ * 
+ * Added some general crystallographic calculations and stuff for conus scans 
+ * 
+ * Mark Koennecke, March 2006
  */
 #ifndef SICSUBFOUR
 #define SICSUBFOUR
@@ -73,4 +78,39 @@ typedef struct {
  */
 int searchIndex(lattice direct, double lambda, double two_theta, double max_deviation,
 		int limit, refIndex index[], int maxIndex);  
+/**
+ * calculate the angle between two reflections, given their miller indices
+ * @param B The B metric matrix
+ * @param r1 first reflection
+ * @param r2 second reflection
+ * @return angle between reflections
+ */        
+double angleBetweenReflections(MATRIX B, reflection r1, reflection r2);
+/**
+ * calculate the length of the scattering vector belonging to r
+ * @param B The B metric matrix to use
+ * @param r The reflction for wihic to calculate
+ * @return The length of the scattering vector
+ */
+double scatteringVectorLength(MATRIX B, reflection r);
+/**
+ * build the conversion MATRIX from the Busing Levy psi system (z = UB*h)
+ * to the coordinate system of a given center reflection for a cone
+ * @param r The reflection around which the cone is situated
+ * @param lambda The wavelength
+ * @return A sutiable conversion matrix or NULL in case of errors
+ */
+MATRIX makeInstToConeVectorMatrix(reflection r, double lambda);
+/**
+ * calculate a scattering vector on a cone around a given reflection at a
+ * specified cone rotation angle and a cone opening angle. The center 
+ * reflection of the cone is hidden in the conversion matrix.
+ * @param openingAngle The opening angle of the cone
+ * @param coneAngle The angle on the cone
+ * @param coneToPsi The matrix for the conversion from the cone coordinate
+ * system to the psi instrument coordinate system.
+ * @return a scattering vector on the cone
+ */
+MATRIX calcConeVector(double openingAngle, double coneAngle, 
+    double length, MATRIX coneToPsi);
 #endif 

vector.c

@@ -12,6 +12,8 @@
 #include <assert.h>
 #include <math.h>
 #include "vector.h"
+#include "trigd.h"
+
 /*----------------------------------------------------------------------*/
 MATRIX makeVector(){
   return mat_creat(3,1,ZERO_MATRIX);
@@ -143,3 +145,26 @@ MATRIX matFromTwoVectors(MATRIX v1, MATRIX v2){
   killVector(a3);
   return result;
 }
+/*-----------------------------------------------------------------------*/
+double angleBetween(MATRIX v1, MATRIX v2){
+    double angle, angles;
+    MATRIX v3 = NULL;
+    
+    angle = vectorDotProduct(v1,v2)/(vectorLength(v1) * vectorLength(v2));
+    v3 = vectorCrossProduct(v1,v2);
+    if(v3 != NULL){
+        angles = vectorLength(v3)/(vectorLength(v1) * vectorLength(v2));
+        angle  = Atan2d(angles,angle);
+    } else {
+        angle = Acosd(angle);
+    }
+    return angle;
+}
+/*----------------------------------------------------------------------*/
+void scaleVector(MATRIX v, double scale){
+    int i;
+    
+    for(i = 0; i < 3; i++){
+        v[i][0]  *= scale;
+    }
+}

vector.h

@@ -90,4 +90,18 @@ MATRIX vectorCrossProduct(MATRIX v1, MATRIX v2);
  * @return A matrix as descibed above.
  */
 MATRIX matFromTwoVectors(MATRIX v1, MATRIX v2);
+/**
+ * calculate the nagle between two vectors
+ * @param v1 first vector
+ * @param v2 second vector
+ * @return angle in degree
+ */
+double angleBetween(MATRIX v1, MATRIX v2);
+
+/**
+ * scale the vector v with scale
+ * @param v The vector to scale
+ * @param scale the scaling factor apply
+ */
+void scaleVector(MATRIX v, double scale);
 #endif