From 5c30a7ea7b5e6d04829f21d9262be17869d224cd Mon Sep 17 00:00:00 2001 From: koennecke Date: Fri, 1 Apr 2005 13:48:25 +0000 Subject: [PATCH] - Added brute force indexing support to ubcalc - Added calculation of UB from 3 reflections to ubcalc - Added calculation of lattice constants from UB to ubcalc - Some fixes in stdscan in order to make the scripted scan work --- cell.c | 27 ++++++ cell.h | 9 +- fourlib.c | 2 +- fourlib.h | 6 ++ histsim.c | 24 +++-- interface.c | 30 +++++- stdscan.c | 37 +++++++- ubcalc.c | 267 +++++++++++++++++++++++++++++++++++++++++++++++----- ubfour.c | 238 +++++++++++++++++++++++++++++++++++++++++++++- ubfour.h | 42 ++++++++- 10 files changed, 637 insertions(+), 45 deletions(-) diff --git a/cell.c b/cell.c index 9779d656..68732f46 100644 --- a/cell.c +++ b/cell.c @@ -147,4 +147,31 @@ int calculateBMatrix(lattice direct, MATRIX B) { return 1; } +/*--------------------------------------------------------------------------*/ +int cellFromUB(MATRIX UB, plattice direct){ + MATRIX UBTRANS, GINV, G; + + UBTRANS = mat_tran(UB); + if(UBTRANS == NULL){ + return CELLNOMEMORY; + } + GINV = mat_mul(UBTRANS,UB); + if(GINV == NULL){ + mat_free(UBTRANS); + return CELLNOMEMORY; + } + G = mat_inv(GINV); + if(G == NULL){ + mat_free(UBTRANS); + mat_free(GINV); + return CELLNOMEMORY; + } + direct->a = sqrt(G[0][0]); + direct->b = sqrt(G[1][1]); + direct->c = sqrt(G[2][2]); + direct->alpha = Acosd(G[1][2]/(direct->b * direct->c)); + direct->beta = Acosd(G[2][0]/(direct->a * direct->c)); + direct->gamma = Acosd(G[0][1]/(direct->a * direct->c)); + return 1; +} diff --git a/cell.h b/cell.h index 8accfffc..4b1e042a 100644 --- a/cell.h +++ b/cell.h @@ -14,7 +14,7 @@ * error codes */ #define REC_NO_VOLUME -100 - +#define CELLNOMEMORY -101 /** * lattice parameters: either reciprocal or direct */ @@ -43,4 +43,11 @@ * @return 1 on success, an negative error code else */ int calculateBMatrix(lattice direct, MATRIX B); +/** + * calculate the cell constants from a UB matrix + * @param UB The input UB matrix. + * @param direct A pointer to a structure holding the new cell constants + * @return 1 on success, an error c ode < 0 on failure + */ +int cellFromUB(MATRIX UB, plattice direct); #endif diff --git a/fourlib.c b/fourlib.c index 69421c67..90314f73 100644 --- a/fourlib.c +++ b/fourlib.c @@ -213,7 +213,7 @@ static void turnEquatorial(MATRIX z1, double *chi, double *phi){ /*---------------------------------------------------------------------- calculate d-spacing and theta from z1 -----------------------------------------------------------------------*/ -static int calcTheta(double lambda, MATRIX z1, double *d, double *theta){ +int calcTheta(double lambda, MATRIX z1, double *d, double *theta){ double dstar, sintheta; dstar = sqrt(z1[0][0]*z1[0][0] + z1[1][0]*z1[1][0] + z1[2][0]*z1[2][0]); diff --git a/fourlib.h b/fourlib.h index 55cc7e5a..d7d3db70 100644 --- a/fourlib.h +++ b/fourlib.h @@ -162,6 +162,12 @@ void chimat(MATRIX chim, double chi); * phi rotation matrix for phi into phim */ void phimat(MATRIX phim, double phi); +/** + * calcTheta calculates theta for a z1 + * returns 1 on success, 0 else + */ +int calcTheta(double lambda, MATRIX z1, double *d, double *theta); + #endif diff --git a/histsim.c b/histsim.c index 20952150..d3bca918 100644 --- a/histsim.c +++ b/histsim.c @@ -50,12 +50,6 @@ #include "histsim.h" -#define TESTVAL 128 -/* - - define TESTVAL to set the simulated histogram to a fixed value for - testing -*/ static int iSet = 0; static HistInt iSetVal = 0; static HistMode eHistMode; @@ -86,6 +80,16 @@ free(self->pPriv); self->pPriv = NewSIMCounter("HistoSim",fFail); } + + /* + configured test value + */ + status = StringDictGet(pOption,"testval",pData,131); + if(status) + { + iSet = 1; + iSetVal = atoi(pData); + } return 1; } @@ -172,12 +176,7 @@ SCWrite(pCon,"ERROR: histogram out of range",eError); return 0; } -#ifdef TESTVAL - for(ii = iStart; ii < iEnd; ii++) - { - lData[ii-iStart] = TESTVAL; - } -#else + if(iSet == 1) { for(ii = iStart; ii < iEnd; ii++) @@ -192,7 +191,6 @@ lData[ii-iStart] = random(); } } -#endif return 1; } /*------------------------------------------------------------------------*/ diff --git a/interface.c b/interface.c index 4055e4b4..5c77abbf 100644 --- a/interface.c +++ b/interface.c @@ -42,7 +42,30 @@ #include "fortify.h" #include "sics.h" #include "motor.h" - +/*========================================================================= + Empty driveable interface functions + ==========================================================================*/ +static int EmptyHalt(void *self){ + return OKOK; +} +/*-----------------------------------------------------------------------*/ +static int EmptyLimits(void *self, float fVal, char *error, int errLen){ + return 1; +} +/*-----------------------------------------------------------------------*/ +static long EmptyValue(void *self, SConnection *pCon, float fVal){ + SCWrite(pCon,"WARNING: empty SetValue",eWarning); + return OKOK; +} +/*-----------------------------------------------------------------------*/ +static int EmptyStatus(void *self, SConnection *pCon){ + return HWIdle; +} +/*------------------------------------------------------------------------*/ +static float EmptyGet(void *self, SConnection *pCon){ + SCWrite(pCon,"WARNING: empty GetValue",eWarning); + return 555.55; +} /*-------------------------------------------------------------------------*/ pIDrivable CreateDrivableInterface(void) { @@ -55,6 +78,11 @@ } memset(pRes,0,sizeof(IDrivable)); pRes->ID = DRIVEID; + pRes->Halt = EmptyHalt; + pRes->CheckLimits = EmptyLimits; + pRes->SetValue = EmptyValue; + pRes->CheckStatus = EmptyStatus; + pRes->GetValue = EmptyGet; return pRes; } /*-------------------------------------------------------------------------*/ diff --git a/stdscan.c b/stdscan.c index 80197749..1be616a0 100644 --- a/stdscan.c +++ b/stdscan.c @@ -47,12 +47,15 @@ static char *fixExtension(char *filename) int iLen, iNum, iYear; char pNumText[10]; CommandList *pCom = NULL; + char simName[255]; /* make a simulated filename if in simulation mode */ - if(pServ->simMode) - return strdup("sim001001901.sim"); + if(pServ->simMode){ + snprintf(simName,255,"%s/tmp/sim001001901.sim",getenv("HOME")); + return strdup(simName); + } pRes = makeFilename(pSics,pCon); if(pRes == NULL) @@ -756,6 +759,10 @@ int ScriptWriteHeader(pScanData self){ } /*---------------------------------------------------------------------*/ int ScriptPrepareScan(pScanData self){ + /* configure counter */ + SetCounterMode((pCounter)self->pCounterData,self->iMode); + SetCounterPreset((pCounter)self->pCounterData, self->fPreset); + self->iCounts = 0; return StandardScriptInvoke(self,"prepare"); } /*-----------------------------------------------------------------------*/ @@ -764,7 +771,31 @@ int ScriptScanDrive(pScanData self, int iPoint){ } /*------------------------------------------------------------------------*/ int ScriptScanCount(pScanData self, int iPoint){ - return StandardScriptInvokeWithPoint(self,"count",iPoint); + pDynString command = NULL; + int status; + char pNumber[50]; + + command = GetStandardInvocation(self,"count"); + if(command == NULL){ + return 0; + } + snprintf(pNumber,49,"%d",iPoint); + DynStringConcat(command,pNumber); + if(self->iMode == eTimer){ + DynStringConcat(command," timer "); + } else { + DynStringConcat(command," monitor "); + } + snprintf(pNumber,49," %f ", self->fPreset); + DynStringConcat(command,pNumber); + + status = InterpExecute(self->pSics, self->pCon, + GetCharArray(command)); + DeleteDynString(command); + if(status != 1) { + return 0; + } + return status; } /*-------------------------------------------------------------------------*/ int ScriptScanCollect(pScanData self, int iPoint){ diff --git a/ubcalc.c b/ubcalc.c index 3d66027f..5757eb2f 100644 --- a/ubcalc.c +++ b/ubcalc.c @@ -20,9 +20,13 @@ /*---------------------------------------------------------------------*/ typedef struct { pObjectDescriptor pDes; + pHKL hkl; lattice direct; - reflection r1, r2; + reflection r1, r2, r3; MATRIX UB; + double allowedDeviation; + int indexSearchLimit; + int maxSuggestions; } UBCALC, *pUBCALC; /*----------------------------------------------------------------------*/ static void killUBCALC(void *pData){ @@ -53,10 +57,16 @@ static int SaveUBCalc(void *data, char *name, FILE *fd){ fprintf(fd,"%s ref2 %f %f %f %f %f %f %f\n",name, self->r2.h, self->r2.k, self->r2.l, self->r2.s2t, self->r2.om, self->r2.chi, self->r2.phi); + fprintf(fd,"%s ref3 %f %f %f %f %f %f %f\n",name, + self->r3.h, self->r3.k, self->r3.l, + self->r3.s2t, self->r3.om, self->r3.chi, self->r3.phi); + fprintf(fd,"%s difftheta %f\n", name, self->allowedDeviation); + fprintf(fd, "%s maxindex %d\n", name, self->indexSearchLimit); + fprintf(fd ,"%s maxlist %d\n", name, self->maxSuggestions); return 1; } /*---------------------------------------------------------------------*/ -static pUBCALC makeUBCALC(){ +static pUBCALC makeUBCALC(pHKL hkl){ pUBCALC pNew = NULL; pNew = (pUBCALC)malloc(sizeof(UBCALC)); @@ -70,6 +80,10 @@ static pUBCALC makeUBCALC(){ return NULL; } pNew->pDes->SaveStatus = SaveUBCalc; + pNew->hkl = hkl; + pNew->allowedDeviation = .3; + pNew->indexSearchLimit = 5; + pNew->maxSuggestions = 10; return pNew; } /*----------------------------------------------------------------------*/ @@ -77,17 +91,25 @@ int MakeUBCalc(SConnection *pCon, SicsInterp *pSics, void *pData, int argc, char *argv[]){ pUBCALC pNew = NULL; int status; + pHKL hkl = NULL; - pNew = makeUBCALC(); + if(argc < 3){ + SCWrite(pCon,"ERROR: missing argument to MakeUBCalc: MakeUBCalc name hklobject",eError); + return 0; + } + + hkl = FindCommandData(pSics,argv[2],"4-Circle-Calculus"); + if(hkl == NULL){ + SCWrite(pCon,"ERROR: HKL object not found or wrong type",eError); + return 0; + } + + pNew = makeUBCALC(hkl); if(pNew == NULL){ SCWrite(pCon,"ERROR: out of memory creating UBCALC",eError); return 0; } - if(argc > 1){ - status = AddCommand(pSics,argv[1],UBCalcWrapper,killUBCALC,pNew); - } else { - status = AddCommand(pSics,"ubcalc",UBCalcWrapper,killUBCALC,pNew); - } + status = AddCommand(pSics,argv[1],UBCalcWrapper,killUBCALC,pNew); if(status != 1){ SCWrite(pCon,"ERROR: failed to create duplicate UBCALC module",eError); } @@ -319,18 +341,12 @@ static void listReflection(SConnection *pCon, char *name, } /*---------------------------------------------------------------------*/ static int sendUBToHKL(SConnection *pCon, SicsInterp *pSics, - char *name, MATRIX UB){ - pHKL hkl = NULL; - char pBueffel[256]; + pHKL hkl, MATRIX UB){ 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; - } + assert(hkl != NULL); + for(i = 0; i < 3; i++){ ub[i] = UB[0][i]; ub[i+3] = UB[1][i]; @@ -340,6 +356,198 @@ static int sendUBToHKL(SConnection *pCon, SicsInterp *pSics, SCSendOK(pCon); return 1; } +/*---------------------------------------------------------------------*/ +static int setUBCalcParameters(pUBCALC self, SConnection *pCon, char *name, char *value){ + + if(strcmp(name,"difftheta") == 0){ + if(!SCMatchRights(pCon,usUser)){ + return 0; + } + self->allowedDeviation = atof(value); + SCparChange(pCon); + SCSendOK(pCon); + return 1; + } else if(strcmp(name,"maxindex") == 0){ + if(!SCMatchRights(pCon,usUser)){ + return 0; + } + self->indexSearchLimit = atoi(value); + SCparChange(pCon); + SCSendOK(pCon); + return 1; + } else if(strcmp(name,"maxlist") == 0){ + if(!SCMatchRights(pCon,usUser)){ + return 0; + } + self->maxSuggestions = atoi(value); + SCparChange(pCon); + SCSendOK(pCon); + return 1; + } else { + SCWrite(pCon,"ERROR: subcommand not recognized",eError); + return 0; + } +} +/*---------------------------------------------------------------------*/ +static int getUBCalcParameters(pUBCALC self, SConnection *pCon, char *name){ + char pBueffel[255]; + int ret = 0; + + if(strcmp(name,"difftheta") == 0){ + snprintf(pBueffel,255,"ubcalc.difftheta = %f",self->allowedDeviation); + ret = 1; + } else if(strcmp(name,"maxindex") == 0){ + snprintf(pBueffel,255,"ubcalc.maxindex = %d", self->indexSearchLimit); + ret = 1; + } else if(strcmp(name,"maxlist") == 0){ + snprintf(pBueffel,255,"ubcalc.maxindex = %d", self->maxSuggestions); + ret = 1; + } else { + snprintf(pBueffel,255,"ERROR: subcommand not known"); + } + SCWrite(pCon,pBueffel,eValue); + return ret; +} +/*---------------------------------------------------------------------*/ +static void listPar(pUBCALC self, char *name, SConnection *pCon){ + char pBueffel[255]; + + snprintf(pBueffel,255,"%s.difftheta = %f", name, self->allowedDeviation); + SCWrite(pCon,pBueffel,eValue); + snprintf(pBueffel,255,"%s.maxindex = %d", name, self->indexSearchLimit); + SCWrite(pCon,pBueffel,eValue); + snprintf(pBueffel,255,"%s.maxlist = %d", name, self->maxSuggestions); + SCWrite(pCon,pBueffel,eValue); + +} +/*---------------------------------------------------------------------*/ +static int findIndex(pUBCALC self, SConnection *pCon, SicsInterp *pSics, + int argc, char *argv[]){ + float two_theta, lambda; + pMotor pMot = NULL; + int status, numRef, i; + refIndex *index = NULL; + Tcl_DString list; + char pLine[255]; + + if(argc > 2){ + two_theta = atof(argv[2]); + } else { + pMot = FindMotor(pSics,"stt"); + if(pMot == NULL){ + SCWrite(pCon,"ERROR: cannot find stt motor",eError); + return 0; + } + MotorGetSoftPosition(pMot,pCon,&two_theta); + } + + GetLambda(self->hkl,&lambda); + + numRef = self->maxSuggestions; + index = (refIndex *)malloc(numRef*sizeof(refIndex)); + if(index == NULL){ + SCWrite(pCon,"ERROR: out of memory allocating index list",eError); + return 0; + } + memset(index,0,numRef*sizeof(refIndex)); + + status = searchIndex(self->direct,lambda, two_theta,self->allowedDeviation, + self->indexSearchLimit, index, numRef); + + if(status < 0){ + if(status == UBNOMEMORY){ + SCWrite(pCon,"ERROR: out of memory searching indices",eError); + return 0; + } else if(status == REC_NO_VOLUME){ + SCWrite(pCon,"ERROR: bad cell parameters",eError); + return 0; + } else { + SCWrite(pCon,"ERROR: unknown error code",eError); + return 0; + } + } + + if(status < numRef) { + numRef = status; + } + Tcl_DStringInit(&list); + for(i = 0; i < numRef; i++){ + snprintf(pLine,255," %3d %3d %3d %8.2f %8.2f %8.2f\r\n", + (int)index[i].h, (int)index[i].k, (int)index[i].l, + two_theta, index[i].t2calc,index[i].t2diff); + Tcl_DStringAppend(&list,pLine,-1); + } + if(numRef == 0){ + Tcl_DStringAppend(&list,"No suitable reflections found",-1); + } + SCWrite(pCon,Tcl_DStringValue(&list),eValue); + + Tcl_DStringFree(&list); + free(index); + return 1; +} +/*---------------------------------------------------------------------*/ +static int calcUB3Ref(pUBCALC self, SConnection *pCon){ + float lambda; + MATRIX newUB = NULL; + int errCode = 1; + char pBueffel[256]; + + GetLambda(self->hkl,&lambda); + + newUB = calcUBFromThreeReflections(self->r1, self->r2, self->r3, + lambda, &errCode); + if(newUB == NULL){ + switch(errCode){ + case UBNOMEMORY: + SCWrite(pCon,"ERROR: out of memory calculating UB",eError); + break; + case INVALID_LAMBDA: + SCWrite(pCon,"ERROR: bad value for wavelength",eError); + break; + case NOTRIGHTHANDED: + SCWrite(pCon,"ERROR: reflections are coplanar or do not form a right handed system", + eError); + break; + default: + SCWrite(pCon,"ERROR: unknown error code from UB calculation",eError); + break; + } + return 0; + } else { + if(self->UB != NULL){ + mat_free(self->UB); + } + self->UB = newUB; + SCSendOK(pCon); + } + return 1; +} +/*--------------------------------------------------------------------*/ +static int cellFromUBWrapper(pUBCALC self, SConnection *pCon){ + int status; + char pBueffel[256]; + lattice direct; + + status = cellFromUB(self->UB,&direct); + if(status < 0){ + switch(status){ + case CELLNOMEMORY: + SCWrite(pCon,"ERROR: out of memory while calculating cell",eError); + break; + default: + SCWrite(pCon,"ERROR: unknown error calculating cell",eError); + break; + } + return 0; + } else { + snprintf(pBueffel,255," %8.4f %8.4f %8.4f %6.2f %6.2f %6.2f", + direct.a, direct.b, direct.c, + direct.alpha, direct.beta, direct.gamma); + SCWrite(pCon,pBueffel,eValue); + } + return 1; +} /*----------------------------------------------------------------------*/ int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData, int argc, char *argv[]){ @@ -358,11 +566,17 @@ int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData, return readReflection(pCon,pSics, &self->r1,argc,argv); } else if(strcmp(argv[1],"ref2") ==0){ return readReflection(pCon,pSics, &self->r2,argc,argv); + } else if(strcmp(argv[1],"ref3") == 0){ + return readReflection(pCon,pSics, &self->r3,argc,argv); } else if(strcmp(argv[1],"listub") == 0){ listUB(pCon,self->UB); return 1; } else if(strcmp(argv[1],"calcub") == 0){ return calcUB(self,pCon); + } else if(strcmp(argv[1],"calcub3ref") == 0){ + return calcUB3Ref(self,pCon); + } else if(strcmp(argv[1],"cellfromub") == 0){ + return cellFromUBWrapper(self,pCon); } else if(strcmp(argv[1],"listcell") == 0){ listCell(pCon,argv[0],self->direct); return 1; @@ -372,18 +586,27 @@ int UBCalcWrapper(SConnection *pCon, SicsInterp *pSics, void *pData, } else if(strcmp(argv[1],"listref2") == 0){ listReflection(pCon,argv[0],"ref2",self->r2); return 1; + } else if(strcmp(argv[1],"listref3") == 0){ + listReflection(pCon,argv[0],"ref3",self->r3); + return 1; } else if(strcmp(argv[1],"list") == 0){ listCell(pCon,argv[0],self->direct); listReflection(pCon,argv[0],"ref1",self->r1); listReflection(pCon,argv[0],"ref2",self->r2); + listReflection(pCon,argv[0],"ref3",self->r3); listUB(pCon,self->UB); + listPar(self,argv[0],pCon); return 1; - } else if(strcmp(argv[1],"sendto") == 0){ - if(argc < 3){ - SCWrite(pCon,"ERROR: Nothing to send UB too",eError); - return 0; + } else if(strcmp(argv[1],"activate") == 0){ + return sendUBToHKL(pCon,pSics,self->hkl,self->UB); + } else if(strcmp(argv[1],"index") == 0){ + return findIndex(self,pCon,pSics,argc, argv); + } else { + if(argc > 2){ + return setUBCalcParameters(self,pCon,argv[1],argv[2]); + } else { + return getUBCalcParameters(self,pCon,argv[1]); } - return sendUBToHKL(pCon,pSics,argv[2],self->UB); } return 1; } diff --git a/ubfour.c b/ubfour.c index ac5ec6e1..dd52be23 100644 --- a/ubfour.c +++ b/ubfour.c @@ -1,11 +1,12 @@ /** * This is a library for calculating UB matrices for four circle diffraction. - * The algorithm and setting definitions is from: + * The algorithm and settings definition is from: * * Busing & Levy, Acta Cryst. (1967), 22, 457ff * * Implemented: * - UB from cell cell constants and two reflections. + * - Brute force index search * * Mark Koennecke, March 2005 */ @@ -13,6 +14,9 @@ #include #include "vector.h" #include "trigd.h" +#include "fourlib.h" +#include "lld.h" +#define ABS(x) (x < 0 ? -(x) : (x)) /*--------------------------------------------------------------------------------------*/ static MATRIX calcUVectorFromAngles(reflection r){ MATRIX u; @@ -50,6 +54,7 @@ MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1, MATRIX u1, u2, h1, h2; double ud[3]; int status; + reflection r; *errCode = 1; @@ -134,3 +139,234 @@ MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1, return UB; } +/*-----------------------------------------------------------------------------------*/ +static void storeReflection(reflection r, double two_theta_obs, + double two_theta_calc, int list){ + refIndex ri, test; + int count = 0, status, pos = 0; + + ri.h = r.h; + ri.k = r.k; + ri.l = r.l; + ri.t2obs = two_theta_obs; + ri.t2calc = two_theta_calc; + ri.t2diff = ABS(two_theta_obs - two_theta_calc); + + /* + locate the last entry bigger then us + */ + status = LLDnodePtr2First(list); + while(status == 1){ + LLDnodeDataTo(list,&test); + count++; + if(test.t2diff == ri.t2diff){ + LLDnodeDataFrom(list,&ri); + return; + } + if(test.t2diff > ri.t2diff){ + break; + } + status = LLDnodePtr2Next(list); + } + /* + special case: empty list + */ + if(count == 0){ + LLDnodeAppendFrom(list,&ri); + return; + } + /* + special case: append after last + */ + LLDnodePtr2Last(list); + LLDnodeDataTo(list,&test); + if(ri.t2diff > test.t2diff){ + LLDnodeAppendFrom(list,&ri); + return; + } + + status = LLDnodePtr2First(list); + pos = 0; + while(status == 1){ + LLDnodeDataTo(list,&test); + pos++; + if(pos == count){ + LLDnodeInsertFrom(list,&ri); + return; + } + status = LLDnodePtr2Next(list); + } +} +/*---------------------------------------------------------------------------- + u_transform(i) = u(i)*(2*sin(theta)/lambda) + -----------------------------------------------------------------------------*/ +static void uToScatteringVector(MATRIX u, double theta, double lambda){ + double scale; + int i; + + scale = (2. * Sind(theta))/lambda; + for(i = 0; i < 3; i++){ + u[i][0] *= scale; + } +} +/*----------------------------------------------------------------------------*/ +static MATRIX buildHCHIMatrix(MATRIX u1, MATRIX u2, MATRIX u3){ + int i; + MATRIX HCHI; + + HCHI = mat_creat(3,3,ZERO_MATRIX); + if(HCHI == NULL){ + return NULL; + } + for(i = 0; i < 3; i++){ + HCHI[i][0] = u1[i][0]; + HCHI[i][1] = u2[i][0]; + HCHI[i][2] = u3[i][0]; + } + return HCHI; +} +/*----------------------------------------------------------------------------*/ +static MATRIX buildIndexMatrix(reflection r1, reflection r2, reflection r3){ + MATRIX HI; + int i; + + HI = mat_creat(3,3,ZERO_MATRIX); + if(HI == NULL){ + return NULL; + } + HI[0][0] = r1.h; + HI[1][0] = r1.k; + HI[2][0] = r1.l; + + HI[0][1] = r2.h; + HI[1][1] = r2.k; + HI[2][1] = r2.l; + + HI[0][2] = r3.h; + HI[1][2] = r3.k; + HI[2][2] = r3.l; + + return HI; +} +/*-----------------------------------------------------------------------------*/ +MATRIX calcUBFromThreeReflections(reflection r1, reflection r2, reflection r3, + double lambda, int *errCode){ + MATRIX u1, u2, u3, HCHI, HI, HIINV, UB; + double det; + + if(lambda <= .1){ + *errCode = INVALID_LAMBDA; + return NULL; + } + + *errCode = 1; + + u1 = calcUVectorFromAngles(r1); + u2 = calcUVectorFromAngles(r2); + u3 = calcUVectorFromAngles(r3); + uToScatteringVector(u1,r1.s2t/2.,lambda); + uToScatteringVector(u2,r2.s2t/2.,lambda); + uToScatteringVector(u3,r3.s2t/2.,lambda); + + HCHI = buildHCHIMatrix(u1,u2,u3); + HI = buildIndexMatrix(r1,r2,r3); + if(HCHI == NULL || HI == NULL){ + *errCode = UBNOMEMORY; + killVector(u1); + killVector(u2); + killVector(u3); + return NULL; + } + + HIINV = mat_inv(HI); + if(HIINV == NULL){ + *errCode = UBNOMEMORY; + killVector(u1); + killVector(u2); + killVector(u3); + mat_free(HI); + mat_free(HCHI); + return NULL; + } + UB = mat_mul(HCHI,HIINV); + + det = mat_det(UB); + if(det < .0){ + mat_free(UB); + UB = NULL; + *errCode = NOTRIGHTHANDED; + } + + killVector(u1); + killVector(u2); + killVector(u3); + mat_free(HI); + mat_free(HCHI); + mat_free(HIINV); + + return UB; +} +/*---------------------------------------------------------------------------------*/ +static int copyReflections(int list, refIndex index[], int maxIndex){ + int count = 0, status; + refIndex ri; + + status = LLDnodePtr2First(list); + while(status == 1 && count < maxIndex){ + LLDnodeDataTo(list,&ri); + index[count] = ri; + status = LLDnodePtr2Next(list); + count++; + } + return count; +} +/*----------------------------------------------------------------------------------- + - matching reflections will be entered in to a list in a sorted way. This list + is copied into the index array. + - There is some waste here in allocating and deallocating the HC vector in the + inner loop. I'am to lazy to resolve this.... May be I'am spared..... + -----------------------------------------------------------------------------------*/ +int searchIndex(lattice direct, double lambda, double two_theta, double max_deviation, + int limit, refIndex index[], int maxIndex){ + int status, i, j, k, list; + MATRIX B, HC; + double theta, d; + reflection r; + + B = mat_creat(3,3,UNIT_MATRIX); + if(B == NULL) { + return UBNOMEMORY; + } + + status = calculateBMatrix(direct,B); + if(status < 0) { + return status; + } + + list = LLDcreate(sizeof(refIndex)); + if(list <0) { + return UBNOMEMORY; + } + + for(i = -limit; i < limit; i++){ + r.h = (double)i; + for(j = -limit; j < limit; j++){ + r.k = (double)j; + for(k = -limit; k < limit; k++){ + r.l = (double)k; + HC = reflectionToHC(r,B); + status = calcTheta(lambda, HC, &d, &theta); + if(status == 1){ + if(ABS(two_theta - 2. * theta) <= max_deviation){ + storeReflection(r,two_theta, theta * 2., list); + } + } + killVector(HC); + } + } + } + mat_free(B); + status = copyReflections(list,index,maxIndex); + LLDdelete(list); + return status; +} diff --git a/ubfour.h b/ubfour.h index 84134b68..9a0da3ea 100644 --- a/ubfour.h +++ b/ubfour.h @@ -3,8 +3,9 @@ * The algorithm and setting definitions is from: * Busing & Levy, Acta Cryst. (1967), 22, 457ff * - * This initial version only supports the calculation of the UB matrix from - * cell constants and two reflections. + * Implemented: + * - UB from cell cell constants and two reflections. + * - Brute force index search * * Mark Koennecke, march 2005 */ @@ -19,7 +20,8 @@ * error codes: also see cell.h */ #define UBNOMEMORY -200 - +#define INVALID_LAMBDA -201 +#define NOTRIGHTHANDED -202 /** * a reflection data structure holding the indices h,k,l and * the magic four circle angles two_theta, om, chi and phi. @@ -37,4 +39,38 @@ typedef struct{ * @return The resulting UB matrix or NULL on errors */ MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1, reflection r2, int *errCode); +/** + * calculate the UB matrix from three reflections. The three reflections must not be + * coplanar and must reflect a right handed + * @param r1 The first reflection + * @param r2 The second reflection + * @param r3 The third reflection. + * @param errCode a code indictang errors which happened. + * @return A UB matrix on success or NULL on errors. Then errcode will indicate + * the type of teh error. + */ +MATRIX calcUBFromThreeReflections(reflection r1, reflection r2, reflection r3, + double lambda, int *errCode); +/** + * a data structure holding an indexing suggestion + */ +typedef struct { + double h, k, l; /* suggested index */ + double t2obs, t2calc, t2diff; /* 2 theta observed and calculated and the difference */ +}refIndex, *prefIndex; +/** + * search for possible indexes matching the two theta value given. This is a brute force search + * @param direct The lattice constants of the crystal + * @param lambda The wavelength used + * @param two_theta The two theta value of the reflection + * @param max_deviation maximum allowed diviation of calculated two thetas from tw-Theta given + * @param limit Index limit for the search + * @param index Preallocated array for storing index suggestions. The procedure will + * sort the entries in this array according to the difference in two theta + * @param maxIndex The number of entries allowed for index. + * @return The number of indexes in index or a negative error code when an error + * occurs. + */ +int searchIndex(lattice direct, double lambda, double two_theta, double max_deviation, + int limit, refIndex index[], int maxIndex); #endif