linse/sics
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

- added makeauxub to tasub

Browse Source 
- fixes for MARS
- extended maximize to honour maxpts and the in360  flag
- added regression histmem driver
This commit is contained in:
koennecke
2006-10-20 14:54:39 +00:00
parent 9f668b7681
commit 4f069341f5
14 changed files with 623 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
counter.c
View File

@ -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",

7
fitcenter.c
View File

@ -468,6 +468,13 @@
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;
}
}
/* print results */

18
histmem.c
View File

@ -68,6 +68,10 @@
/*
#define LOADDEBUG 1
*/
/*
* from histregress.c
*/
extern pHistDriver CreateRegressHM(pStringDict pOpt);
/*------------------------------------------------------------------------*/
static int HistHalt(void *pData)
{
@ -454,6 +458,10 @@
{
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)
{

273
histregress.c Normal file
View File

@ -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 <string.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#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;
}

4
hkl.c
View File

@ -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;

2
make_gen
View File

@ -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

73
maximize.c
View File

@ -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;

49
multicounter.c
View File

@ -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

2
optimise.c
View File

@ -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;
}
/*

2
simcter.c
View File

@ -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 */

2
tasdrive.c
View File

@ -336,9 +336,9 @@ static int startMotors(ptasMot self, tasAngles angles,
return status;
}
}
}
writeMotPos(pCon,silent,"ach",val, curve);
}
}
if(driveQ == 0){
return OKOK;

105
tasub.c
View File

@ -9,6 +9,7 @@
#include <assert.h>
#include "sics.h"
#include "lld.h"
#include "trigd.h"
#include "tasub.h"
#include "tasdrive.h"
/*------------------- motor indexes in motor data structure ---------*/
@ -875,6 +876,98 @@ static void listDiagnostik(ptasUB self, SConnection *pCon){
}
}
/*------------------------------------------------------------------*/
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;
@ -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){

100
tasublib.c
View File

@ -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){

13
tasublib.h
View File

@ -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