From e076a00fde805c23c21bb68f66d0d1e7b77f3323 Mon Sep 17 00:00:00 2001 From: koennecke Date: Tue, 30 Jun 2009 06:43:23 +0000 Subject: [PATCH] - Fixed NB indexing - Fixed an issue with setbin generated time binnings --- histmem.c | 6 +- hmdata.c | 2 +- makefile | 1 + scan.c | 2 +- scriptcontext.c | 4 ++ singlenb.c | 54 ++++++++++++-- ubfour.c | 185 +++++++++++++++++++++++++++++++++++++++++++++++- ubfour.h | 29 +++++++- 8 files changed, 269 insertions(+), 14 deletions(-) diff --git a/histmem.c b/histmem.c index ed715f9e..cf10da9d 100644 --- a/histmem.c +++ b/histmem.c @@ -70,6 +70,8 @@ #include "event.h" #include "status.h" #include "site.h" +#define ABS(x) (x < 0 ? -(x) : (x)) + /* #define LOADDEBUG 1 */ @@ -1070,9 +1072,9 @@ static pDynString formatTOF(pHistMem self) * add the extra one */ if(iLength > 3){ - delta = timebin[iLength -2] - timebin[iLength -1]; + delta = ABS(timebin[iLength -2] - timebin[iLength -1]); } - snprintf(number,20," %12d", (int)(timebin[iLength -1] + delta)); + snprintf(number,20," %12d", (int)(timebin[iLength -1] + delta - delay )); DynStringConcat(result,number); return result; } else { diff --git a/hmdata.c b/hmdata.c index 138c3aa4..98492105 100644 --- a/hmdata.c +++ b/hmdata.c @@ -191,7 +191,7 @@ int setTimeBin(pHMdata self, int index, float value) return 0; } self->tofMode = 1; - if (index > self->nTimeChan) { + if (index >= self->nTimeChan) { self->nTimeChan = index + 1; return resizeBuffer(self); } diff --git a/makefile b/makefile index 03dfa3af..adc02ada 100644 --- a/makefile +++ b/makefile @@ -19,3 +19,4 @@ usage: @ ls -1 makefile_* | pr -t -o 4 +# DO NOT DELETE diff --git a/scan.c b/scan.c index f7624ec3..a64eb7b2 100644 --- a/scan.c +++ b/scan.c @@ -1764,7 +1764,7 @@ int ScanWrapper(SConnection * pCon, SicsInterp * pSics, void *pData, /* format them */ sprintf(pPtr, "scan.%s = ", ScanVarName(pVar)); for (i = 0; i < self->iNP; i++) { - sprintf(pItem, "{%12.3f} ", fData[i]); + sprintf(pItem, "{%12.4f} ", fData[i]); strcat(pPtr, pItem); } SCWrite(pCon, pPtr, eValue); diff --git a/scriptcontext.c b/scriptcontext.c index 2404e7b4..c73af0e0 100644 --- a/scriptcontext.c +++ b/scriptcontext.c @@ -967,6 +967,10 @@ void SctQueueNode(SctController * controller, Hdb * node, SctData *data; hdbCallback *cb; + /* + * The test for read below is questionable. If this becomes a problem + * take it out and tell Mark and Markus about the fact. + */ if (!FindHdbCallbackData(node, controller) && strcmp(action,"read") == 0) { cb = MakeHipadabaCallback(SctMainCallback, controller, NULL); assert(cb); diff --git a/singlenb.c b/singlenb.c index 454aedfc..cbecf592 100644 --- a/singlenb.c +++ b/singlenb.c @@ -213,6 +213,31 @@ static int getNBReflection(pSingleDiff self, char *id, reflection * r) } return 1; } +/*---------------------------------------------------------------------*/ +static int getNBNBReflection(pSingleDiff self, char *id, reflection * r) +{ + pSICSOBJ refList; + double hkl[3], angles[4], z1[3]; + double omnb, gamma, nu, stt, om, chi, phi; + + + refList = SXGetReflectionList(); + if (!GetRefIndexID(refList, id, hkl)) { + return 0; + } else { + r->h = hkl[0]; + r->k = hkl[1]; + r->l = hkl[2]; + GetRefAnglesID(refList, id, angles); + gamma = angles[0]; + omnb = angles[1]; + nu = angles[2]; + r->gamma = gamma; + r->omnb = omnb; + r->nu = nu; + } + return 1; +} /*-------------------------------------------------------------------*/ MATRIX calcNBUBFromTwo(pSingleDiff self, @@ -229,16 +254,16 @@ MATRIX calcNBUBFromTwo(pSingleDiff self, direct.beta = self->cell[4]; direct.gamma = self->cell[5]; - if (!getNBReflection(self, refid1, &r1)) { + if (!getNBNBReflection(self, refid1, &r1)) { *err = REFERR; return NULL; } - if (!getNBReflection(self, refid2, &r2)) { + if (!getNBNBReflection(self, refid2, &r2)) { *err = REFERR; return NULL; } - newUB = calcUBFromCellAndReflections(direct, r1, r2, err); + newUB = calcNBUBFromCellAndReflections(direct, r1, r2, err); return newUB; } @@ -250,20 +275,20 @@ MATRIX calcNBFromThree(pSingleDiff self, MATRIX newUB; reflection r1, r2, r3; - if (!getNBReflection(self, refid1, &r1)) { + if (!getNBNBReflection(self, refid1, &r1)) { *err = REFERR; return NULL; } - if (!getNBReflection(self, refid2, &r2)) { + if (!getNBNBReflection(self, refid2, &r2)) { *err = REFERR; return NULL; } - if (!getNBReflection(self, refid3, &r3)) { + if (!getNBNBReflection(self, refid3, &r3)) { *err = REFERR; return NULL; } - newUB = calcUBFromThreeReflections(r1, r2, r3, self->lambda, err); + newUB = calcNBUBFromThreeReflections(r1, r2, r3, self->lambda, err); return newUB; } @@ -271,12 +296,27 @@ MATRIX calcNBFromThree(pSingleDiff self, static int calcNBZ1(pSingleDiff self, char *refid, double z1[3]) { reflection r1; + pSICSOBJ refList; + double angles[4]; + double omnb, gamma, nu; + refList = SXGetReflectionList(); + if (!GetRefAnglesID(refList, refid, angles)) { + return 0; + } else { + gamma = angles[0]; + omnb = angles[1]; + nu = angles[2]; + z1FromNormalBeam(self->lambda, omnb, gamma, nu, z1); + return 1; + } + /* if (!getNBReflection(self, refid, &r1)) { return 0; } z1FromAngles(self->lambda, r1.s2t, r1.om, r1.chi, r1.chi, z1); return 1; + */ } /*--------------------------------------------------------------------*/ diff --git a/ubfour.c b/ubfour.c index ec571162..2e38f220 100644 --- a/ubfour.c +++ b/ubfour.c @@ -41,6 +41,35 @@ static MATRIX calcUVectorFromAngles(reflection r) vectorSet(u, 2, Cosd(om) * Sind(r.chi)); return u; } +/*-------------------------------------------------------------------------------------- + * The code for this was directly lifted from rafnb from ILL, routine UBH in rafnbb.f. + * I have yet to find a proper description of normal beam geometry calculations. + *--------------------------------------------------------------------------------------*/ +static MATRIX calcNBUVectorFromAngles(reflection r) +{ + MATRIX u; + double om, gamma, nu; + + u = makeVector(); + if (u == NULL) { + return NULL; + } + /* + VAB(1,I) = SN(1)*CS(2)*CS(3) + SN(2)*(1.-CS(1)*CS(3)) + VAB(2,I) = SN(1)*SN(2)*CS(3) - CS(2)*(1.-CS(1)*CS(3)) + VAB(3,I) = SN(3) + 1 = gamma, 2 = om, 3 = nu, CS = cos, SN = sin +*/ + gamma = r.gamma; + om = r.omnb; + nu = r.nu; + vectorSet(u, 0, + Sind(gamma)*Cosd(om)*Cosd(nu) + Sind(om)*(1. - Cosd(gamma)*Cosd(nu))); + vectorSet(u, 1, + Sind(gamma)*Sind(om)*Cosd(nu) - Cosd(om)*(1. - Cosd(gamma)*Cosd(nu))); + vectorSet(u, 2, Sind(nu)); + return u; +} /*--------------------------------------------------------------------------------------*/ static MATRIX reflectionToHC(reflection r, MATRIX B) @@ -153,6 +182,99 @@ MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1, return UB; } +/*---------------------------------------------------------------------------------------*/ +MATRIX calcNBUBFromCellAndReflections(lattice direct, reflection r1, + reflection r2, int *errCode) +{ + MATRIX B, HT, UT, U, UB, HTT; + MATRIX u1, u2, h1, h2; + double ud[3]; + int status; + reflection r; + + *errCode = 1; + + /* + calculate the B matrix and the HT matrix + */ + B = mat_creat(3, 3, ZERO_MATRIX); + status = calculateBMatrix(direct, B); + if (status < 0) { + *errCode = status; + return NULL; + } + h1 = reflectionToHC(r1, B); + h2 = reflectionToHC(r2, B); + if (h1 == NULL || h2 == NULL) { + *errCode = UBNOMEMORY; + return NULL; + } + HT = matFromTwoVectors(h1, h2); + if (HT == NULL) { + *errCode = UBNOMEMORY; + return NULL; + } + + /* + calculate U vectors and UT matrix + */ + u1 = calcNBUVectorFromAngles(r1); + u2 = calcNBUVectorFromAngles(r2); + if (u1 == NULL || u2 == NULL) { + *errCode = UBNOMEMORY; + return NULL; + } + UT = matFromTwoVectors(u1, u2); + if (UT == NULL) { + *errCode = UBNOMEMORY; + return NULL; + } + + /* + debugging output + printf("B-matrix\n"); + mat_dump(B); + printf("HT-matrix\n"); + mat_dump(HT); + printf("UT-matrix\n"); + mat_dump(UT); + */ + + /* + UT = U * HT + */ + HTT = mat_tran(HT); + if (HTT == NULL) { + *errCode = UBNOMEMORY; + return NULL; + } + U = mat_mul(UT, HTT); + if (U == NULL) { + *errCode = UBNOMEMORY; + return NULL; + } + UB = mat_mul(U, B); + if (UB == NULL) { + *errCode = UBNOMEMORY; + } + + /* + clean up + */ + killVector(h1); + killVector(h2); + mat_free(HT); + mat_free(HTT); + + killVector(u1); + killVector(u2); + mat_free(UT); + + mat_free(U); + mat_free(B); + + return UB; +} /*-----------------------------------------------------------------------------------*/ static void storeReflection(reflection r, double two_theta_obs, @@ -277,7 +399,7 @@ MATRIX calcUBFromThreeReflections(reflection r1, reflection r2, int *errCode) { MATRIX u1, u2, u3, HCHI, HI, HIINV, UB; - double det; + double det, z1[3]; if (lambda <= .1) { *errCode = INVALID_LAMBDA; @@ -292,7 +414,68 @@ MATRIX calcUBFromThreeReflections(reflection r1, reflection r2, uToScatteringVector(u1, r1.s2t / 2., lambda); uToScatteringVector(u2, r2.s2t / 2., lambda); uToScatteringVector(u3, r3.s2t / 2., lambda); + + HCHI = buildHCHIMatrix(u1, u2, u3); + HI = buildIndexMatrix(r1, r2, r3); + if (HCHI == NULL || HI == NULL) { + *errCode = UBNOMEMORY; + killVector(u1); + killVector(u2); + killVector(u3); + return NULL; + } + HIINV = mat_inv(HI); + if (HIINV == NULL) { + *errCode = NOTRIGHTHANDED; + killVector(u1); + killVector(u2); + killVector(u3); + mat_free(HI); + mat_free(HCHI); + return NULL; + } + UB = mat_mul(HCHI, HIINV); + + det = mat_det(UB); + if (det < .0) { + mat_free(UB); + UB = NULL; + *errCode = NOTRIGHTHANDED; + } + + killVector(u1); + killVector(u2); + killVector(u3); + mat_free(HI); + mat_free(HCHI); + mat_free(HIINV); + + return UB; +} +/*-----------------------------------------------------------------------------*/ +MATRIX calcNBUBFromThreeReflections(reflection r1, reflection r2, + reflection r3, double lambda, + int *errCode) +{ + MATRIX u1, u2, u3, HCHI, HI, HIINV, UB; + double det, z1[3]; + + if (lambda <= .1) { + *errCode = INVALID_LAMBDA; + return NULL; + } + + *errCode = 1; + + z1FromNormalBeam(lambda,r1.gamma,r1.omnb, r1.nu,z1); + u1 = vectorToMatrix(z1); + z1FromNormalBeam(lambda,r2.gamma,r2.omnb, r2.nu,z1); + u2 = vectorToMatrix(z1); + z1FromNormalBeam(lambda,r3.gamma,r3.omnb, r3.nu,z1); + u3 = vectorToMatrix(z1); + + HCHI = buildHCHIMatrix(u1, u2, u3); HI = buildIndexMatrix(r1, r2, r3); if (HCHI == NULL || HI == NULL) { diff --git a/ubfour.h b/ubfour.h index 90e595b2..b87b4558 100644 --- a/ubfour.h +++ b/ubfour.h @@ -33,7 +33,8 @@ */ typedef struct { double h, k, l; - double s2t, om, chi, phi; + double s2t, om, chi, phi; /* bisecting */ + double gamma, omnb, nu; /* normal beam */ } reflection; /** * calculate a UB matrix from cell constants and two reflections @@ -45,6 +46,17 @@ typedef struct { */ MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1, reflection r2, int *errCode); +/** + * calculate a UB matrix for normal beam geometry 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 calcNBUBFromCellAndReflections(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 @@ -58,6 +70,19 @@ MATRIX calcUBFromCellAndReflections(lattice direct, reflection r1, MATRIX calcUBFromThreeReflections(reflection r1, reflection r2, reflection r3, double lambda, int *errCode); +/** + * calculate the UB matrix from three reflections for normal beam geomtry. + * The three reflections must not be coplanar and must reflect a right handed set. + * @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 calcNBUBFromThreeReflections(reflection r1, reflection r2, + reflection r3, double lambda, + int *errCode); /** * a data structure holding an indexing suggestion */ @@ -92,7 +117,7 @@ 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 + * @param r The reflction for which to calculate * @return The length of the scattering vector */ double scatteringVectorLength(MATRIX B, reflection r);