From 4f069341f5fce54275571fe5869dd91ed2fdb516 Mon Sep 17 00:00:00 2001 From: koennecke Date: Fri, 20 Oct 2006 14:54:39 +0000 Subject: [PATCH] - added makeauxub to tasub - fixes for MARS - extended maximize to honour maxpts and the in360 flag - added regression histmem driver --- counter.c | 1 + fitcenter.c | 7 ++ histmem.c | 22 +++- histregress.c | 273 +++++++++++++++++++++++++++++++++++++++++++++++++ hkl.c | 4 +- make_gen | 2 +- maximize.c | 73 ++++++++++--- multicounter.c | 49 +++++++++ optimise.c | 2 +- simcter.c | 2 +- tasdrive.c | 2 +- tasub.c | 111 +++++++++++++++++++- tasublib.c | 100 +++++++++++++++--- tasublib.h | 13 +++ 14 files changed, 623 insertions(+), 38 deletions(-) create mode 100644 histregress.c diff --git a/counter.c b/counter.c index 1af90206..7f585787 100644 --- a/counter.c +++ b/counter.c @@ -914,6 +914,7 @@ SCSendOK(pCon); return 1; case 11: /* status */ + self->pCountInt->TransferData(self,pCon); if(GetCounterMode(self) == ePreset) { sprintf(pBueffel,"%s.CountStatus = %d %d Beam: %ld E6", diff --git a/fitcenter.c b/fitcenter.c index dbfeb30d..c39bfff0 100644 --- a/fitcenter.c +++ b/fitcenter.c @@ -467,6 +467,13 @@ sprintf(pBueffel,"%f", self->fCenter); SCWrite(pCon,pBueffel,eValue); return 1; + } + if(strcmp(argv[1],"data") == 0) + { + snprintf(pBueffel,255,"%f,%f,%ld", + self->fCenter, self->FWHM, self->lPeak); + SCWrite(pCon,pBueffel,eValue); + return 1; } } diff --git a/histmem.c b/histmem.c index d29de6a1..07de915f 100644 --- a/histmem.c +++ b/histmem.c @@ -68,6 +68,10 @@ /* #define LOADDEBUG 1 */ +/* + * from histregress.c + */ +extern pHistDriver CreateRegressHM(pStringDict pOpt); /*------------------------------------------------------------------------*/ static int HistHalt(void *pData) { @@ -452,8 +456,12 @@ } else if(strcmp(driver,"mcstas") == 0) { - pNew->pDriv = NewMcStasHM(pOption); - } + pNew->pDriv = NewMcStasHM(pOption); + } + else if(strcmp(driver,"regress") == 0) + { + pNew->pDriv = CreateRegressHM(pOption); + } else { site = getSite(); @@ -1264,6 +1272,16 @@ static int checkHMEnd(pHistMem self, char *text){ SCSendOK(pCon); return 1; } + /* pause */ + else if(strcmp(argv[1],"pause") == 0) + { if(!SCMatchRights(pCon,usUser)) + { + return 0; + } + self->pDriv->Pause(self->pDriv,pCon); + SCSendOK(pCon); + return 1; + } /* normal counting*/ else if(strcmp(argv[1],"count") == 0) { diff --git a/histregress.c b/histregress.c new file mode 100644 index 00000000..8003e27f --- /dev/null +++ b/histregress.c @@ -0,0 +1,273 @@ +/*---------------------------------------------------------------------------- + + H I S T S I M + + A simulated histogram memory for regression tests. + + All the counting error stuff is redirected to a regression counter; see + documentation there. This just adds data handling. + + copyright: see file COPYRIGHT + + Mark Koennecke, October 2006 + ----------------------------------------------------------------------------*/ +#include +#include +#include +#include +#include "fortify.h" +#include "sics.h" +#include "countdriv.h" +#include "counter.h" +#include "stringdict.h" +#include "HistMem.h" +#include "HistDriv.i" +#include "histsim.h" + +static int iSet = 0; +static HistInt iSetVal = 0; +static HistMode eHistMode; +/*--------------------------------------------------------------------------*/ +static int RegressConfig(pHistDriver self, SConnection *pCon, + pStringDict pOption, SicsInterp *pSics) +{ + int i, iLength = 1, status; + char pData[132]; + float fFail; + pCounterDriver count; + + count = (pCounterDriver)self->pPriv; + + if(eHistMode == eHTOF) + { + for(i = 0; i < self->data->rank; i++) + { + iLength *= self->data->iDim[i]; + } + iLength *= self->data->nTimeChan; + } + + /* + deal with error settings + */ + status = StringDictGet(pOption,"errortype",pData,131); + if(status) + { + + fFail = atof(pData); + count->Set(count,"errortype",1,fFail); + } + status = StringDictGet(pOption,"recover",pData,131); + if(status) + { + + fFail = atof(pData); + count->Set(count,"recover",1,fFail); + } + status = StringDictGet(pOption,"finish",pData,131); + if(status) + { + + fFail = atof(pData); + count->Set(count,"finish",1,fFail); + } + + /* + configured test value + */ + status = StringDictGet(pOption,"testval",pData,131); + if(status) + { + iSet = 1; + iSetVal = atoi(pData); + } + + return 1; + } +/*-------------------------------------------------------------------------*/ + static int RegressStart(pHistDriver self, SConnection *pCon) + { + pCounterDriver pDriv; + + pDriv = (pCounterDriver)self->pPriv; + pDriv->fPreset = self->fCountPreset; + pDriv->eMode = self->eCount; + return pDriv->Start(pDriv); + } +/*-------------------------------------------------------------------------*/ + static int RegressPause(pHistDriver self, SConnection *pCon) + { + pCounterDriver pDriv; + + pDriv = (pCounterDriver)self->pPriv; + pDriv->fPreset = self->fCountPreset; + pDriv->eMode = self->eCount; + return pDriv->Pause(pDriv); + } +/*------------------------------------------------------------------------*/ + static int RegressContinue(pHistDriver self, SConnection *pCon) + { + pCounterDriver pDriv; + + pDriv = (pCounterDriver)self->pPriv; + pDriv->fPreset = self->fCountPreset; + pDriv->eMode = self->eCount; + return pDriv->Continue(pDriv); + } +/*-------------------------------------------------------------------------*/ + static int RegressHalt(pHistDriver self) + { + pCounterDriver pDriv; + + pDriv = (pCounterDriver)self->pPriv; + return pDriv->Halt(pDriv); + } +/*-------------------------------------------------------------------------*/ + static int RegressGetCountStatus(pHistDriver self, SConnection *pCon) + { + pCounterDriver pDriv; + float fControl; + + pDriv = (pCounterDriver)self->pPriv; + return pDriv->GetStatus(pDriv,&fControl); + } +/*-------------------------------------------------------------------------*/ + static int RegressGetError(pHistDriver self, int *iCode, char *pError, int iLen) + { + pCounterDriver pDriv; + + pDriv = (pCounterDriver)self->pPriv; + return pDriv->GetError(pDriv, iCode,pError,iLen); + } +/*-------------------------------------------------------------------------*/ + static int RegressTryAndFixIt(pHistDriver self, int iCode) + { + pCounterDriver pDriv; + + pDriv = (pCounterDriver)self->pPriv; + return pDriv->TryAndFixIt(pDriv, iCode); + } +/*--------------------------------------------------------------------------*/ + static int RegressGetData(pHistDriver self, SConnection *pCon) + { + pCounterDriver pDriv; + + pDriv = (pCounterDriver)self->pPriv; + + return pDriv->ReadValues(pDriv); + } +/*--------------------------------------------------------------------------*/ + static int RegressGetHistogram(pHistDriver self, SConnection *pCon, + int i, int iStart, int iEnd, HistInt *lData) + { + int ii; + + if(i < 0) + { + SCWrite(pCon,"ERROR: histogram out of range",eError); + return 0; + } + + if(iSet == 1) + { + for(ii = iStart; ii < iEnd; ii++) + { + lData[ii-iStart] = iSetVal; + } + } + else + { + for(ii = iStart; ii < iEnd; ii++) + { + lData[ii-iStart] = random(); + } + } + return 1; + } +/*------------------------------------------------------------------------*/ + static int RegressSetHistogram(pHistDriver self, SConnection *pCon, + int i, int iStart, int iEnd, HistInt *lData) + { + iSet = 1; + iSetVal = lData[0]; + return 1; + } + +/*-------------------------------------------------------------------------*/ + static int RegressPreset(pHistDriver self, SConnection *pCon, HistInt iVal) + { + iSet = 1; + iSetVal = iVal; + return 1; + } +/*------------------------------------------------------------------------*/ + static int RegressFreePrivate(pHistDriver self) + { + pCounterDriver pDriv; + + pDriv = (pCounterDriver)self->pPriv; + DeleteCounterDriver(pDriv); + return 1; + } +/*------------------------------------------------------------------------*/ + static long RegressGetMonitor(pHistDriver self, int i, SConnection *pCon) + { + pCounterDriver pDriv; + long lVal; + + pDriv = (pCounterDriver)self->pPriv; + return pDriv->lCounts[i]; + + } +/*------------------------------------------------------------------------*/ + static float RegressGetTime(pHistDriver self, SConnection *pCon) + { + pCounterDriver pDriv; + long lVal; + + pDriv = (pCounterDriver)self->pPriv; + return pDriv->fTime; + + } +/*-------------------------------------------------------------------------*/ + pHistDriver CreateRegressHM(pStringDict pOpt) + { + pHistDriver pNew = NULL; + + /* create the general driver */ + pNew = CreateHistDriver(pOpt); + if(!pNew) + { + return NULL; + } + + /* put a Regresscounter in */ + pNew->pPriv = (void *)NewRegressCounter("HistoRegress"); + if(!pNew->pPriv) + { + DeleteHistDriver(pNew); + return NULL; + } + + /* configure all those functions */ + pNew->Configure = RegressConfig; + pNew->Start = RegressStart; + pNew->Halt = RegressHalt; + pNew->GetCountStatus = RegressGetCountStatus; + pNew->GetError = RegressGetError; + pNew->TryAndFixIt = RegressTryAndFixIt; + pNew->GetData = RegressGetData; + pNew->GetHistogram = RegressGetHistogram; + pNew->SetHistogram = RegressSetHistogram; + pNew->GetMonitor = RegressGetMonitor; + pNew->GetTime = RegressGetTime; + pNew->Preset = RegressPreset; + pNew->FreePrivate = RegressFreePrivate; + pNew->Pause = RegressPause; + pNew->Continue = RegressContinue; + StringDictAddPair(pOpt,"errortype","0"); + StringDictAddPair(pOpt,"recover","1"); + StringDictAddPair(pOpt,"testval","0"); + + return pNew; + } diff --git a/hkl.c b/hkl.c index 64c54a6b..d40176d9 100644 --- a/hkl.c +++ b/hkl.c @@ -794,7 +794,7 @@ static int calculateNormalBeam(MATRIX z1, pHKL self, SConnection *pCon, status = z1mToNormalBeam(self->fLambda, z3, &gamma, &omnb, &nu); - omnb += 180.; + /* omnb += 180.; */ mat_free(z3); if(status != 1) { @@ -803,7 +803,7 @@ static int calculateNormalBeam(MATRIX z1, pHKL self, SConnection *pCon, if(checkNormalBeam(omnb, &gamma, nu,fSet,pCon,self)){ return 1; } else { - if(checkNormalBeam(360. - omnb, &gamma, nu, fSet,pCon,self)){ + if(checkNormalBeam(omnb + 360., &gamma, nu, fSet,pCon,self)){ return 1; } else { return 0; diff --git a/make_gen b/make_gen index 6af7f61f..e50953e1 100644 --- a/make_gen +++ b/make_gen @@ -31,7 +31,7 @@ SOBJ = network.o ifile.o conman.o SCinter.o splitter.o passwd.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 cone.o hipadaba.o sicshipadaba.o statistics.o \ - moregress.o hdbcommand.o multicounter.o regresscter.o + moregress.o hdbcommand.o multicounter.o regresscter.o histregress.o MOTOROBJ = motor.o simdriv.o COUNTEROBJ = countdriv.o simcter.o counter.o diff --git a/maximize.c b/maximize.c index 4857e45f..be65ddf9 100644 --- a/maximize.c +++ b/maximize.c @@ -55,6 +55,7 @@ pObjectDescriptor pDes; pCounter pCount; int i360; + int maxpts; }Maxxii; /*----------------------------------------------------------------------- @@ -167,10 +168,10 @@ /* search to the left until out of space or lCts < lMax/2. */ SCWrite(pCon,"Searching for low angle boundary..",eWarning); - for(i = MAXPTS/2; i >= 0; i--) + for(i = self->maxpts/2; i >= 0; i--) { /* drive motor */ - fPos = fStart - (MAXPTS/2 - i)*fStep; + fPos = fStart - (self->maxpts/2 - i)*fStep; fPos = in360(self,fPos); if(maxDrive(pVar,pVarName,fPos,pCon) != 1) { @@ -218,7 +219,7 @@ goto start; } /* no peak found or normal peak: continue at other side */ - if( (i < 1) || (y[MAXPTS/2] > lMax/2) ) + if( (i < 1) || (y[self->maxpts/2] > lMax/2) ) { iSkip = 0; } @@ -227,7 +228,7 @@ /* next case: all of the peak in measured half: find max value and skip the right half */ - for(i = MAXPTS/2; i > 0; i--) + for(i = self->maxpts/2; i > 0; i--) { if(y[i] > lMax/2) { @@ -246,10 +247,10 @@ lMin = 100000; lMax = -100000; SCWrite(pCon,"Searching for high angle boundary..",eWarning); - for(i = MAXPTS/2; i < MAXPTS; i++) + for(i = self->maxpts/2; i < self->maxpts; i++) { /* drive motor */ - fPos = fStart + (i - MAXPTS/2) * fStep; + fPos = fStart + (i - self->maxpts/2) * fStep; fPos = in360(self,fPos); if(maxDrive(pVar,pVarName,fPos,pCon) != 1) { @@ -292,18 +293,18 @@ iTop = i; iTop++; /* first case: peak is at high angle side */ - if( (i > MAXPTS-2) && (lMax*0.5 > lMin) ) + if( (i > self->maxpts-2) && (lMax*0.5 > lMin) ) { goto start; } /* second case: no peak */ - if( (iTop > MAXPTS-2) ) + if( (iTop > self->maxpts-2) ) { SCWrite(pCon,"ERROR: no peak found!",eError); return 0; } /* third case: normal peak */ - if(y[MAXPTS/2] >= 0.5*lMax) + if(y[self->maxpts/2] >= 0.5*lMax) { iTop--; } @@ -313,7 +314,7 @@ */ else { - for(i = MAXPTS/2; i < MAXPTS; i++) + for(i = self->maxpts/2; i < self->maxpts; i++) { if(y[i] > lMax/2) { @@ -325,7 +326,7 @@ } } /* end of iSkip */ - if( (iBot < 2) || (iTop > MAXPTS-2) || (lMax < lMin*2) ) + if( (iBot < 2) || (iTop > self->maxpts-2) || (lMax < lMin*2) ) { SCWrite(pCon,"ERROR: no peak found!",eError); return 0; @@ -420,6 +421,7 @@ pNew->pDes = CreateDescriptor("Maximizer"); pNew->pCount = pCom->pData; pNew->i360 = 0; + pNew->maxpts = 100; AddCommand(pSics,"max",MaximizeAction,MaxKill,pNew); return 1; @@ -434,7 +436,7 @@ double dVal; float fStep, fPreset; CounterMode eCount; - int iRet; + int iRet, iVal; self = (pMax)pData; assert(self); @@ -446,9 +448,56 @@ return 1; } + /* enough arguments ?*/ if(argc < 5) { + if(argc > 1) + { + strtolower(argv[1]); + if(strcmp(argv[1],"in360") == 0) + { + if(argc > 2) + { + iVal = atoi(argv[2]); + if(iVal != 0 && iVal != 1) { + SCWrite(pCon,"ERROR: only 0, 1 allowed for in360",eError); + return 0; + } + self->i360 = iVal; + SCSendOK(pCon); + return 1; + } + else + { + snprintf(pBueffel,255,"max.in360 = %d", self->i360); + SCWrite(pCon,pBueffel,eValue); + return 1; + } + } + if(strcmp(argv[1],"maxpts") == 0) + { + if(argc > 2) + { + iVal = atoi(argv[2]); + if(iVal < 10 || iVal > 100) { + SCWrite(pCon,"ERROR: maxpst must be between 10 and 100", + eError); + return 0; + } + self->maxpts = iVal; + SCSendOK(pCon); + return 1; + } + else + { + snprintf(pBueffel,255,"max.maxpts = %d", self->maxpts); + SCWrite(pCon,pBueffel,eValue); + return 1; + } + } + + } SCWrite(pCon,"ERROR: Insufficient number of arguments to max", eError); return 0; diff --git a/multicounter.c b/multicounter.c index a0ccdde6..8ab9a351 100644 --- a/multicounter.c +++ b/multicounter.c @@ -26,6 +26,7 @@ #include "splitter.h" #define MAXSLAVE 16 +#define NOCOUNTERS -2727 /*=============== code for the driver ======================================*/ typedef struct { void *slaveData[MAXSLAVE]; @@ -103,6 +104,11 @@ static int MMCCStatus(void *pData, SConnection *pCon){ } assert(self); + if(self->nSlaves == 0) { + pCount->pDriv->iErrorCode = NOCOUNTERS; + return HWFault; + } + status = self->slaves[0]->CheckCountStatus(self->slaveData[0],pCon); if(status == HWIdle || status == HWFault){ /* @@ -256,6 +262,43 @@ static void MMCCParameter(void *pData, float fPreset, CounterMode eMode ){ eMode); } } +/*======================= Driver Interface ==============================*/ +static int MultiCounterSet(struct __COUNTER *self, char *name, + int iCter, float fVal){ + + return 0; +} +/*-----------------------------------------------------------------------*/ +static int MultiCounterGet(struct __COUNTER *self, char *name, + int iCter, float *fVal){ + + return 0; +} +/*-----------------------------------------------------------------------*/ +static int MultiCounterSend(struct __COUNTER *self, char *pText, + char *reply, int replylen){ + + strncpy(reply,"NOT Implemented",replylen); + return 0; +} +/*---------------------------------------------------------------------*/ +static int MultiCounterError(struct __COUNTER *pData, int *iCode, + char *error, int errlen){ + pCounter pCount = NULL; + + pCount = (pCounter)pData; + + if(pCount->pDriv->iErrorCode == NOCOUNTERS){ + strncpy(error,"NO counters configured!",errlen); + } else { + strncpy(error,"Not Implemented", errlen); + } + return COTERM; +} +/*----------------------------------------------------------------------*/ +static int MultiCounterFix(struct __COUNTER *self, int iCode){ + return COTERM; +} /*=============== Interpreter Interface ================================ */ int MultiCounterAction(SConnection *pCon, SicsInterp *pSics, void *pData, int argc, char *argv[]){ @@ -327,6 +370,12 @@ int MakeMultiCounter(SConnection *pCon, SicsInterp *pSics, SCWrite(pCon,"ERROR: out of memory in MakeMultiCounter",eError); return 0; } + pDriv->Get = MultiCounterGet; + pDriv->GetError = MultiCounterError; + pDriv->TryAndFixIt = MultiCounterFix; + pDriv->Set = MultiCounterSet; + pDriv->Send += MultiCounterSend; /* assign interface functions diff --git a/optimise.c b/optimise.c index 226f5001..89a3ad89 100644 --- a/optimise.c +++ b/optimise.c @@ -596,7 +596,7 @@ static int findDirection(pOptimise self, pOVarEntry pOvar, SConnection *pCon) self->pScanner->fPreset = self->fPreset; direction = findDirection(self,pOvar, pCon); - if(direction < 0){ + if(direction < -1){ return direction; } /* diff --git a/simcter.c b/simcter.c index 1b501e67..9963c191 100644 --- a/simcter.c +++ b/simcter.c @@ -337,7 +337,7 @@ static float FAILRATE; pData->lEnd = 0; pData->iPause = 0; pRes->pData = (void *)pData; - pRes->iNoOfMonitors = 2; + pRes->iNoOfMonitors = 8; pRes->fTime = 0; /* assign functions */ diff --git a/tasdrive.c b/tasdrive.c index 1c0bf6b7..de3fff85 100644 --- a/tasdrive.c +++ b/tasdrive.c @@ -336,8 +336,8 @@ static int startMotors(ptasMot self, tasAngles angles, return status; } } + writeMotPos(pCon,silent,"ach",val, curve); } - writeMotPos(pCon,silent,"ach",val, curve); } if(driveQ == 0){ diff --git a/tasub.c b/tasub.c index e498011b..462e96f0 100644 --- a/tasub.c +++ b/tasub.c @@ -9,6 +9,7 @@ #include #include "sics.h" #include "lld.h" +#include "trigd.h" #include "tasub.h" #include "tasdrive.h" /*------------------- motor indexes in motor data structure ---------*/ @@ -875,8 +876,100 @@ static void listDiagnostik(ptasUB self, SConnection *pCon){ } } /*------------------------------------------------------------------*/ -static int calcUB(ptasUB self, SConnection *pCon, SicsInterp *pSics, +static int calcAuxUB(ptasUB self, SConnection *pCon, SicsInterp *pSics, int argc, char *argv[]){ + int status; + tasReflection r1, r2; + char pBueffel[256]; + MATRIX UB = NULL, B = NULL; + + if(argc < 5){ + SCWrite(pCon, + "ERROR: not enough arguments for UB calculation, need HKJL of second plane vector", + eError); + return 0; + } + + if(!SCMatchRights(pCon,usUser)){ + return 0; + } + + status = Tcl_GetDouble(InterpGetTcl(pSics),argv[2],&r2.qe.qh); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[2]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(InterpGetTcl(pSics),argv[3],&r2.qe.qk); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[3]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + status = Tcl_GetDouble(InterpGetTcl(pSics),argv[4],&r2.qe.ql); + if(status != TCL_OK){ + snprintf(pBueffel,255,"ERROR: failed to convert %s to number",argv[4]); + SCWrite(pCon,pBueffel,eError); + return 0; + } + + status = findReflection(self->reflectionList, 0,&r1); + if(status != 1){ + snprintf(pBueffel,255,"ERROR: cannot find first reflection"); + SCWrite(pCon,pBueffel,eError); + return 0; + } + B = mat_creat(3,3,ZERO_MATRIX); + if(B == NULL){ + SCWrite(pCon,"ERROR: out of memory creating B matrix",eError); + return 0; + } + status = calculateBMatrix(self->cell,B); + if(status < 0){ + SCWrite(pCon,"ERROR: bad cell constants, no volume",eError); + mat_free(B); + return 0; + } + + status = makeAuxReflection(B, r1, &r2,self->machine.ss_sample); + mat_free(B); + if(status < 0){ + SCWrite(pCon,"ERROR: out of memory in makeAuxUB",eError); + return 0; + } + + UB = calcTasUBFromTwoReflections(self->cell,r1,r2,&status); + if(UB == NULL){ + switch(status){ + case UBNOMEMORY: + SCWrite(pCon,"ERROR: out of memory calculating UB matrix",eError); + break; + case REC_NO_VOLUME: + SCWrite(pCon,"ERROR: bad cell constants, no volume",eError); + break; + } + return 0; + } + if(mat_det(UB) < .000001){ + SCWrite(pCon,"ERROR: invalid UB matrix, check reflections",eError); + return 0; + } + if(self->machine.UB != NULL){ + mat_free(self->machine.UB); + } + if(self->machine.planeNormal != NULL){ + mat_free(self->machine.planeNormal); + } + self->machine.UB = UB; + self->machine.planeNormal = calcPlaneNormal(r1,r2); + self->ubValid = 1; + SCparChange(pCon); + SCSendOK(pCon); + return 1; +} +/*------------------------------------------------------------------*/ +static int calcUB(ptasUB self, SConnection *pCon, SicsInterp *pSics, + int argc, char *argv[]){ int idx1, idx2, status; tasReflection r1, r2; char pBueffel[256]; @@ -884,8 +977,8 @@ static int calcUB(ptasUB self, SConnection *pCon, SicsInterp *pSics, if(argc < 4){ SCWrite(pCon, - "ERROR: not enough arguments for UB calculation, need index of two reflections", - eError); + "ERROR: not enough arguments for UB calculation, need index of two reflections", + eError); return 0; } @@ -955,10 +1048,20 @@ static int calcUB(ptasUB self, SConnection *pCon, SicsInterp *pSics, /*-----------------------------------------------------------------*/ static int calcUBFromCell(ptasUB self, SConnection *pCon){ MATRIX B, U, UB; + tasReflection r1; int status; B = mat_creat(3,3,UNIT_MATRIX); U = mat_creat(3,3,UNIT_MATRIX); + status = findReflection(self->reflectionList, 0,&r1); + if(status == 1) { + /* + U[0][0] = Cosd(r1.angles.a3); + U[0][1] = -Sind(r1.angles.a3); + U[1][0] = Sind(r1.angles.a3); + U[1][1] = Cosd(r1.angles.a3); + */ + } if(B == NULL || U == NULL){ SCWrite(pCon,"ERROR: out of memory in calcUBFromCell",eError); return 0; @@ -1459,6 +1562,8 @@ int TasUBWrapper(SConnection *pCon,SicsInterp *pSics, void *pData, return 1; } else if(strcmp(argv[1],"makeub") == 0){ return calcUB(self,pCon,pSics,argc,argv); + } else if(strcmp(argv[1],"makeauxub") == 0){ + return calcAuxUB(self,pCon,pSics,argc,argv); } else if(strcmp(argv[1],"makeubfromcell") == 0){ return calcUBFromCell(self,pCon); } else if(strcmp(argv[1],"calcang") == 0){ diff --git a/tasublib.c b/tasublib.c index 4d8a45f3..07af0d1b 100644 --- a/tasublib.c +++ b/tasublib.c @@ -139,6 +139,38 @@ static double calcTheta(double ki, double kf, double two_theta){ return rtan(ABS(ki) - ABS(kf)*Cosd(two_theta), ABS(kf)*Sind(two_theta))/DEGREE_RAD; } +/*-------------------------------------------------------------------------*/ +double tasAngleBetweenReflections(MATRIX B, tasReflection r1, tasReflection r2){ + MATRIX chi1, chi2, h1, h2; + double angle; + + h1 = makeVector(); + if(h1 == NULL){ + return -9999.99; + } + h1[0][0] = r1.qe.qh; + h1[1][0] = r1.qe.qk; + h1[2][0] = r1.qe.ql; + + h2 = makeVector(); + if(h2 == NULL){ + return -999.99; + } + h2[0][0] = r2.qe.qh; + h2[1][0] = r2.qe.qk; + h2[2][0] = r2.qe.ql; + + 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; +} /*--------------------------------------------------------------------*/ static MATRIX uFromAngles(double om, double sgu, double sgl){ MATRIX u; @@ -161,6 +193,59 @@ static MATRIX calcTasUVectorFromAngles(tasReflection r){ om = r.angles.a3 - theta; return uFromAngles(om,r.angles.sgu, r.angles.sgl); } +/*-----------------------------------------------------------------------------*/ +static MATRIX tasReflectionToQC(tasQEPosition r, MATRIX UB){ + MATRIX Q, QC; + + Q = makeVector(); + if(Q == NULL){ + return NULL; + } + vectorSet(Q,0,r.qh); + vectorSet(Q,1,r.qk); + vectorSet(Q,2,r.ql); + QC = mat_mul(UB,Q); + killVector(Q); + return QC; +} +/*-----------------------------------------------------------------*/ +int makeAuxReflection(MATRIX B, tasReflection r1, tasReflection *r2, + int ss){ + double theta, om, q, cos2t; + MATRIX QC; + + r2->angles = r1.angles; + r2->qe.ki = r1.qe.ki; + r2->qe.kf= r1.qe.kf; + + theta = calcTheta(r1.qe.ki,r1.qe.kf,r1.angles.sample_two_theta); + om = r1.angles.a3 - theta; + om += tasAngleBetweenReflections(B,r1,*r2); + + QC = tasReflectionToHC(r2->qe,B); + if(QC == NULL){ + return UBNOMEMORY; + } + + q = vectorLength(QC); + q = 2.*PI*vectorLength(QC); + cos2t = (r1.qe.ki*r1.qe.ki + r1.qe.kf*r1.qe.kf - q*q)/ + (2. * ABS(r1.qe.ki) * ABS(r1.qe.kf)); + if(ABS(cos2t) > 1.){ + killVector(QC); + return TRIANGLENOTCLOSED; + } + r2->angles.sample_two_theta = ss*Acosd(cos2t); + theta = calcTheta(r1.qe.ki,r1.qe.kf,r2->angles.sample_two_theta); + r2->angles.a3 = om + theta; + r2->angles.a3 -= 180.; + if(r2->angles.a3 < -180.){ + r2->angles.a3 += 360.; + } + mat_free(QC); + + return 1; +} /*-------------------------------------------------------------------*/ MATRIX calcPlaneNormal(tasReflection r1, tasReflection r2){ MATRIX u1 = NULL, u2 = NULL, planeNormal = NULL; @@ -326,21 +411,6 @@ static MATRIX buildTVMatrix(MATRIX U1V, MATRIX U2V){ killVector(T3V); return T; } -/*-----------------------------------------------------------------------------*/ -static MATRIX tasReflectionToQC(tasQEPosition r, MATRIX UB){ - MATRIX Q, QC; - - Q = makeVector(); - if(Q == NULL){ - return NULL; - } - vectorSet(Q,0,r.qh); - vectorSet(Q,1,r.qk); - vectorSet(Q,2,r.ql); - QC = mat_mul(UB,Q); - killVector(Q); - return QC; -} /*----------------------------------------------------------------------------*/ static MATRIX buildRMatrix(MATRIX UB, MATRIX planeNormal, tasQEPosition qe, int *errorCode){ diff --git a/tasublib.h b/tasublib.h index 1bbfe463..35fb59ff 100644 --- a/tasublib.h +++ b/tasublib.h @@ -124,6 +124,19 @@ double maCalcHorizontalCurvature(maCrystal data, double two_theta); */ double maCalcK(maCrystal data, double two_theta); /*======================= reciprocal space =============================*/ +/** + * make an auxiliary reflection which has the same sgu and sgl as r1, but + * an omega which is adjusted to the angle difference between r1 and + * the target. This is useful for generating an auxilairy UB during + * alignment. + * @param B The B matrix + * @param r1 The first reflection + * @param r2 a pointer to the second reflection, QH, QK, QL initialized + * @param ss The scattering sense at the sample + * @return 1 on success, a negative error code on failure. + */ +int makeAuxReflection(MATRIX B, tasReflection r1, + tasReflection *r2, int ss); /** * calculate a UB from two reflections and the cell. * @param cell The lattice constant of the crystal