diff --git a/reflist.c b/reflist.c new file mode 100644 index 00000000..d9cd95c2 --- /dev/null +++ b/reflist.c @@ -0,0 +1,878 @@ +/** + * Reflection list: a configurable list of reflections and + * their setting angles. To be used for four circle and possibly + * TAS diffraction. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, July 2008 + */ +#include +#include +#include "reflist.h" +#include "sicshipadaba.h" +#include "stptok.h" + +#define INAME "index" +#define ANAME "angles" +#define FLAG "flag" +static char undef[] = "Undefined"; +#define IDXFMT " %8.4f" +#define ANGFMT " %8.2f" +/*----------------------------------------------------------------------*/ +typedef struct { + int idxCount; + int angCount; + int count; +} RLPriv, *pRLPriv; +/*---------------------------------------------------------------------- + * The point of the code below and the callback is to update the + * configuration counts whenever the corresponding parameters + * are changed. + * ---------------------------------------------------------------------*/ +static void UpdateConfiguration(pSICSOBJ refl) +{ + pRLPriv priv = refl->pPrivate; + hdbValue v; + char *pPtr = NULL, token[50]; + + + SICSHdbGetPar(refl, NULL, "indexheader", &v); + pPtr = v.v.text; + priv->idxCount = 0; + while ((pPtr = stptok(pPtr, token, 50, ",")) != NULL) { + priv->idxCount++; + } + + SICSHdbGetPar(refl, NULL, "anglesheader", &v); + pPtr = v.v.text; + priv->angCount = 0; + while ((pPtr = stptok(pPtr, token, 50, ",")) != NULL) { + priv->angCount++; + } + ClearReflectionList(refl); +} + +/*-----------------------------------------------------------------------*/ +static hdbCallbackReturn CalcConfiguration(pHdb node, void *userData, + pHdbMessage mm) +{ + + if (GetHdbSetMessage(mm) != NULL) { + /* + * This is after the normal SICS parameter SetUpdateCallback + * which already updated the node + */ + UpdateConfiguration((pSICSOBJ) userData); + } + return hdbContinue; +} + +/*-----------------------------------------------------------------------*/ +static int SaveRefList(void *data, char *name, FILE * fd) +{ + int i; + double ang[4], hkl[3]; + pSICSOBJ self = data; + + SaveSICSOBJ(data, name, fd); + for (i = 0; i < ReflectionListCount(self); i++) { + memset(hkl, 0, 3 * sizeof(double)); + memset(ang, 0, 4 * sizeof(double)); + GetRefIndex(self, i, hkl); + GetRefAngles(self, i, ang); + fprintf(fd, "%s addax %f %f %f %f %f %f %f\n", name, + hkl[0], hkl[1], hkl[2], ang[0], ang[1], ang[2], ang[3]); + } + return 1; +} + +/*-----------------------------------------------------------------------*/ +pSICSOBJ MakeReflectionListInt(char *name) +{ + pSICSOBJ pNew = NULL; + pHdb node = NULL; + hdbValue v; + pRLPriv priv = NULL; + + pNew = MakeSICSOBJ(name, "ReflectionList"); + priv = malloc(sizeof(RLPriv)); + if (pNew == NULL || priv == NULL) { + return NULL; + } + + pNew->pPrivate = priv; + pNew->KillPrivate = DefaultKill; + pNew->pDes->SaveStatus = SaveRefList; + + v = MakeHdbText("h,k,l"); + node = MakeSICSHdbPar("indexheader", usMugger, v); + AppendHipadabaCallback(node, MakeHipadabaCallback(CalcConfiguration, + pNew, NULL)); + AddHipadabaChild(pNew->objectNode, node, NULL); + priv->idxCount = 3; + + v = MakeHdbText("stt,om,chi,phi"); + node = MakeSICSHdbPar("anglesheader", usMugger, v); + AppendHipadabaCallback(node, MakeHipadabaCallback(CalcConfiguration, + pNew, NULL)); + AddHipadabaChild(pNew->objectNode, node, NULL); + priv->angCount = 4; + priv->count = 0; + + v = makeHdbValue(HIPNONE, 0); + node = MakeSICSHdbPar("list", usMugger, v); + AddHipadabaChild(pNew->objectNode, node, NULL); + + return pNew; +} + +/*-----------------------------------------------------------------------*/ +static int MakeCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + int status; + + status = NewListReflection(self); + if (status >= 0) { + SCPrintf(pCon, eValue, "newref = %4.4d", status); + return 1; + } else { + SCWrite(pCon, "ERROR: failed to create reflection", eError); + return 0; + } + return 0; +} + +/*-----------------------------------------------------------------------*/ +static int ClearCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + ClearReflectionList(self); + return 1; +} + +/*----------------------------------------------------------------------*/ +static int AddIndexCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + char path[80]; + int status; + pHdb node = NULL; + hdbValue v; + + if (nPar < 1) { + SCWrite(pCon, "ERROR: need indexes as arguments", eError); + return 0; + } + + status = NewListReflection(self); + snprintf(path, 80, "list/%4.4d/%s", status, INAME); + node = GetHipadabaNode(self->objectNode, path); + assert(node != NULL); + + cloneHdbValue(&node->value, &v); + if (!readHdbValue(&v, par[0]->value.v.text, path, 80)) { + SCPrintf(pCon, eError, "ERROR: failed to convert %s to indices", + par[0]->value.v.text); + return 0; + } + SetHipadabaPar(node, v, pCon); + ReleaseHdbValue(&v); + return 1; +} + +/*----------------------------------------------------------------------*/ +void AddRefIdx(pSICSOBJ refl, double hkl[]) +{ + int status; + hdbValue v; + char path[80]; + pHdb node; + + status = NewListReflection(refl); + snprintf(path, 80, "list/%4.4d/%s", status, INAME); + node = GetHipadabaNode(refl->objectNode, path); + assert(node != NULL); + + v = MakeHdbFloatArray(3, hkl); + SetHipadabaPar(node, v, NULL); +} + +/*----------------------------------------------------------------------*/ +void AddRefIdxAng(pSICSOBJ refl, double hkl[], double ang[]) +{ + int status; + hdbValue v; + char path[80]; + pHdb node; + + status = NewListReflection(refl); + snprintf(path, 80, "list/%4.4d/%s", status, INAME); + node = GetHipadabaNode(refl->objectNode, path); + assert(node != NULL); + + v = MakeHdbFloatArray(3, hkl); + SetHipadabaPar(node, v, NULL); + + snprintf(path, 80, "list/%4.4d/%s", status, ANAME); + node = GetHipadabaNode(refl->objectNode, path); + assert(node != NULL); + + v = MakeHdbFloatArray(3, ang); + SetHipadabaPar(node, v, NULL); +} + +/*----------------------------------------------------------------------*/ +static int AddAnglesCmd(pSICSOBJ self, SConnection * pCon, + pHdb commandNode, pHdb par[], int nPar) +{ + char path[80]; + int status; + pHdb node = NULL; + hdbValue v; + + if (nPar < 1) { + SCWrite(pCon, "ERROR: need angles as arguments", eError); + return 0; + } + + status = NewListReflection(self); + snprintf(path, 80, "list/%4.4d/%s", status, ANAME); + node = GetHipadabaNode(self->objectNode, path); + assert(node != NULL); + + cloneHdbValue(&node->value, &v); + if (!readHdbValue(&v, par[0]->value.v.text, path, 80)) { + SCPrintf(pCon, eError, "ERROR: failed to convert %s to angles", + par[0]->value.v.text); + return 0; + } + SetHipadabaPar(node, v, pCon); + ReleaseHdbValue(&v); + return 1; +} + +/*----------------------------------------------------------------------*/ +static int AddIndexesAnglesCmd(pSICSOBJ self, SConnection * pCon, + pHdb commandNode, pHdb par[], int nPar) +{ + char path[80]; + int status, i; + pHdb node = NULL; + hdbValue v; + char *pPtr = NULL; + pRLPriv priv = self->pPrivate; + + if (nPar < 1) { + SCWrite(pCon, "ERROR: need indexes/angles as arguments", eError); + return 0; + } + + status = NewListReflection(self); + snprintf(path, 80, "list/%4.4d/%s", status, INAME); + node = GetHipadabaNode(self->objectNode, path); + assert(node != NULL); + + cloneHdbValue(&node->value, &v); + if (!readHdbValue(&v, par[0]->value.v.text, path, 80)) { + SCPrintf(pCon, eError, "ERROR: failed to convert %s to indices", + par[0]->value.v.text); + return 0; + } + SetHipadabaPar(node, v, pCon); + ReleaseHdbValue(&v); + + /* advance pointer beyond indices */ + pPtr = par[0]->value.v.text; + for (i = 0; i < priv->idxCount; i++) { + pPtr = stptok(pPtr, path, 80, " "); + } + + snprintf(path, 80, "list/%4.4d/%s", status, ANAME); + node = GetHipadabaNode(self->objectNode, path); + assert(node != NULL); + + cloneHdbValue(&node->value, &v); + if (!readHdbValue(&v, pPtr, path, 80)) { + SCPrintf(pCon, eError, "ERROR: failed to convert %s to angles", + par[0]->value.v.text); + return 0; + } + SetHipadabaPar(node, v, pCon); + ReleaseHdbValue(&v); + + return 1; +} + +/*-----------------------------------------------------------------------*/ +static int DelCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + pHdb node = NULL; + char path[132]; + + if (nPar < 1) { + SCWrite(pCon, "ERROR: need id to delete as argument", eError); + return 0; + } + snprintf(path, 132, "list/%s", par[0]->value.v.text); + node = GetHipadabaNode(self->objectNode, path); + if (node != NULL) { + DeleteHipadabaNode(node, pCon); + return 1; + } + SCWrite(pCon, "ERROR: reflection not found", eError); + return 0; +} + +/*-----------------------------------------------------------------------*/ +static int FlagCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + pHdb node = NULL; + char path[132]; + + if (nPar < 1) { + SCWrite(pCon, "ERROR: need id to flag for", eError); + return 0; + } + snprintf(path, 132, "list/%s/%s", par[0]->value.v.text, FLAG); + node = GetHipadabaNode(self->objectNode, path); + if (node == NULL) { + SCPrintf(pCon, eError, "ERROR: bad ID %s", par[0]->value.v.text); + return 0; + } + if (nPar < 2) { + SCPrintf(pCon, eValue, "%s.%s.flag = %d", self->pDes->name, + par[0]->value.v.text, node->value.v.intValue); + return 1; + } else { + SetHipadabaPar(node, par[1]->value, pCon); + return 1; + } + return 0; +} + +/*-----------------------------------------------------------------------*/ +static int ShowCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + pHdb node = NULL; + char path[132], num[20]; + char data[1024]; + hdbValue v; + int i; + + if (nPar < 1) { + SCWrite(pCon, "ERROR: need id for reflection to show", eError); + return 0; + } + snprintf(path, 132, "list/%s/%s", par[0]->value.v.text, FLAG); + node = GetHipadabaNode(self->objectNode, path); + if (node == NULL) { + SCPrintf(pCon, eError, "ERROR: bad ID %s", par[0]->value.v.text); + return 0; + } + GetHipadabaPar(node, &v, pCon); + snprintf(data, 1024, "%s %d", par[0]->value.v.text, v.v.intValue); + + snprintf(path, 132, "list/%s/%s", par[0]->value.v.text, INAME); + node = GetHipadabaNode(self->objectNode, path); + assert(node != NULL); + GetHipadabaPar(node, &v, pCon); + for (i = 0; i < v.arrayLength; i++) { + snprintf(num, 20, " %f", v.v.floatArray[i]); + strncat(data, num, 1024); + } + + snprintf(path, 132, "list/%s/%s", par[0]->value.v.text, ANAME); + node = GetHipadabaNode(self->objectNode, path); + assert(node != NULL); + GetHipadabaPar(node, &v, pCon); + for (i = 0; i < v.arrayLength; i++) { + snprintf(num, 20, " %f", v.v.floatArray[i]); + strncat(data, num, 1024); + } + SCWrite(pCon, data, eValue); + return 1; +} + +/*-----------------------------------------------------------------------*/ +static int ListCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + char buffer[132], path[132], *pPtr, token[20], entry[20]; + pRLPriv priv = self->pPrivate; + int i; + pHdb node = NULL, data; + pDynString txt; + + SCStartBuffering(pCon); + /* + * build header + */ + strcpy(buffer, " ID FLAG"); + node = GetHipadabaNode(self->objectNode, "indexheader"); + assert(node != NULL); + pPtr = node->value.v.text; + while ((pPtr = stptok(pPtr, token, 20, ",")) != NULL) { + snprintf(entry, 20, " %8s", token); + strncat(buffer, entry, 132); + } + node = GetHipadabaNode(self->objectNode, "anglesheader"); + assert(node != NULL); + pPtr = node->value.v.text; + while ((pPtr = stptok(pPtr, token, 20, ",")) != NULL) { + snprintf(entry, 20, " %8s", token); + strncat(buffer, entry, 132); + } + SCWrite(pCon, buffer, eValue); + + node = GetHipadabaNode(self->objectNode, "list"); + node = node->child; + while (node != NULL) { + data = GetHipadabaNode(node, FLAG); + assert(data != NULL); + snprintf(buffer, 132, " %5s %4d", node->name, data->value.v.intValue); + data = GetHipadabaNode(node, INAME); + for (i = 0; i < priv->idxCount; i++) { + sprintf(entry, IDXFMT, data->value.v.floatArray[i]); + strncat(buffer, entry, 132); + } + data = GetHipadabaNode(node, ANAME); + for (i = 0; i < priv->angCount; i++) { + sprintf(entry, ANGFMT, data->value.v.floatArray[i]); + strncat(buffer, entry, 132); + } + SCWrite(pCon, buffer, eValue); + node = node->next; + } + + txt = SCEndBuffering(pCon); + SCWrite(pCon, GetCharArray(txt), eValue); + return 1; +} + +/*-----------------------------------------------------------------------*/ +static int NamesCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + char buffer[132]; + pHdb node = NULL; + pDynString txt; + + SCStartBuffering(pCon); + + node = GetHipadabaNode(self->objectNode, "list"); + node = node->child; + while (node != NULL) { + snprintf(buffer, 132, "%s", node->name); + SCWrite(pCon, buffer, eValue); + node = node->next; + } + + txt = SCEndBuffering(pCon); + SCWrite(pCon, GetCharArray(txt), eValue); + return 1; +} + +/*----------------------------------------------------------------------*/ +static int SetIndexCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + char path[80]; + int status; + pHdb node = NULL; + hdbValue v; + char *pPtr = NULL, refl[20]; + + if (nPar < 1) { + SCWrite(pCon, "ERROR: need id and indexes as arguments", eError); + return 0; + } + + pPtr = par[0]->value.v.text; + pPtr = stptok(pPtr, refl, 20, " "); + + snprintf(path, 80, "list/%s/%s", refl, INAME); + node = GetHipadabaNode(self->objectNode, path); + assert(node != NULL); + + cloneHdbValue(&node->value, &v); + if (!readHdbValue(&v, pPtr, path, 80)) { + SCPrintf(pCon, eError, "ERROR: failed to convert %s to indices", + par[0]->value.v.text); + return 0; + } + SetHipadabaPar(node, v, pCon); + ReleaseHdbValue(&v); + return 1; +} + +/*----------------------------------------------------------------------*/ +static int SetAnglesCmd(pSICSOBJ self, SConnection * pCon, + pHdb commandNode, pHdb par[], int nPar) +{ + char path[80]; + int status; + pHdb node = NULL; + hdbValue v; + char *pPtr = NULL, refl[20]; + + if (nPar < 1) { + SCWrite(pCon, "ERROR: need id and angles as arguments", eError); + return 0; + } + + pPtr = par[0]->value.v.text; + pPtr = stptok(pPtr, refl, 20, " "); + + snprintf(path, 80, "list/%s/%s", refl, ANAME); + node = GetHipadabaNode(self->objectNode, path); + if (node == NULL) { + SCPrintf(pCon, eError, "ERROR: reflection with ID %s not found", refl); + return 0; + } + + cloneHdbValue(&node->value, &v); + if (!readHdbValue(&v, pPtr, path, 80)) { + SCPrintf(pCon, eError, "ERROR: failed to convert %s to angles", + par[0]->value.v.text); + return 0; + } + SetHipadabaPar(node, v, pCon); + ReleaseHdbValue(&v); + return 1; +} + +/*-----------------------------------------------------------------------*/ +int MakeReflectionList(SConnection * pCon, SicsInterp * pSics, + void *data, int argc, char *argv[]) +{ + pSICSOBJ pNew = NULL; + pHdb cmd; + + if (argc < 2) { + SCWrite(pCon, "ERROR: need name of reflection list", eError); + return 0; + } + pNew = CreateReflectionList(pCon, pSics, argv[1]); + if (pNew == NULL) { + return 0; + } + return 1; +} + +/*------------------------------------------------------------------------*/ +pSICSOBJ CreateReflectionList(SConnection * pCon, SicsInterp * pSics, + char *name) +{ + pSICSOBJ pNew = NULL; + pHdb cmd; + + pNew = MakeReflectionListInt(name); + if (pNew == NULL) { + SCWrite(pCon, "ERROR: failed to create reflection list", eError); + return 0; + } + + /* command initialization */ + AddSICSHdbPar(pNew->objectNode, "make", usUser, MakeSICSFunc(MakeCmd)); + + AddSICSHdbPar(pNew->objectNode, "clear", usUser, MakeSICSFunc(ClearCmd)); + + AddSICSHdbPar(pNew->objectNode, "print", usUser, MakeSICSFunc(ListCmd)); + + AddSICSHdbPar(pNew->objectNode, "names", usUser, MakeSICSFunc(NamesCmd)); + + cmd = + AddSICSHdbPar(pNew->objectNode, "del", usUser, MakeSICSFunc(DelCmd)); + AddSICSHdbPar(cmd, "id", usUser, MakeHdbText("")); + + cmd = + AddSICSHdbPar(pNew->objectNode, "addx", usUser, + MakeSICSFunc(AddIndexCmd)); + AddSICSHdbPar(cmd, "args", usUser, MakeHdbText("")); + + cmd = + AddSICSHdbPar(pNew->objectNode, "setx", usUser, + MakeSICSFunc(SetIndexCmd)); + AddSICSHdbPar(cmd, "args", usUser, MakeHdbText("")); + + cmd = + AddSICSHdbPar(pNew->objectNode, "adda", usUser, + MakeSICSFunc(AddAnglesCmd)); + AddSICSHdbPar(cmd, "args", usUser, MakeHdbText("")); + + cmd = + AddSICSHdbPar(pNew->objectNode, "seta", usUser, + MakeSICSFunc(SetAnglesCmd)); + AddSICSHdbPar(cmd, "args", usUser, MakeHdbText("")); + + cmd = AddSICSHdbPar(pNew->objectNode, "addax", usUser, + MakeSICSFunc(AddIndexesAnglesCmd)); + AddSICSHdbPar(cmd, "args", usUser, MakeHdbText("")); + + cmd = AddSICSHdbPar(pNew->objectNode, "flag", usUser, + MakeSICSFunc(FlagCmd)); + AddSICSHdbPar(cmd, "id", usUser, MakeHdbText("")); + AddSICSHdbPar(cmd, "value", usUser, MakeHdbInt(0)); + + cmd = AddSICSHdbPar(pNew->objectNode, "show", usUser, + MakeSICSFunc(ShowCmd)); + AddSICSHdbPar(cmd, "id", usUser, MakeHdbText("")); + + AddCommand(pSics, name, InterInvokeSICSOBJ, KillSICSOBJ, pNew); + return pNew; +} + +/*------------------------------------------------------------------------*/ +void ClearReflectionList(pSICSOBJ refl) +{ + pHdb node = NULL; + hdbValue v; + pRLPriv priv = refl->pPrivate; + + node = GetHipadabaNode(refl->objectNode, "list"); + assert(node); + DeleteHipadabaNode(node, NULL); + v = makeHdbValue(HIPNONE, 0); + node = MakeSICSHdbPar("list", usMugger, v); + AddHipadabaChild(refl->objectNode, node, NULL); + priv->count = 0; +} + +/*-------------------------------------------------------------------------*/ +int ReflectionListCount(pSICSOBJ refl) +{ + pHdb node = NULL; + + node = GetHipadabaNode(refl->objectNode, "list"); + assert(node); + return CountHdbChildren(node); +} + +/*-------------------------------------------------------------------------*/ +int NewListReflection(pSICSOBJ refl) +{ + int count; + char path[20]; + pHdb listNode = NULL, node = NULL, refNode = NULL; + hdbValue v; + pRLPriv priv = refl->pPrivate; + + listNode = GetHipadabaNode(refl->objectNode, "list"); + assert(listNode); + + count = priv->count; + priv->count++; + snprintf(path, 20, "%4.4d", count); + node = MakeHipadabaNode(path, HIPNONE, 0); + AddHipadabaChild(listNode, node, NULL); + + v = makeHdbValue(HIPFLOATAR, priv->idxCount); + AddSICSHdbPar(node, INAME, usUser, v); + + v = makeHdbValue(HIPFLOATAR, priv->angCount); + AddSICSHdbPar(node, ANAME, usUser, v); + + v = makeHdbValue(HIPINT, 0); + AddSICSHdbPar(node, FLAG, usUser, v); + return count; +} + +/*--------------------------------------------------------------------------*/ +static pHdb findReflection(pSICSOBJ refl, int idx) +{ + pHdb node = NULL; + int count = 0; + + if (idx < 0) { + return node; + } + + node = GetHipadabaNode(refl->objectNode, "list"); + assert(node); + + node = node->child; + while (node != NULL && count != idx) { + count++; + node = node->next; + } + return node; +} + +/*---------------------------------------------------------------------------*/ +void DelListReflection(pSICSOBJ refl, int idx) +{ + pHdb node = NULL; + + node = findReflection(refl, idx); + if (node != NULL) { + DeleteHipadabaNode(node, NULL); + } +} + +/*-----------------------------------------------------------------------------*/ +int SetRefIndex(pSICSOBJ refl, int idx, double hkl[]) +{ + pHdb node = NULL; + hdbValue v; + pRLPriv priv = refl->pPrivate; + + + node = findReflection(refl, idx); + if (node != NULL) { + node = GetHipadabaNode(node, INAME); + assert(node); + v = MakeHdbFloatArray(priv->idxCount, hkl); + SetHipadabaPar(node, v, NULL); + return 1; + } else { + return 0; + } +} + +/*-----------------------------------------------------------------------------*/ +int SetRefAngles(pSICSOBJ refl, int idx, double ang[]) +{ + pHdb node = NULL; + hdbValue v; + pRLPriv priv = refl->pPrivate; + + + node = findReflection(refl, idx); + if (node != NULL) { + node = GetHipadabaNode(node, ANAME); + assert(node); + v = MakeHdbFloatArray(priv->angCount, ang); + SetHipadabaPar(node, v, NULL); + return 1; + } else { + return 0; + } +} + +/*-----------------------------------------------------------------------------*/ +int SetRefFlag(pSICSOBJ refl, int idx, int val) +{ + pHdb node = NULL; + hdbValue v; + pRLPriv priv = refl->pPrivate; + + + node = findReflection(refl, idx); + if (node != NULL) { + node = GetHipadabaNode(node, FLAG); + assert(node); + v = MakeHdbInt(val); + SetHipadabaPar(node, v, NULL); + return 1; + } else { + return 0; + } +} + +/*-------------------------------------------------------------------------------*/ +int GetRefIndex(pSICSOBJ refl, int idx, double hkl[]) +{ + pHdb node = NULL; + pRLPriv priv = refl->pPrivate; + + node = findReflection(refl, idx); + if (node != NULL) { + node = GetHipadabaNode(node, INAME); + assert(node); + memcpy(hkl, node->value.v.floatArray, priv->idxCount * sizeof(double)); + return 1; + } else { + return 0; + } +} + +/*-------------------------------------------------------------------------------*/ +int GetRefIndexID(pSICSOBJ refl, char *id, double hkl[]) +{ + pHdb node = NULL; + pRLPriv priv = refl->pPrivate; + char path[132]; + + snprintf(path, 132, "list/%s", id); + node = GetHipadabaNode(refl->objectNode, path); + if (node != NULL) { + node = GetHipadabaNode(node, INAME); + assert(node); + memcpy(hkl, node->value.v.floatArray, priv->idxCount * sizeof(double)); + return 1; + } else { + return 0; + } +} + +/*-------------------------------------------------------------------------------*/ +int GetRefAngles(pSICSOBJ refl, int idx, double ang[]) +{ + pHdb node = NULL; + pRLPriv priv = refl->pPrivate; + + + node = findReflection(refl, idx); + if (node != NULL) { + node = GetHipadabaNode(node, ANAME); + assert(node); + memcpy(ang, node->value.v.floatArray, priv->angCount * sizeof(double)); + return 1; + } else { + return 0; + } +} + +/*-------------------------------------------------------------------------------*/ +int GetRefAnglesID(pSICSOBJ refl, char *id, double hkl[]) +{ + pHdb node = NULL; + pRLPriv priv = refl->pPrivate; + char path[132]; + + snprintf(path, 132, "list/%s", id); + node = GetHipadabaNode(refl->objectNode, path); + if (node != NULL) { + node = GetHipadabaNode(node, ANAME); + assert(node); + memcpy(hkl, node->value.v.floatArray, priv->angCount * sizeof(double)); + return 1; + } else { + return 0; + } +} + +/*-------------------------------------------------------------------------------*/ +int GetRefFlag(pSICSOBJ refl, int idx) +{ + pHdb node = NULL; + pRLPriv priv = refl->pPrivate; + + + node = findReflection(refl, idx); + if (node != NULL) { + node = GetHipadabaNode(node, FLAG); + assert(node); + return node->value.v.intValue; + } else { + return -1; + } +} + +/*--------------------------------------------------------------------------------*/ +char *GetRefName(pSICSOBJ refl, int idx) +{ + pHdb node = NULL; + + node = findReflection(refl, idx); + if (node != NULL) { + return node->name; + } else { + return undef; + } +} diff --git a/reflist.h b/reflist.h new file mode 100644 index 00000000..ea591e7c --- /dev/null +++ b/reflist.h @@ -0,0 +1,60 @@ +/** + * Reflection list: a configurable list of reflections and + * their setting angles. To be used for four circle and possibly + * TAS diffraction. + * + * copyright: see File COPYRIGHT + * + * Mark Koennecke, July 2008 + */ +#ifndef REFLIST_H_ +#define REFLIST_H_ +#include +#include +/** + * This is an internal creation command: it only creates the data + * structure but does not add the commands + * \param name The name of the reflection list + * \return A SICSOBJ representing a new reflection list + */ +pSICSOBJ MakeReflectionListInt(char *name); +/** + * This is the standard SICS creation function for reflection lists + */ +int MakeReflectionList(SConnection * pCon, SicsInterp * pSics, + void *data, int argc, char *argv[]); +/** + * This creates a full reflection list with commands and adds it to the + * interpreter. + * \param pCon for error messages + * \param pSics the interpreter to add the command to + * \param The name of the reflection list + * \return A SICSOBJ representing the reflection list + */ +pSICSOBJ CreateReflectionList(SConnection * pCon, SicsInterp * pSics, + char *name); + +void ConfigureReflectionListIndex(pSICSOBJ refl, char *header); +void ConfigureReflectionListSettings(pSICSOBJ refl, char *header); + +void ClearReflectionList(pSICSOBJ refl); +int NewListReflection(pSICSOBJ refl); +void DelListReflection(pSICSOBJ refl, int idx); +int ReflectionListCount(pSICSOBJ refl); + +void AddRefIdx(pSICSOBJ refl, double hkl[]); +void AddRefIdxAng(pSICSOBJ refl, double hkl[], double ang[]); + +int SetRefIndex(pSICSOBJ refl, int idx, double hkl[]); +int SetRefAngles(pSICSOBJ refl, int idx, double ang[]); +int SetRefFlag(pSICSOBJ refl, int idx, int val); + +int GetRefIndex(pSICSOBJ refl, int idx, double hkl[]); +int GetRefIndexID(pSICSOBJ refl, char *id, double hkl[]); +int GetRefAngles(pSICSOBJ refl, int idx, double ang[]); +int GetRefAnglesID(pSICSOBJ refl, char *id, double ang[]); +int GetRefFlag(pSICSOBJ refl, int idx); + +char *GetRefName(pSICSOBJ refl, int idx); + +#endif /*REFLIST_H_ */ diff --git a/singlebi.c b/singlebi.c new file mode 100644 index 00000000..2ed0c732 --- /dev/null +++ b/singlebi.c @@ -0,0 +1,489 @@ +/** + * This is an implementation of the polymorphic single crystal calculation + * system defined in singlediff.h for a bisecting four circle diffractometer with + * an eulerian cradle. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, August 2008 + */ +#include +#include +#include +#include "singlediff.h" +#include "fourlib.h" +#include "ubfour.h" +#include "motor.h" +#include "singlex.h" +#include "motorlist.h" +#include "lld.h" + +/* + the tolerance in chi we give before we allow to fix omega with phi +*/ +#define CHITOLERANCE 3. +#define ABS(x) (x < 0 ? -(x) : (x)) +/*--------------------------------------------------------------------*/ +static int chiVertical(double chi) +{ + if (ABS(chi - .0) < CHITOLERANCE) { + return 1; + } + if (ABS(chi - 180.0) < CHITOLERANCE) { + return 1; + } + return 0; +} + +/*-----------------------------------------------------------------------*/ +static int checkTheta(pSingleDiff self, double *stt) +{ + char pError[132]; + int iTest; + float fHard; + pMotor pTheta; + + pTheta = SXGetMotor(TwoTheta); + if (pTheta == NULL) { + return 0; + } + + iTest = MotorCheckBoundary(pTheta, (float) *stt, &fHard, pError, 131); + if (!iTest) { + return -1; + } + return iTest; +} + +/*-----------------------------------------------------------------------*/ +static int checkBisecting(pSingleDiff self, + double *stt, double om, double chi, double phi) +{ + int iTest; + float fHard, fLimit; + char pError[132]; + pMotor pMot; + + /* check two theta */ + iTest = checkTheta(self, stt); + if (iTest != 1) { + return -1; + } + + /* for omega check against the limits +- SCANBORDER in order to allow for + a omega scan + */ + pMot = SXGetMotor(Omega); + if (pMot == NULL) { + return 0; + } + MotorGetPar(pMot, "softlowerlim", &fLimit); + if ((float) om < fLimit) { + iTest = 0; + } else { + iTest = 1; + MotorGetPar(pMot, "softupperlim", &fLimit); + if ((float) om > fLimit) { + iTest = 0; + } else { + iTest = 1; + } + } + + /* check chi and phi */ + pMot = SXGetMotor(Chi); + if (pMot == NULL) { + return 0; + } + iTest += MotorCheckBoundary(pMot, (float) chi, &fHard, pError, 131); + pMot = SXGetMotor(Phi); + if (pMot == NULL) { + return 0; + } + iTest += MotorCheckBoundary(pMot, (float) phi, &fHard, pError, 131); + if (iTest == 3) { /* none of them burns */ + return 1; + } + return 0; +} + +/*----------------------------------------------------------------------- + tryOmegaTweak tries to calculate a psi angle in order to put an + offending omega back into range. + + This routine also handles the special case when chi ~ 0 or chi ~ 180. + Then it is possible to fix a omega problem by turing in phi. + -----------------------------------------------------------------------*/ +static int tryOmegaTweak(pSingleDiff self, MATRIX z1, double *stt, + double *om, double *chi, double *phi) +{ + int status; + float fLower, fUpper, omTarget, omOffset, phiSign; + double dumstt, offom, offchi, offphi, lambda; + pMotor pOmega, pPhi; + + pOmega = SXGetMotor(Omega); + pPhi = SXGetMotor(Phi); + lambda = self->lambda; + if (pOmega == NULL) { + return 0; + } + + status = checkBisecting(self, stt, *om, *chi, *phi); + if (status < 0) { + return 0; /* stt is burning */ + } else if (status == 1) { + return 1; + } + + /* + Is omega really the problem? + */ + omTarget = -9999; + MotorGetPar(pOmega, "softlowerlim", &fLower); + MotorGetPar(pOmega, "softupperlim", &fUpper); + if (*om < fLower) { + omTarget = fLower; + } + if (*om > fUpper) { + omTarget = fUpper; + } + if (omTarget < -7000) { + return 0; + } + + /* + calculate omega offset + */ + omOffset = *om - omTarget; + omOffset = -omOffset; + + /* + check for the special case of chi == 0 or chi == 180 + */ + if (chiVertical(*chi)) { + dumstt = *stt; + offom = omTarget; + offchi = *chi; + MotorGetPar(pPhi, "sign", &phiSign); + offphi = *phi - omOffset * phiSign; + if (checkBisecting(self, &dumstt, offom, offchi, offphi)) { + *om = offom; + *chi = offchi; + *phi = offphi; + return 1; + } + } + + /* + calculate angles with omega offset + */ + status = z1ToAnglesWithOffset(lambda, z1, omOffset, &dumstt, + &offom, &offchi, &offphi); + if (!status) { + return 0; + } + + if (checkBisecting(self, &dumstt, offom, offchi, offphi)) { + *om = offom; + *chi = offchi; + *phi = offphi; + return 1; + } + return 0; +} + +/*----------------------------------------------------------------------*/ +static int hklInRange(void *data, double fSet[4], int mask[4]) +{ + pSingleDiff self = (pSingleDiff) data; + float fHard, fLimit; + char pError[132]; + int i, test; + double dTheta; + pMotor pOmega, pChi, pPhi; + + pOmega = SXGetMotor(Omega); + pChi = SXGetMotor(Chi); + pPhi = SXGetMotor(Phi); + if (pOmega == NULL || pChi == NULL || pPhi == NULL) { + return 0; + } + + /* 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(pOmega, "softlowerlim", &fLimit); + if ((float) fSet[1] < fLimit) { + mask[1] = 1; + } else { + mask[1] = 0; + MotorGetPar(pOmega, "softupperlim", &fLimit); + if ((float) fSet[1] > fLimit) { + mask[1] = 0; + } else { + mask[1] = 1; + } + } + + /* check chi and phi */ + mask[2] = MotorCheckBoundary(pChi, fSet[2], &fHard, pError, 131); + mask[3] = MotorCheckBoundary(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 calculateBiSettings(pSingleDiff self, + double *hkl, double *settings) +{ + MATRIX z1; + double stt, om, chi, phi, psi, ompsi, chipsi, phipsi, myPsi; + int i, test, mask[4], iRetry; + double lambda; + + myPsi = hkl[3]; + if (myPsi > 0.1) { + iRetry = 1; + } else { + iRetry = 699; + } + + z1 = calculateScatteringVector(self, hkl); + + /* + just the plain angle calculation + */ + if (!z1mToBisecting(self->lambda, z1, &stt, &om, &chi, &phi)) { + return 0; + } + + settings[0] = stt; + settings[1] = om; + settings[2] = chi; + settings[3] = phi; + if (iRetry == 1) { + rotatePsi(om, chi, phi, myPsi, &ompsi, &chipsi, &phipsi); + settings[1] = ompsi; + settings[2] = circlify(chipsi); + settings[3] = circlify(phipsi); + return 1; + } else { + if (hklInRange(self, settings, mask) == 1) { + return 1; + } else { + if (tryOmegaTweak(self, z1, &stt, &om, &chi, &phi) == 1) { + settings[0] = stt; + settings[1] = om; + settings[2] = chi; + settings[3] = phi; + return 1; + } else { + return findAllowedBisecting(self->lambda, z1, settings, hklInRange, + self); + } + } + } + return 0; +} + +/*-------------------------------------------------------------------*/ +static int settingsToBiList(struct __SingleDiff *self, double *settings) +{ + + setNewMotorTarget(self->motList, (char *) SXGetMotorName(TwoTheta), + (float) settings[0]); + setNewMotorTarget(self->motList, (char *) SXGetMotorName(Omega), + (float) settings[1]); + setNewMotorTarget(self->motList, (char *) SXGetMotorName(Chi), + (float) settings[2]); + setNewMotorTarget(self->motList, (char *) SXGetMotorName(Phi), + (float) settings[3]); + return 1; +} + +/*--------------------------------------------------------------------*/ +static int hklFromBiAngles(struct __SingleDiff *self, double *hkl) +{ + pIDrivable pDriv; + MATRIX UBinv, z1m, rez; + double z1[3], stt, om, chi, phi; + int i; + + pDriv = makeMotListInterface(); + pDriv->GetValue(&self->motList, pServ->dummyCon); + + UBinv = mat_inv(self->UB); + if (UBinv == NULL) { + return 0; + } + + stt = getListMotorPosition(self->motList, + (char *) SXGetMotorName(TwoTheta)); + om = getListMotorPosition(self->motList, (char *) SXGetMotorName(Omega)); + chi = getListMotorPosition(self->motList, (char *) SXGetMotorName(Chi)); + phi = getListMotorPosition(self->motList, (char *) SXGetMotorName(Phi)); + z1FromAngles(self->lambda, stt, om, chi, phi, z1); + z1m = vectorToMatrix(z1); + rez = mat_mul(UBinv, z1m); + for (i = 0; i < 3; i++) { + hkl[i] = rez[i][0]; + } + + mat_free(UBinv); + mat_free(z1m); + mat_free(rez); + return 1; +} + +/*--------------------------------------------------------------------*/ +static int hklFromBiAnglesGiven(struct __SingleDiff *self, + double *settings, double *hkl) +{ + MATRIX UBinv, z1m, rez; + double z1[3], stt, om, chi, phi; + int i; + + UBinv = mat_inv(self->UB); + if (UBinv == NULL) { + return 0; + } + + stt = settings[0]; + om = settings[1]; + chi = settings[2]; + phi = settings[3]; + + z1FromAngles(self->lambda, stt, om, chi, phi, z1); + z1m = vectorToMatrix(z1); + rez = mat_mul(UBinv, z1m); + for (i = 0; i < 3; i++) { + hkl[i] = rez[i][0]; + } + + mat_free(UBinv); + mat_free(z1m); + mat_free(rez); + return 1; +} + +/*------------------------------------------------------------------*/ +static int getBiReflection(char *id, reflection * r) +{ + pSICSOBJ refList; + double hkl[3], angles[4]; + + 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); + r->s2t = angles[0]; + r->om = angles[1]; + r->chi = angles[2]; + r->phi = angles[3]; + } + return 1; +} + +/*-------------------------------------------------------------------*/ +MATRIX calcBiUBFromTwo(pSingleDiff self, + char *refid1, char *refid2, int *err) +{ + MATRIX newUB; + reflection r1, r2; + lattice direct; + + direct.a = self->cell[0]; + direct.b = self->cell[1]; + direct.c = self->cell[2]; + direct.alpha = self->cell[3]; + direct.beta = self->cell[4]; + direct.gamma = self->cell[5]; + + if (!getBiReflection(refid1, &r1)) { + *err = REFERR; + return NULL; + } + if (!getBiReflection(refid2, &r2)) { + *err = REFERR; + return NULL; + } + + newUB = calcUBFromCellAndReflections(direct, r1, r2, err); + + return newUB; +} + +/*-------------------------------------------------------------------*/ +MATRIX calcBiUBFromThree(pSingleDiff self, + char *refid1, char *refid2, char *refid3, + int *err) +{ + MATRIX newUB; + reflection r1, r2, r3; + + if (!getBiReflection(refid1, &r1)) { + *err = REFERR; + return NULL; + } + if (!getBiReflection(refid2, &r2)) { + *err = REFERR; + return NULL; + } + if (!getBiReflection(refid3, &r3)) { + *err = REFERR; + return NULL; + } + + newUB = calcUBFromThreeReflections(r1, r2, r3, self->lambda, err); + return newUB; +} + +/*--------------------------------------------------------------------*/ +static int calcBiZ1(pSingleDiff self, char *refid, double z1[3]) +{ + reflection r1; + + if (!getBiReflection(refid, &r1)) { + return 0; + } + z1FromAngles(self->lambda, r1.s2t, r1.om, r1.chi, r1.phi, z1); + return 1; +} + +/*--------------------------------------------------------------------*/ +void initializeSingleBisecting(pSingleDiff diff) +{ + + if (diff->motList >= 0) { + LLDdelete(diff->motList); + } + diff->motList = LLDcreate(sizeof(MotControl)); + addMotorToList(diff->motList, (char *) SXGetMotorName(TwoTheta), .0); + addMotorToList(diff->motList, (char *) SXGetMotorName(Omega), .0); + addMotorToList(diff->motList, (char *) SXGetMotorName(Chi), .0); + addMotorToList(diff->motList, (char *) SXGetMotorName(Phi), .0); + + diff->calculateSettings = calculateBiSettings; + diff->settingsToList = settingsToBiList; + diff->hklFromAngles = hklFromBiAngles; + diff->hklFromAnglesGiven = hklFromBiAnglesGiven; + diff->calcUBFromTwo = calcBiUBFromTwo; + diff->calcUBFromThree = calcBiUBFromThree; + diff->calcZ1 = calcBiZ1; +} diff --git a/singlediff.c b/singlediff.c new file mode 100644 index 00000000..d9fb0318 --- /dev/null +++ b/singlediff.c @@ -0,0 +1,25 @@ +/** + * This is the implemetation of some common support functions for single + * crystal diffraction. + * + * copyright: see file COPYRIGHT + * + * Mark Koenencke, August 2008 +*/ +#include +#include "singlediff.h" +/*-----------------------------------------------------------------------*/ +MATRIX calculateScatteringVector(pSingleDiff pSingle, double fHKL[3]) +{ + MATRIX z1, hkl; + int i; + + hkl = mat_creat(3, 1, ZERO_MATRIX); + for (i = 0; i < 3; i++) { + hkl[i][0] = fHKL[i]; + } + z1 = mat_mul(pSingle->UB, hkl); + mat_free(hkl); + + return z1; +} diff --git a/singlediff.h b/singlediff.h new file mode 100644 index 00000000..26eea571 --- /dev/null +++ b/singlediff.h @@ -0,0 +1,91 @@ +/** + * This header file defines the polymorphic functions required for + * a single crystal diffractmeter, regardless of the geometry used. + * This exists in order to support the multiple modes of the TRICS + * instrument. But can be useful elsewhere. This code assumes that + * all motor driving happens via a motorlist as defined in + * motorlist.* + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, August 2008 + */ +#ifndef SINGLEDIFF_H_ +#define SINGLEDIFF_H_ +#include "matrix/matrix.h" + +#define REFERR -17001 + +/** + * A polymorphic structure for single crystal calculations. + */ +typedef struct __SingleDiff { + int motList; /** A list of associated real motors */ + void *userData; /** a pointer to a user data structure */ + MATRIX UB; /** The UB matrix */ + double lambda; /** The wavelength to use */ + double cell[6]; /** The unit cell */ + /** + * \brief calculates the settings for an input reciprocal space + * position hkl + * \param self A pointer to this data structure + * \param hkl The input reciprocal space position + * \param settings The angles calculated. + * \return 1 on success, 0 on error + */ + int (*calculateSettings) (struct __SingleDiff * self, + double *hkl, double *settings); + /** + * \brief copy setting angles into motor list + * \param self A pointer to this data structure + * \param settings The angles to prime the motor list with + * \return 1 on success, 0 on error + */ + int (*settingsToList) (struct __SingleDiff * self, double *settings); + /** + * \brief read all angles and convert to reciprocal space + * \param self A pointer to this data structure + * \param Result hkl values + * \return 1 on success, 0 on error + */ + int (*hklFromAngles) (struct __SingleDiff * self, double *hkl); + /** + * \brief Take angles given and convert to reciprocal space + * \param self A pointer to this data structure + * \param Result hkl values + * \return 1 on success, 0 on error + */ + int (*hklFromAnglesGiven) (struct __SingleDiff * self, double *settings, + double *hkl); + /** + * \brief calculate a UB matrix from two reflections (and the cell) + * \param self A pointer to this data structure + * \param refid1 first reflection ID + * \param refid2 second reflection ID + * \return NULL on error, a matrix else + */ + MATRIX(*calcUBFromTwo) (struct __SingleDiff * self, + char *refid1, char *refid2, int *err); + /** + * \brief calculate a UB matrix from three reflections + * \param self A pointer to this data structure + * \param refid1 first reflection ID + * \param refid2 second reflection ID + * \param refid3 third reflection ID + * \return NULL on error, a matrix else + */ + MATRIX(*calcUBFromThree) (struct __SingleDiff * self, + char *refid1, char *refid2, char *refid3, + int *err); + /** + * \brief calculate the scattering vector from the angles. + * \param refid The id of the reflection to calculate the + * scattering vector for + * \param z1 output scattering vector + * \return 1 on success, o on fail + */ + int (*calcZ1) (struct __SingleDiff * self, char *refid, double z1[3]); +} SingleDiff, *pSingleDiff; + +MATRIX calculateScatteringVector(pSingleDiff pSingle, double hkl[3]); +#endif /*SINGLEDIFF_H_ */ diff --git a/singlenb.c b/singlenb.c new file mode 100644 index 00000000..936851c1 --- /dev/null +++ b/singlenb.c @@ -0,0 +1,300 @@ +/** + * This is an implementation of the polymorphic single crystal calculation + * system defined in singlediff.h for a diffractometer in normal beam geometry. + * This means the detector tilts out of the instrument plane upwards or + * downwards. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, August 2008 + */ +#include +#include +#include +#include "singlediff.h" +#include "fourlib.h" +#include "ubfour.h" +#include "motor.h" +#include "singlex.h" +#include "motorlist.h" +#include "lld.h" + +/*-----------------------------------------------------------------------*/ +static int checkTheta(pSingleDiff self, double *stt) +{ + char pError[132]; + int iTest; + float fHard; + pMotor pTheta; + + pTheta = SXGetMotor(TwoTheta); + if (pTheta == NULL) { + return 0; + } + + iTest = MotorCheckBoundary(pTheta, (float) *stt, &fHard, pError, 131); + if (!iTest) { + return -1; + } + return iTest; +} + +/*-----------------------------------------------------------------------*/ +static int checkNormalBeam(double om, double *gamma, double nu, + double fSet[4], pSingleDiff self) +{ + int iTest; + char pError[132]; + float fHard; + pMotor pMot; + + fSet[0] = (float) *gamma; + fSet[1] = (float) om; + fSet[2] = (float) nu; + + /* check omega, gamma and nu */ + pMot = SXGetMotor(Omega); + if (pMot == NULL) { + return 0; + } + iTest = MotorCheckBoundary(pMot, (float) om, &fHard, pError, 131); + iTest += checkTheta(self, gamma); + pMot = SXGetMotor(Nu); + if (pMot == NULL) { + return 0; + } + iTest += MotorCheckBoundary(pMot, (float) nu, &fHard, pError, 131); + if (iTest == 3) { + return 1; + } + return 0; +} + +/*-------------------------------------------------------------------*/ +static int calculateNBSettings(pSingleDiff self, + double *hkl, double *settings) +{ + MATRIX z1; + double gamma, om, nu; + int status; + + z1 = calculateScatteringVector(self, hkl); + + status = z1mToNormalBeam(self->lambda, z1, &gamma, &om, &nu); + mat_free(z1); + if (status != 1) { + return 0; + } + if (checkNormalBeam(om, &gamma, nu, settings, self)) { + return 1; + } else { + if (checkNormalBeam(om + 360., &gamma, nu, settings, self)) { + return 1; + } else { + return 0; + } + } + return 0; +} + +/*-------------------------------------------------------------------*/ +static int settingsToNBList(struct __SingleDiff *self, double *settings) +{ + + setNewMotorTarget(self->motList, (char *) SXGetMotorName(TwoTheta), + (float) settings[0]); + setNewMotorTarget(self->motList, (char *) SXGetMotorName(Omega), + (float) settings[1]); + setNewMotorTarget(self->motList, (char *) SXGetMotorName(Nu), + (float) settings[2]); + return 1; +} + +/*--------------------------------------------------------------------*/ +static int hklFromNBAngles(struct __SingleDiff *self, double *hkl) +{ + pIDrivable pDriv; + MATRIX UBinv, z1m, rez; + double z1[3], stt, om, nu; + int i; + + pDriv = makeMotListInterface(); + pDriv->GetValue(&self->motList, pServ->dummyCon); + + UBinv = mat_inv(self->UB); + if (UBinv == NULL) { + return 0; + } + + stt = getListMotorPosition(self->motList, + (char *) SXGetMotorName(TwoTheta)); + om = getListMotorPosition(self->motList, (char *) SXGetMotorName(Omega)); + nu = getListMotorPosition(self->motList, (char *) SXGetMotorName(Nu)); + + z1FromNormalBeam(self->lambda, om, stt, nu, z1); + z1m = vectorToMatrix(z1); + rez = mat_mul(UBinv, z1m); + for (i = 0; i < 3; i++) { + hkl[i] = rez[i][0]; + } + + mat_free(UBinv); + mat_free(z1m); + mat_free(rez); + return 1; +} + +/*--------------------------------------------------------------------*/ +static int hklFromNBAnglesGiven(struct __SingleDiff *self, + double *settings, double *hkl) +{ + pIDrivable pDriv; + MATRIX UBinv, z1m, rez; + double z1[3], stt, om, nu; + int i; + + UBinv = mat_inv(self->UB); + if (UBinv == NULL) { + return 0; + } + + stt = settings[0]; + om = settings[1]; + nu = settings[2]; + + z1FromNormalBeam(self->lambda, om, stt, nu, z1); + z1m = vectorToMatrix(z1); + rez = mat_mul(UBinv, z1m); + for (i = 0; i < 3; i++) { + hkl[i] = rez[i][0]; + } + + mat_free(UBinv); + mat_free(z1m); + mat_free(rez); + return 1; +} + +/*----------------------------------------------------------------------------- + * For UB matrix calculations I use a trick. I do now know how to do + * that properly for normal beam geometry. But I know that the UB for + * bisecting and normal beam is the same. Thus I calculate the scattering + * vector z1 from the normal beam angles, then proceed to calculate + * bisecting angles from the scattering vector and use the bisecting + * geometry UB matrix calculation routines. + * -------------------------------------------------------------------------*/ +static int getNBReflection(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]; + z1FromNormalBeam(self->lambda, omnb, gamma, nu, z1); + if (!z1ToBisecting(self->lambda, z1, &stt, &om, &chi, &phi)) { + return 0; + } + r->s2t = stt; + r->om = om; + r->chi = chi; + r->phi = phi; + } + return 1; +} + +/*-------------------------------------------------------------------*/ +MATRIX calcNBUBFromTwo(pSingleDiff self, + char *refid1, char *refid2, int *err) +{ + MATRIX newUB; + reflection r1, r2; + lattice direct; + + direct.a = self->cell[0]; + direct.b = self->cell[1]; + direct.c = self->cell[2]; + direct.alpha = self->cell[3]; + direct.beta = self->cell[4]; + direct.gamma = self->cell[5]; + + if (!getNBReflection(self, refid1, &r1)) { + *err = REFERR; + return NULL; + } + if (!getNBReflection(self, refid2, &r2)) { + *err = REFERR; + return NULL; + } + + newUB = calcUBFromCellAndReflections(direct, r1, r2, err); + + return newUB; +} + +/*-------------------------------------------------------------------*/ +MATRIX calcNBFromThree(pSingleDiff self, + char *refid1, char *refid2, char *refid3, int *err) +{ + MATRIX newUB; + reflection r1, r2, r3; + + if (!getNBReflection(self, refid1, &r1)) { + *err = REFERR; + return NULL; + } + if (!getNBReflection(self, refid2, &r2)) { + *err = REFERR; + return NULL; + } + if (!getNBReflection(self, refid3, &r3)) { + *err = REFERR; + return NULL; + } + + newUB = calcUBFromThreeReflections(r1, r2, r3, self->lambda, err); + return newUB; +} + +/*--------------------------------------------------------------------*/ +static int calcNBZ1(pSingleDiff self, char *refid, double z1[3]) +{ + reflection r1; + + if (!getNBReflection(self, refid, &r1)) { + return 0; + } + z1FromAngles(self->lambda, r1.s2t, r1.om, r1.chi, r1.chi, z1); + return 1; +} + +/*--------------------------------------------------------------------*/ +void initializeNormalBeam(pSingleDiff diff) +{ + + if (diff->motList != 0) { + LLDdelete(diff->motList); + } + diff->motList = LLDcreate(sizeof(MotControl)); + addMotorToList(diff->motList, (char *) SXGetMotorName(TwoTheta), .0); + addMotorToList(diff->motList, (char *) SXGetMotorName(Omega), .0); + addMotorToList(diff->motList, (char *) SXGetMotorName(Nu), .0); + + diff->calculateSettings = calculateNBSettings; + diff->settingsToList = settingsToNBList; + diff->hklFromAngles = hklFromNBAngles; + diff->hklFromAnglesGiven = hklFromNBAnglesGiven; + diff->calcUBFromTwo = calcNBUBFromTwo; + diff->calcUBFromThree = calcNBFromThree; + diff->calcZ1 = calcNBZ1; +} diff --git a/singlenb.h b/singlenb.h new file mode 100644 index 00000000..f43c2de7 --- /dev/null +++ b/singlenb.h @@ -0,0 +1,16 @@ +/** + * This is an implementation of the polymorphic single crystal calculation + * system defined in singlediff.h for a diifractometer in normal beam geometry. + * This means the detector tilts out of the instrument plane upwards or + * downwards. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, August 2008 + */ +#ifndef SINGLENB_H_ +#define SINGLENB_H_ + +void initializeNormalBeam(pSingleDiff diff); + +#endif /*SINGLENB_H_ */ diff --git a/singlex.c b/singlex.c new file mode 100644 index 00000000..2d5ac844 --- /dev/null +++ b/singlex.c @@ -0,0 +1,868 @@ +/** + * This is the implementation file for the single crystal diffraction module within SICS. + * This is a master module which initializes all the other ones and holds + * data structures and entities particular to single crystal diffraction. + * + * The configuration of lambda is a bit tricky: it can either be a + * drivable (connected to a monochromator) or a SICS variable or we + * use an internal Hdb node. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, July 2008 + * + * Added SymRef command, Mark Koennecke, November 2008 + */ +#include "singlex.h" +#include "sicshipadaba.h" +#include "sicsvar.h" +#include "hkl.h" +#include "hklmot.h" +#include "ubcalc.h" +#include "fourmess.h" +#include "cone.h" +#include "singlediff.h" +#include "singelbi.h" +#include "singlenb.h" +#include "singletas.h" +#include "lld.h" + +extern void InitSimpleIndex(SConnection * pCon, SicsInterp * pSics); /* from simindex.c */ + +/* There can be only one singlex module in any instance of SICS */ +static pSICSOBJ singlex = NULL; +/*-------------------------------------------------------------------------*/ +typedef struct { + pMotor stt; + pMotor om; + pMotor chi; + pMotor phi; + pMotor nu; + pMotor sgu; + pMotor sgl; + pSICSOBJ refl; + double ub[9]; + double cell[6]; + SingleXModes mode; + pIDrivable lambdaDriv; + void *lambdaVar; + T_SgInfo *spaceGroup; + pSingleDiff diffractometer; +} SingleX, *pSingleX; +/*--------------------------------------------------------------------------*/ +static MotorFunction TextToFunc(char *txt) +{ + static char *map[] = { "stt", + "om", + "chi", + "phi", + "nu", + "sgl", + "sgu", + NULL, + }; + int count = 0; + while (map[count] != NULL) { + if (strcmp(map[count], txt) == 0) { + return count; + } + count++; + } + return -1; +} + +/*---------------------------------------------------------------------------*/ +static int ConfigureCmd(pSICSOBJ self, SConnection * pCon, + pHdb commandNode, pHdb par[], int nPar) +{ + pSingleX priv = self->pPrivate; + hdbValue v; + MotorFunction mf; + pMotor mot = NULL; + pDummy pDum = NULL; + CommandList *pCom = NULL; + + assert(priv != NULL); + if (nPar < 2) { + SCWrite(pCon, "ERROR: need two parameters to configure", eError); + return 0; + } + GetHipadabaPar(par[0], &v, pCon); + + if (strcmp(v.v.text, "lambda") == 0) { + GetHipadabaPar(par[1], &v, pCon); + pCom = FindCommand(pServ->pSics, v.v.text); + if (pCom != NULL) { + pDum = (pDummy) pCom->pData; + if ((priv->lambdaDriv = + pDum->pDescriptor->GetInterface(pDum, DRIVEID)) != NULL) { + priv->lambdaVar = pDum; + return 1; + } + if (strcmp(pDum->pDescriptor->name, "SicsVariable") == 0) { + priv->lambdaVar = pDum; + return 1; + } + } + SCWrite(pCon, "WARNING: creating local lambda node", eWarning); + AddSICSHdbPar(self->objectNode, "lambda", usUser, MakeHdbFloat(1.178)); + return 1; + } + + mf = TextToFunc(v.v.text); + if (mf < 0) { + SCPrintf(pCon, eError, + "ERROR: failed to map %s to configuration parameter", + v.v.text); + return 0; + } + + GetHipadabaPar(par[1], &v, pCon); + mot = (pMotor) FindCommandData(pServ->pSics, v.v.text, "Motor"); + if (mot == NULL) { + SCPrintf(pCon, eError, "ERROR: fialed to locate motor %s", v.v.text); + return 0; + } + switch (mf) { + case TwoTheta: + priv->stt = mot; + break; + case Omega: + priv->om = mot; + break; + case Chi: + priv->chi = mot; + break; + case Phi: + priv->phi = mot; + break; + case Nu: + priv->nu = mot; + break; + case Sgl: + priv->sgl = mot; + break; + case Sgu: + priv->sgu = mot; + break; + } + + return 1; +} + +/*---------------------------------------------------------------------------*/ +static int MotorCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + pSingleX priv = self->pPrivate; + MotorFunction mf; + pMotor pMot = NULL; + float fval; + + if (nPar < 1) { + SCWrite(pCon, "ERROR: missing argument to motor cmd ", eError); + return 0; + } + + mf = TextToFunc(par[0]->value.v.text); + if (mf < 0) { + SCPrintf(pCon, eError, "ERROR: %s is not a four circle motor", + par[0]->value.v.text); + return 0; + } + + pMot = SXGetMotor(mf); + if (pMot == NULL) { + SCPrintf(pCon, eError, "ERROR: %s motor is not configured", + par[0]->value.v.text); + return 0; + } + MotorGetSoftPosition(pMot, pCon, &fval); + SCPrintf(pCon, eValue, "%f", fval); + return 1; +} + +/*---------------------------------------------------------------------------*/ +static int MotorNamCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + pSingleX priv = self->pPrivate; + MotorFunction mf; + pMotor pMot = NULL; + + if (nPar < 1) { + SCWrite(pCon, "ERROR: missing argument to motor cmd ", eError); + return 0; + } + + mf = TextToFunc(par[0]->value.v.text); + if (mf < 0) { + SCPrintf(pCon, eError, "ERROR: %s is not a four circle motor", + par[0]->value.v.text); + return 0; + } + + pMot = SXGetMotor(mf); + if (pMot == NULL) { + SCPrintf(pCon, eError, "ERROR: %s motor is not configured", + par[0]->value.v.text); + return 0; + } + SCPrintf(pCon, eValue, "%s", pMot->name); + return 1; +} + +/*-------------------------------------------------------------------------- + * check the validity of the mode and set our internal field to + * the right value + * ------------------------------------------------------------------------*/ +extern char *trim(char *txt); +/*-------------------------------------------------------------------------*/ +static int findModeIndex(char *mode) +{ + static char *modes[] = { + "bi", + "nb", + "tas", + NULL + }; + int count = 0; + + while (modes[count] != NULL) { + if (strcmp(mode, modes[count]) == 0) { + return count; + } + count++; + } + return -1; +} + +/*-------------------------------------------------------------------------*/ +static hdbCallbackReturn SetModeCB(pHdb node, void *userData, + pHdbMessage mm) +{ + pHdbDataMessage set = NULL; + pSingleX priv = (pSingleX) singlex->pPrivate; + SConnection *pCon = NULL; + char *pTest; + int modeIdx; + + if ((set = GetHdbSetMessage(mm)) != NULL) { + pTest = trim(set->v->v.text); + pCon = set->callData; + modeIdx = findModeIndex(pTest); + switch (modeIdx) { + case Bisecting: + if (priv->chi == NULL || priv->om == NULL || priv->phi == NULL + || priv->stt == NULL) { + if (pCon != NULL) { + SCWrite(pCon, + "ERROR: required motor for bisecting not configured", + eError); + } + return hdbAbort; + } + priv->mode = Bisecting; + initializeSingleBisecting(priv->diffractometer); + return hdbContinue; + break; + case NB: + if (priv->nu == NULL || priv->om == NULL || priv->stt == NULL) { + if (pCon != NULL) { + SCWrite(pCon, + "ERROR: required motor for normal beam not configured", + eError); + } + return hdbAbort; + } + priv->mode = NB; + initializeNormalBeam(priv->diffractometer); + return hdbContinue; + break; + case Tas: + if (priv->sgu == NULL || priv->sgl == NULL || priv->om == NULL + || priv->stt == NULL) { + if (pCon != NULL) { + SCWrite(pCon, "ERROR: required motor for tas not configured", + eError); + } + return hdbAbort; + } + priv->mode = Tas; + initializeSingleTas(priv->diffractometer); + return hdbContinue; + default: + if (set->callData != NULL) { + pCon = set->callData; + SCPrintf(pCon, eError, "ERROR: %s is no valid singlex mode", + set->v->v.text); + } + return hdbAbort; + break; + } + } + return hdbContinue; +} + +/*-------------------------------------------------------------------------*/ +static hdbCallbackReturn GetLatticeCB(pHdb node, void *userData, + pHdbMessage mm) +{ + pHdbDataMessage get = NULL; + pSingleX priv = (pSingleX) singlex->pPrivate; + + + if ((get = GetHdbGetMessage(mm)) != NULL) { + node->value.v.intValue = priv->spaceGroup->XtalSystem; + return hdbContinue; + } + + return hdbContinue; +} + +/*------------------------------------------------------------------------------------*/ +static hdbCallbackReturn SetSpaceGroupCB(pHdb node, void *userData, + pHdbMessage mm) +{ + pHdbDataMessage set = NULL; + pSingleX priv = (pSingleX) singlex->pPrivate; + SConnection *pCon = NULL; + const T_TabSgName *tsgn = NULL; + const char *sgName; + + if ((set = GetHdbSetMessage(mm)) != NULL) { + pCon = set->callData; + tsgn = FindTabSgNameEntry(set->v->v.text, 'A'); + if (tsgn == NULL && set->callData != NULL) { + SCPrintf(pCon, eError, "ERROR: %s is no valid singlex spacegroup", + set->v->v.text); + return hdbAbort; + } + sgName = tsgn->HallSymbol; + InitSgInfo(priv->spaceGroup); + priv->spaceGroup->TabSgName = tsgn; + ParseHallSymbol(sgName, priv->spaceGroup); + if (SgError != NULL) { + if (pCon != NULL) { + SCPrintf(pCon, eError, "ERROR: %s parsing space group symbol %s", + SgError, set->v->v.text); + } + SgError = NULL; + return hdbAbort; + } + CompleteSgInfo(priv->spaceGroup); + return hdbContinue; + } + return hdbContinue; +} + +/*-------------------------------------------------------------------------*/ +static hdbCallbackReturn SetUBCB(pHdb node, void *userData, pHdbMessage mm) +{ + pHdbDataMessage set = NULL; + SConnection *pCon = NULL; + hdbValue v; + pHdb oldUB; + + if ((set = GetHdbSetMessage(mm)) != NULL) { + oldUB = (pHdb) userData; + assert(oldUB != NULL); + GetHipadabaPar(node, &v, set->callData); + SetHipadabaPar(oldUB, v, set->callData); + return hdbContinue; + } + return hdbContinue; +} + +/*-----------------------------------------------------------------------------*/ +static int RecoverUBCmd(pSICSOBJ self, SConnection * pCon, + pHdb commandNode, pHdb par[], int nPar) +{ + pHdb oldub, ub; + hdbValue v; + + oldub = GetHipadabaNode(self->objectNode, "oldub"); + ub = GetHipadabaNode(self->objectNode, "ub"); + GetHipadabaPar(oldub, &v, pCon); + UpdateHipadabaPar(ub, v, pCon); + SCSendOK(pCon); + return 1; +} + +/*-----------------------------------------------------------------------------*/ +static int SymRefCmd(pSICSOBJ self, SConnection * pCon, pHdb commandNode, + pHdb par[], int nPar) +{ + double hkl[3], settings[6]; + char buffer[132]; + pSingleX priv = (pSingleX) singlex->pPrivate; + T_Eq_hkl equiv; + int i; + + + if (nPar < 3) { + SCWrite(pCon, "ERROR: need start h,k,l for symref ", eError); + return 0; + } + hkl[0] = par[0]->value.v.intValue; + hkl[1] = par[1]->value.v.intValue; + hkl[2] = par[2]->value.v.intValue; + + BuildEq_hkl(priv->spaceGroup, &equiv, (int) hkl[0], (int) hkl[1], + (int) hkl[2]); + for (i = 0; i < equiv.N; i++) { + hkl[0] = equiv.h[i]; + hkl[1] = equiv.k[i]; + hkl[2] = equiv.l[i]; + if (priv->diffractometer-> + calculateSettings(priv->diffractometer, hkl, settings) == 1) { + snprintf(buffer, 512, "%d,%d,%d", equiv.h[i], equiv.k[i], + equiv.l[i]); + SCWrite(pCon, buffer, eValue); + return 1; + } + } + SCWrite(pCon, "ERROR: no reachable symmetry equivalent found", eError); + return 0; +} + +/*-----------------------------------------------------------------------------*/ +static int SingleXAction(SConnection * pCon, SicsInterp * pSics, + void *data, int argc, char *argv[]) +{ + pSICSOBJ self = data; + pSingleX priv = self->pPrivate; + double lambda; + + assert(self != NULL); + assert(priv != NULL); + + if (argc > 1) { + if (strcmp(argv[1], "lambda") == 0 && priv->lambdaVar != NULL) { + if (argc < 3) { + lambda = SXGetLambda(); + SCPrintf(pCon, eValue, "singlex.lambda = %f", lambda); + return 1; + } else { + if (priv->lambdaDriv == NULL) { + VarSetFloat((pSicsVariable) priv->lambdaVar, atof(argv[2]), + SCGetRights(pCon)); + return 1; + } else { + priv->lambdaDriv->SetValue(priv->lambdaVar, pCon, atof(argv[2])); + return 1; + } + } + } + } + return InterInvokeSICSOBJ(pCon, pSics, data, argc, argv); +} + +/*---------------------------------------------------------------------------*/ +static void SinglexKill(void *data) +{ + pSingleX priv = data; + + if (priv == NULL) { + return; + } + + if (priv->spaceGroup->ListSeitzMx != NULL) { + free(priv->spaceGroup->ListSeitzMx); + } + if (priv->spaceGroup->ListRotMxInfo != NULL) { + free(priv->spaceGroup->ListRotMxInfo); + } + if (priv->spaceGroup != NULL) { + free(priv->spaceGroup); + } + if (priv->diffractometer != NULL) { + if (priv->diffractometer->motList > 0) { + LLDdelete(priv->diffractometer->motList); + } + if (priv->diffractometer->UB != NULL) { + mat_free(priv->diffractometer->UB); + } + } + free(priv); +} + +/*---------------------------------------------------------------------------*/ +int MakeSingleX(SConnection * pCon, SicsInterp * pSics, + void *data, int argc, char *argv[]) +{ + pSICSOBJ pNew = NULL; + pSingleX priv = NULL; + pHdb cmd = NULL, par = NULL, oldub; + + pNew = MakeSICSOBJ("singlex", "SingleCrystalDiffraction"); + priv = calloc(1, sizeof(SingleX)); + if (pNew == NULL || priv == NULL) { + SCWrite(pCon, "ERROR: out of memory creating singlex", eError); + return 0; + } + priv->refl = CreateReflectionList(pCon, pSics, "ref"); + priv->diffractometer = calloc(1, sizeof(SingleDiff)); + priv->mode = Bisecting; + if (priv->refl == NULL || priv->diffractometer == NULL) { + SCWrite(pCon, "ERROR: out of memory creating singlex", eError); + return 0; + } + priv->diffractometer->motList = -1; + priv->diffractometer->UB = mat_creat(3, 3, UNIT_MATRIX); + + priv->spaceGroup = calloc(1, sizeof(T_SgInfo)); + if (priv->spaceGroup == NULL) { + return 0; + } + priv->spaceGroup->MaxList = 192; + priv->spaceGroup->ListSeitzMx = + malloc(192 * sizeof(*priv->spaceGroup->ListSeitzMx)); + priv->spaceGroup->ListRotMxInfo = + malloc(192 * sizeof(*priv->spaceGroup->ListRotMxInfo)); + if (priv->spaceGroup->ListSeitzMx == NULL + || priv->spaceGroup->ListRotMxInfo == NULL) { + return 0; + } + InitSgInfo(priv->spaceGroup); + singlex = pNew; + + pNew->pPrivate = priv; + pNew->KillPrivate = SinglexKill; + + cmd = + AddSICSHdbPar(pNew->objectNode, "configure", usMugger, + MakeSICSFunc(ConfigureCmd)); + AddSICSHdbPar(cmd, "device", usMugger, MakeHdbText("")); + AddSICSHdbPar(cmd, "name", usMugger, MakeHdbText("")); + + cmd = + AddSICSHdbPar(pNew->objectNode, "motval", usSpy, + MakeSICSFunc(MotorCmd)); + AddSICSHdbPar(cmd, "name", usMugger, MakeHdbText("")); + + cmd = + AddSICSHdbPar(pNew->objectNode, "motnam", usSpy, + MakeSICSFunc(MotorNamCmd)); + AddSICSHdbPar(cmd, "name", usMugger, MakeHdbText("")); + + cmd = + AddSICSHdbPar(pNew->objectNode, "recoverub", usUser, + MakeSICSFunc(RecoverUBCmd)); + + cmd = + AddSICSHdbPar(pNew->objectNode, "cell", usUser, + makeHdbValue(HIPFLOATAR, 6)); + SetHdbProperty(cmd, "__save", "true"); + + cmd = + AddSICSHdbPar(pNew->objectNode, "lattice", usInternal, + MakeHdbInt(0)); + AppendHipadabaCallback(cmd, + MakeHipadabaCallback(GetLatticeCB, NULL, NULL)); + + + cmd = + AddSICSHdbPar(pNew->objectNode, "oldub", usUser, + makeHdbValue(HIPFLOATAR, 9)); + SetHdbProperty(cmd, "__save", "true"); + oldub = cmd; + cmd = + AddSICSHdbPar(pNew->objectNode, "ub", usUser, + makeHdbValue(HIPFLOATAR, 9)); + SetHdbProperty(cmd, "__save", "true"); + PrependHipadabaCallback(cmd, MakeHipadabaCallback(SetUBCB, oldub, NULL)); + + cmd = + AddSICSHdbPar(pNew->objectNode, "planenormal", usUser, + makeHdbValue(HIPFLOATAR, 3)); + SetHdbProperty(cmd, "__save", "true"); + + cmd = AddSICSHdbPar(pNew->objectNode, "mode", usUser, MakeHdbText("bi")); + SetHdbProperty(cmd, "__save", "true"); + PrependHipadabaCallback(cmd, + MakeHipadabaCallback(SetModeCB, NULL, NULL)); + + cmd = + AddSICSHdbPar(pNew->objectNode, "spacegroup", usUser, + MakeHdbText("P")); + SetHdbProperty(cmd, "__save", "true"); + PrependHipadabaCallback(cmd, + MakeHipadabaCallback(SetSpaceGroupCB, NULL, + NULL)); + SetHipadabaPar(cmd, MakeHdbText("P"), pCon); + + cmd = + AddSICSHdbPar(pNew->objectNode, "symref", usUser, + MakeSICSFunc(SymRefCmd)); + AddSICSHdbPar(cmd, "h", usUser, MakeHdbInt(0)); + AddSICSHdbPar(cmd, "k", usUser, MakeHdbInt(0)); + AddSICSHdbPar(cmd, "l", usUser, MakeHdbInt(0)); + + cmd = + AddSICSHdbPar(pNew->objectNode, "peaksearch", usUser, + makeHdbValue(HIPNONE, 1)); + SetHdbProperty(cmd, "__save", "true"); + par = AddSICSHdbPar(cmd, "min2t", usUser, MakeHdbFloat(5.)); + SetHdbProperty(par, "__save", "true"); + par = AddSICSHdbPar(cmd, "step2t", usUser, MakeHdbFloat(1.)); + SetHdbProperty(par, "__save", "true"); + par = AddSICSHdbPar(cmd, "max2t", usUser, MakeHdbFloat(15.)); + SetHdbProperty(par, "__save", "true"); + par = AddSICSHdbPar(cmd, "stepchi", usUser, MakeHdbFloat(2.)); + SetHdbProperty(par, "__save", "true"); + par = AddSICSHdbPar(cmd, "stepphi", usUser, MakeHdbFloat(.5)); + SetHdbProperty(par, "__save", "true"); + par = AddSICSHdbPar(cmd, "stepom", usUser, MakeHdbFloat(.5)); + SetHdbProperty(par, "__save", "true"); + par = AddSICSHdbPar(cmd, "stepnu", usUser, MakeHdbFloat(.5)); + SetHdbProperty(par, "__save", "true"); + par = AddSICSHdbPar(cmd, "phimin", usUser, MakeHdbFloat(.0)); + SetHdbProperty(par, "__save", "true"); + par = AddSICSHdbPar(cmd, "phimax", usUser, MakeHdbFloat(180.)); + SetHdbProperty(par, "__save", "true"); + par = AddSICSHdbPar(cmd, "chimin", usUser, MakeHdbFloat(90.)); + SetHdbProperty(par, "__save", "true"); + par = AddSICSHdbPar(cmd, "chimax", usUser, MakeHdbFloat(180.)); + SetHdbProperty(par, "__save", "true"); + + if (AddCommand(pSics, "singlex", SingleXAction, KillSICSOBJ, pNew)) { + HKLFactory(pCon, pSics, data, argc, argv); + HKLMotInstall(pCon, pSics, data, argc, argv); + CreateUBCalc(pCon, pSics, "ubcalc", "hkl"); + InstallFourMess(pCon, pSics); + MakeCone(pCon, pSics, data, argc, argv); + InitSimpleIndex(pCon, pSics); + } else { + SCWrite(pCon, "ERROR: duplicate initialization of singlex module", + eError); + return 0; + } + return 1; +} + +/*---------------------------------------------------------------------------*/ +pMotor SXGetMotor(MotorFunction m) +{ + pSingleX priv = (pSingleX) singlex->pPrivate; + switch (m) { + case TwoTheta: + return priv->stt; + break; + case Omega: + return priv->om; + break; + case Chi: + return priv->chi; + break; + case Phi: + return priv->phi; + break; + case Nu: + return priv->nu; + break; + case Sgl: + return priv->sgl; + break; + case Sgu: + return priv->sgu; + break; + default: + return NULL; + } + return NULL; +} + +/*---------------------------------------------------------------------------*/ +const char *SXGetMotorName(MotorFunction m) +{ + pMotor mot = SXGetMotor(m); + + if (mot == NULL) { + return NULL; + } else { + return (const char *) mot->name; + } +} + +/*-------------------------------------------------------------------------- + * The technique used here is a bit funny. I do not wish to return something + * from a hdbValue data structure because either because it will be deleted + * in some stage or I introduce a memory leak. Thus I clean up everything here, + * make a copy in our private data structure and return that. Because that + * data structure will persist to the end of the program. Same for cell. + * -------------------------------------------------------------------------*/ +const double *SXGetUB() +{ + hdbValue v; + pHdb node = NULL; + pSingleX priv = (pSingleX) singlex->pPrivate; + + node = GetHipadabaNode(singlex->objectNode, "ub"); + assert(node != NULL); + GetHipadabaPar(node, &v, NULL); + memcpy(priv->ub, v.v.floatArray, 9 * sizeof(double)); + ReleaseHdbValue(&v); + return priv->ub; +} + +/*--------------------------------------------------------------------------*/ +void SXSetUB(double ub[9]) +{ + pHdb node = NULL; + pSingleX priv = (pSingleX) singlex->pPrivate; + + node = GetHipadabaNode(singlex->objectNode, "ub"); + assert(node != NULL); + SetHipadabaPar(node, MakeHdbFloatArray(9, ub), NULL); +} + +/*------------------------------------------------------------------------*/ +double *SXGetPlanenormal() +{ + hdbValue v; + pHdb node = NULL; + pSingleX priv = (pSingleX) singlex->pPrivate; + double *normal = NULL; + + node = GetHipadabaNode(singlex->objectNode, "planenormal"); + assert(node != NULL); + GetHipadabaPar(node, &v, NULL); + normal = malloc(3 * sizeof(double)); + if (normal == NULL) { + return NULL; + } + memcpy(normal, v.v.floatArray, 3 * sizeof(double)); + ReleaseHdbValue(&v); + return normal; +} + +/*--------------------------------------------------------------------------*/ +void SXSetPlanenormal(double normal[3]) +{ + pHdb node = NULL; + pSingleX priv = (pSingleX) singlex->pPrivate; + + node = GetHipadabaNode(singlex->objectNode, "planenormal"); + assert(node != NULL); + SetHipadabaPar(node, MakeHdbFloatArray(3, normal), NULL); +} + +/*--------------------------------------------------------------------------*/ +const double *SXGetCell() +{ + hdbValue v; + pHdb node = NULL; + pSingleX priv = (pSingleX) singlex->pPrivate; + + node = GetHipadabaNode(singlex->objectNode, "cell"); + assert(node != NULL); + GetHipadabaPar(node, &v, NULL); + memcpy(priv->cell, v.v.floatArray, 6 * sizeof(double)); + ReleaseHdbValue(&v); + return priv->cell; +} + +/*--------------------------------------------------------------------------*/ +void SXSetCell(double cell[6]) +{ + pHdb node = NULL; + pSingleX priv = (pSingleX) singlex->pPrivate; + + node = GetHipadabaNode(singlex->objectNode, "cell"); + assert(node != NULL); + SetHipadabaPar(node, MakeHdbFloatArray(6, cell), NULL); +} + +/*---------------------------------------------------------------------------*/ +double SXGetLambda() +{ + pSingleX priv = (pSingleX) singlex->pPrivate; + pSicsVariable lam = NULL; + float val; + hdbValue v; + pHdb node = NULL; + + if (priv->lambdaDriv != NULL) { + return priv->lambdaDriv->GetValue(priv->lambdaVar, pServ->dummyCon); + } else if (priv->lambdaVar != NULL) { + lam = (pSicsVariable) priv->lambdaVar; + VarGetFloat(lam, &val); + return val; + } + node = GetHipadabaNode(singlex->objectNode, "lambda"); + assert(node != NULL); + GetHipadabaPar(node, &v, NULL); + return v.v.doubleValue; +} + +/*---------------------------------------------------------------------------*/ +pSICSOBJ SXGetReflectionList() +{ + pSingleX priv = (pSingleX) singlex->pPrivate; + return priv->refl; +} + +/*---------------------------------------------------------------------------*/ +SingleXModes SXGetMode() +{ + pSingleX priv = (pSingleX) singlex->pPrivate; + return priv->mode; +} + +/*---------------------------------------------------------------------------*/ +T_SgInfo *SXGetSpaceGroup() +{ + pSingleX priv = (pSingleX) singlex->pPrivate; + return priv->spaceGroup; +} + +/*---------------------------------------------------------------------------*/ +pSingleDiff SXGetDiffractometer() +{ + pSingleX priv = (pSingleX) singlex->pPrivate; + const double *cell, *ub; + double lambda; + int i; + + /* Not initialized. + * Trouble is motors need to be configured first, the we can set a diffractometer. + * If this is ever triggered, the most likely cause is that no mode was selected in the + * initialisation file. + */ + if (priv->diffractometer->calculateSettings == NULL) { + return NULL; + } + + cell = SXGetCell(); + memcpy(priv->diffractometer->cell, cell, 6 * sizeof(double)); + priv->diffractometer->lambda = SXGetLambda(); + ub = SXGetUB(); + for (i = 0; i < 3; i++) { + priv->diffractometer->UB[0][i] = ub[i]; + priv->diffractometer->UB[1][i] = ub[i + 3]; + priv->diffractometer->UB[2][i] = ub[i + 6]; + } + return priv->diffractometer; +} + +/*---------------------------------------------------------------------------*/ +void SXSetMode(SingleXModes m) +{ + pHdb node = NULL; + + node = GetHipadabaNode(singlex->objectNode, "mode"); + assert(node != NULL); + + switch (m) { + case Bisecting: + SetHipadabaPar(node, MakeHdbText(strdup("bi")), NULL); + break; + case NB: + SetHipadabaPar(node, MakeHdbText(strdup("nb")), NULL); + break; + case Tas: + SetHipadabaPar(node, MakeHdbText(strdup("tas")), NULL); + break; + } +} diff --git a/singlex.h b/singlex.h new file mode 100644 index 00000000..30646319 --- /dev/null +++ b/singlex.h @@ -0,0 +1,52 @@ +/** + * This is the header file for the single crystal diffraction module within SICS. + * This is a master module which initializes all the other ones and holds + * data structures and entities particular to single crystal diffraction. + * + * copyright: see file COPYRIGHT + * + * Mark Koennecke, July 2008 + */ +#ifndef SINGLEX_H_ +#define SINGLEX_H_ +#include +#include +#include "reflist.h" +#include "sginfo.h" +#include "singlediff.h" + +typedef enum { + TwoTheta, Omega, Chi, Phi, Nu, Sgu, Sgl +} MotorFunction; +pMotor SXGetMotor(MotorFunction m); +const char *SXGetMotorName(MotorFunction m); + +const double *SXGetUB(); +void SXSetUB(double ub[9]); + +double *SXGetPlanenormal(); +void SXSetPlanenormal(double planeNormal[3]); + +const double *SXGetCell(); +void SXSetCell(double c[6]); + +double SXGetLambda(); + +pSICSOBJ SXGetReflectionList(); + +typedef enum { + Bisecting, NB, Tas +} SingleXModes; + +SingleXModes SXGetMode(); +void SXSetMode(SingleXModes m); + +pSingleDiff SXGetDiffractometer(); + +T_SgInfo *SXGetSpaceGroup(); + + +int MakeSingleX(SConnection * pCon, SicsInterp * pSics, + void *data, int argc, char *argv[]); + +#endif /*SINGLEX_H_ */