From f9121c1afe057876a9b3cced49fd1315d03ce728 Mon Sep 17 00:00:00 2001 From: koennecke Date: Tue, 17 Feb 2009 08:40:55 +0000 Subject: [PATCH] - Added missing files --- reflist.c | 878 +++++++++++++++++++++++++++++++++++++++++++++++++++ reflist.h | 60 ++++ singlebi.c | 489 ++++++++++++++++++++++++++++ singlediff.c | 25 ++ singlediff.h | 91 ++++++ singlenb.c | 300 ++++++++++++++++++ singlenb.h | 16 + singlex.c | 868 ++++++++++++++++++++++++++++++++++++++++++++++++++ singlex.h | 52 +++ 9 files changed, 2779 insertions(+) create mode 100644 reflist.c create mode 100644 reflist.h create mode 100644 singlebi.c create mode 100644 singlediff.c create mode 100644 singlediff.h create mode 100644 singlenb.c create mode 100644 singlenb.h create mode 100644 singlex.c create mode 100644 singlex.h 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_ */