linse/sics
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
main
sics/fourmess.c
Koennecke Mark 96114fc803 Added the generation of incomensurate reflections for 0,0,0 in fourmess 
Applied Douglas  task deletion patch
Fixed a bug with repeated motor failures in interface.c
2016-04-13 09:18:13 +02:00

1254 lines
37 KiB
C
Raw Permalink Blame History

/**
* This is a new version of the four circle reflection measurement module. Its task is
* to provide functionality required for measuring a list of reflections. This is done
* with the help of some scripts binding things together. This module provides:
* - a table of scan parameters for measuring reflections at different angles.
* - the functions required to produce the scan data files for single detector
* mode.
* - a list of reflections to measure.
*
* In contrast to the old mesure module, this only provides some help to
* scripts. Thus, through scripts, a more flexible measurement is possible.
*
* copyright: see file COPYRIGHT
*
* Mark Koennecke, July 2008
*/
#include <stdio.h>
#include <assert.h>
#include "singlex.h"
#include "sicsobj.h"
#include "fourtable.h"
#include "sicshipadaba.h"
#include "sicsvar.h"
#include <sics.h>
#include "scan.h"
#include "stdscan.h"
#include "evcontroller.h"
#include "integrate.h"
#include "sginfo.h"
#include "matrix/matrix.h"
#include "cell.h"
#include "fourlib.h"
#include "counter.h"
#include "scanvar.h"
#include "scan.i"
#include "sdynar.h"
#include "lld.h"
extern char *trim(char *);
extern void SNXFormatTime(char *pBueffel, int iLen);
#define ABS(x) (x < 0 ? -(x) : (x))
/*---------------------------------------------------------------------------*/
typedef struct {
FILE *profFile; /* file with reflection profiles, ccl */
FILE *hklFile; /* file with integrated intensities */
pSICSOBJ stepTable; /* table with the scan parameters */
char *currentFileRoot;
pSICSOBJ messList;
pHdb currentRefl; /* the current reflection being measured */
int count;
pScanData pScanner;
int masterCount; /* the number of master reflection as craeted by indgen */
} FourMess, *pFourMess;
/*---------------------------------------------------------------------------*/
static void KillFourMess(void *data)
{
pFourMess priv = (pFourMess) data;
if (priv == NULL) {
return;
}
if (priv->profFile != NULL) {
fclose(priv->profFile);
}
if (priv->hklFile != NULL) {
fclose(priv->hklFile);
}
if (priv->currentFileRoot != NULL) {
free(priv->currentFileRoot);
}
if (priv->stepTable >= 0) {
DeleteFourCircleTable(priv->stepTable);
}
free(priv);
}
/*----------------------------------------------------------------------------*/
static int FourMessAction(SConnection * pCon, SicsInterp * pSics,
void *data, int argc, char *argv[])
{
pFourMess priv = NULL;
pSICSOBJ self = (pSICSOBJ) data;
int err, status;
assert(self != NULL);
priv = self->pPrivate;
assert(priv != NULL);
if (argc < 2) {
SCWrite(pCon, "ERROR: insufficient number of arguments to fmess",
eError);
return 0;
}
if (strcmp(argv[1], "table") == 0) {
return HandleFourCircleCommands(priv->stepTable, pCon,
argc, argv, &err);
}
status = InvokeSICSOBJ(pCon, pSics, data, argc, argv);
if (status < 0) {
SCPrintf(pCon, eError, "ERROR: %s no subcommand or parameter",
argv[1]);
return 0;
}
return status;
}
/*-----------------------------------------------------------------------------*/
static int FourMessClose(pSICSOBJ self, SConnection * pCon,
pHdb commandNode, pHdb par[], int nPar)
{
pFourMess priv = NULL;
priv = self->pPrivate;
if (priv->hklFile != NULL) {
fclose(priv->hklFile);
priv->hklFile = NULL;
}
if (priv->profFile != NULL) {
fclose(priv->profFile);
priv->profFile = NULL;
}
return 1;
}
/*-----------------------------------------------------------------------------*/
static int FourMessStart(pSICSOBJ self, SConnection * pCon,
pHdb commandNode, pHdb par[], int nPar)
{
pFourMess priv = NULL;
char pFilename[1024], pRoot[512];
char pBueffel[1024];
double lambda;
const double *dUB;
pSicsVariable pVar = NULL;
char *pFile = NULL, *pPtr;
FILE *temp = NULL;
pHdb node;
hdbValue v;
assert(self);
assert(pCon);
priv = self->pPrivate;
if (nPar < 1) {
SCWrite(pCon, "ERROR: need file name parameter to start", eError);
return 0;
}
/* close open files if so */
if (priv->hklFile != NULL) {
FourMessClose(self, pCon, commandNode, par, nPar);
}
/* create filename root */
GetHipadabaPar(par[0], &v, pCon);
pFile = strdup(v.v.text);
pPtr = strrchr(pFile, (int) '.');
pPtr++;
*pPtr = '\0';
if (priv->currentFileRoot != NULL) {
free(priv->currentFileRoot);
}
priv->currentFileRoot = strdup(pFile);
free(pFile);
strlcpy(pRoot, priv->currentFileRoot, 511);
/* open the reflection file */
sprintf(pBueffel, "Writing to %s.ccl, .rfl", pRoot);
SCWrite(pCon, pBueffel, eLog);
strcpy(pFilename, pRoot);
strcat(pFilename, "ccl");
priv->profFile = fopen(pFilename, "w");
if (!priv->profFile) {
sprintf(pBueffel, "ERROR: SERIOUS TROUBLE: cannot open %s!",
pFilename);
SCWrite(pCon, pBueffel, eError);
SCSetInterrupt(pCon, eAbortBatch);
return 0;
}
node = GetHipadabaNode(self->objectNode, "template");
assert(node != NULL);
GetHipadabaPar(node, &v, pCon);
temp = fopen(v.v.text, "r");
if (temp == NULL) {
SCWrite(pCon, "ERROR: failed to open header template", eError);
}
if (temp != NULL) {
WriteTemplate(priv->profFile, temp, pFilename, NULL,
pCon, pServ->pSics);
fclose(temp);
}
/* open hkl-data file */
strcpy(pFilename, pRoot);
strcat(pFilename, "rfl");
priv->hklFile = fopen(pFilename, "w");
if (!priv->hklFile) {
sprintf(pBueffel, "ERROR: SERIOUS TROUBLE: cannot open %s!",
pFilename);
SCWrite(pCon, pBueffel, eError);
SCSetInterrupt(pCon, eAbortBatch);
return 0;
}
fputs(pFilename, priv->hklFile);
fputs("\n", priv->hklFile);
/* write some header data */
SNXFormatTime(pBueffel, 1024);
fprintf(priv->hklFile, "filetime = %s\n", pBueffel);
lambda = SXGetLambda();
fprintf(priv->hklFile, "lambda = %f Angstroem\n", lambda);
dUB = SXGetUB();
fprintf(priv->hklFile,
"UB = %7.6f %7.6f %7.6f %7.6f %7.6f %7.6f %7.6f %7.6f %7.6f\n",
dUB[0], dUB[1], dUB[2], dUB[3], dUB[4], dUB[5], dUB[6], dUB[7],
dUB[8]);
/* write sample & user info */
strcpy(pBueffel, "CCL, Instr=TRICS, ");
pVar = FindVariable(pServ->pSics, "sample");
if (pVar) {
fprintf(priv->hklFile, "sample = %s\n", pVar->text);
}
pVar = FindVariable(pServ->pSics, "user");
if (pVar) {
fprintf(priv->hklFile, "user = %s \n", pVar->text);
}
priv->count = 0;
return 1;
}
/*-----------------------------------------------------------------------------*/
static int FourMessReopen(pSICSOBJ self, SConnection * pCon,
pHdb commandNode, pHdb par[], int nPar)
{
pFourMess priv = NULL;
char pFilename[1024];
priv = self->pPrivate;
if(priv->hklFile != NULL){
return 1; /* still open! */
}
if(priv->currentFileRoot == NULL){
SCWrite(pCon,"ERROR: nothing to append too", eLogError);
return 0;
}
strcpy(pFilename, priv->currentFileRoot);
strcat(pFilename, "ccl");
priv->profFile = fopen(pFilename, "a");
if(priv->profFile == NULL){
SCPrintf(pCon,eLogError,"ERROR: failed to reopen %s", pFilename);
return 0;
}
strcpy(pFilename, priv->currentFileRoot);
strcat(pFilename, "rfl");
priv->hklFile = fopen(pFilename, "a");
if(priv->hklFile == NULL){
SCPrintf(pCon,eLogError,"ERROR: failed to reopen %s", pFilename);
return 0;
}
SCPrintf(pCon,eValue,"Reopened %s.ccl,.rfl", priv->currentFileRoot);
return 1;
}
/*----------------------------------------------------------------------------*/
static int FourMessScanPar(pSICSOBJ self, SConnection * pCon,
pHdb commandNode, pHdb par[], int nPar)
{
pFourMess priv = self->pPrivate;
char *scanvar;
double stt, dVal;
int np, preset;
pDynString data;
assert(priv != NULL);
if (nPar < 1) {
SCWrite(pCon, "ERROR: need two theta parameter to scanpar", eError);
return 0;
}
stt = par[0]->value.v.doubleValue;
scanvar = GetFourCircleScanVar(priv->stepTable, stt);
dVal = GetFourCircleStep(priv->stepTable, stt);
np = GetFourCircleScanNP(priv->stepTable, stt);
preset = GetFourCirclePreset(priv->stepTable, stt);
if (strcmp(scanvar, "Not found") == 0) {
SCPrintf(pCon, eValue, "%s,%f,%d,%d", "om", dVal, np, preset);
} else {
SCPrintf(pCon, eValue, "%s,%f,%d,%d", scanvar, dVal, np, preset);
}
return 1;
}
/*--------------------------------------------------------------------------*/
static int FourMessPrepareScan(pScanData self)
{
pVarEntry pVar = NULL;
void *pDings;
int i, iRet;
float fVal;
char pBueffel[512];
char pMessage[1024];
assert(self);
assert(self->pCon);
/* check boundaries of scan variables and allocate storage */
for (i = 0; i < self->iScanVar; i++) {
DynarGet(self->pScanVar, i, &pDings);
pVar = (pVarEntry) pDings;
if (pVar) {
iRet = CheckScanVar(pVar, self->pCon, self->iNP - 1);
if (!iRet) {
return 0;
}
InitScanVar(pVar);
} else {
SCWrite(self->pCon,
"WARNING: Internal error, no scan variable, I try to continue",
eLog);
}
pVar = NULL;
} /* end for */
/* configure counter */
SetCounterMode((pCounter) self->pCounterData, self->iMode);
SetCounterPreset((pCounter) self->pCounterData, self->fPreset);
self->iCounts = 0;
return 1;
}
/*-----------------------------------------------------------------------------*/
static int FourMessPrepareCmd(pSICSOBJ self, SConnection * pCon,
pHdb commandNode, pHdb par[], int nPar)
{
pFourMess priv = self->pPrivate;
return FourMessPrepareScan(priv->pScanner);
}
/*----------------------------------------------------------------------------*/
static hdbCallbackReturn SetScannerCB(pHdb node, void *userData,
pHdbMessage mm)
{
pHdbDataMessage set = NULL;
pFourMess priv = (pFourMess) userData;
SConnection *pCon = NULL;
pScanData old;
if ((set = GetHdbSetMessage(mm)) != NULL) {
old = priv->pScanner;
priv->pScanner =
FindCommandData(pServ->pSics, set->v->v.text, "ScanObject");
if (priv->pScanner == NULL) {
priv->pScanner = old;
pCon = set->callData;
if (pCon != NULL) {
SCWrite(pCon, "ERROR: scan object not found", eError);
}
return hdbAbort;
}
}
return hdbContinue;
}
/*---------------------------------------------------------------------------*/
static double getProtonAverage(pFourMess self)
{
int np, i;
long *lData = NULL, lSum = 0;
np = GetScanNP(self->pScanner);
lData = (long *) malloc((np + 1) * sizeof(long));
if (lData == NULL || np == 0) {
return 0.;
}
memset(lData, 0, (np + 1) * sizeof(long));
GetScanMonitor(self->pScanner, 2, lData, np);
for (i = 0; i < np; i++) {
lSum += lData[i];
}
free(lData);
return (double) lSum / (double) np;
}
/*---------------------------------------------------------------------------*/
static int FourMessStoreIntern(pSICSOBJ self, SConnection * pCon,
double fHkl[3], double fPosition[4], char *extra)
{
pFourMess priv = self->pPrivate;
float fSum, fSigma, fTemp, fStep = .0, fPreset =.0, fMF;
int i, iLF, iRet, iNP, ii;
long *lCounts = NULL;
pEVControl pEva = NULL;
pDummy pPtr = NULL;
pIDrivable pDriv = NULL;
char pBueffel[512], pTime[512], pNum[10];
double prot;
if (priv->pScanner == NULL) {
SCWrite(pCon, "ERROR: store: scan not configured", eLogError);
return 0;
}
if (priv->hklFile == NULL) {
SCWrite(pCon, "ERROR: store: no files open", eLogError);
return 0;
}
priv->count++;
/* get necessary data */
fSum = 0.;
fSigma = 0.;
iRet = ScanIntegrate(priv->pScanner, &fSum, &fSigma);
if (iRet != 1) {
switch (iRet) {
case INTEGLEFT:
sprintf(pBueffel,
"WARNING: integration failed --> no left side to: %f %f %f",
fHkl[0], fHkl[1], fHkl[2]);
break;
case INTEGRIGHT:
sprintf(pBueffel,
"WARNING: integration failed -->no right side to: %f %f %f",
fHkl[0], fHkl[1], fHkl[2]);
break;
case INTEGNOPEAK:
sprintf(pBueffel,
"WARNING: integration failed -->no peak found: %f %f %f",
fHkl[0], fHkl[1], fHkl[2]);
break;
case INTEGFUNNYBACK:
sprintf(pBueffel,
"WARNING: integration problem, asymmetric background: %f %f %f",
fHkl[0], fHkl[1], fHkl[2]);
break;
}
SCWrite(pCon, pBueffel, eLog);
}
iNP = GetScanNP(priv->pScanner);
lCounts = malloc(iNP * sizeof(long));
if (lCounts == NULL) {
SCWrite(pCon, "ERROR: out of memory in store", eLogError);
return 0;
}
GetScanCounts(priv->pScanner, lCounts, iNP);
/* write it */
if (priv->profFile) {
fprintf(priv->profFile,
"%4d %7.3f %7.3f %7.3f %7.2f %7.2f %7.2f %7.2f %7.0f %7.2f\n",
priv->count, fHkl[0], fHkl[1], fHkl[2], fPosition[0],
fPosition[1], fPosition[2], fPosition[3], fSum, fSigma);
}
if (priv->hklFile) {
fprintf(priv->hklFile,
"%5d %6.2f %6.2f %6.2f %7.2f %7.2f %7.2f %7.2f %7.0f %7.2f\n",
priv->count, fHkl[0], fHkl[1], fHkl[2], fPosition[0],
fPosition[1], fPosition[2], fPosition[3], fSum, fSigma);
}
sprintf(pBueffel,
"%5d %6.2f %6.2f %6.2f %7.2f %7.2f %7.2f %7.2f %7.0f %7.2f\n",
priv->count, fHkl[0], fHkl[1], fHkl[2], fPosition[0],
fPosition[1], fPosition[2], fPosition[3], fSum, fSigma);
SCWrite(pCon, pBueffel, eLog);
/* get temperature */
fTemp = -777.77;
pEva = (pEVControl) FindCommandData(pServ->pSics, "temperature",
"Environment Controller");
if (pEva == NULL) {
pPtr = (pDummy) FindCommandData(pServ->pSics, "temperature",
"RemObject");
if (pPtr != NULL) {
pDriv = pPtr->pDescriptor->GetInterface(pPtr, DRIVEID);
if (pDriv != NULL) {
fTemp = pDriv->GetValue(pPtr, pCon);
}
}
} else {
iRet = EVCGetPos(pEva, pCon, &fTemp);
}
/* get mf */
fMF = -0.00;
pEva = (pEVControl) FindCommandData(pServ->pSics, "mf",
"Environment Controller");
if (pEva == NULL) {
pPtr = (pDummy) FindCommandData(pServ->pSics, "mf",
"RemObject");
if (pPtr != NULL) {
pDriv = pPtr->pDescriptor->GetInterface(pPtr, DRIVEID);
if (pDriv != NULL) {
fMF = pDriv->GetValue(pPtr, pCon);
}
}
} else {
iRet = EVCGetPos(pEva, pCon, &fMF);
}
/* write profile */
if (priv->profFile) {
/* collect data */
SNXFormatTime(pBueffel, 512);
GetScanVarStep(priv->pScanner, 0, &fStep);
fPreset = GetScanPreset(priv->pScanner);
prot = getProtonAverage(priv);
/*
They rather wanted fMF in place of the proton average which fell victim to the
monitor assignement chaos anyway.
*/
if(extra == NULL){
fprintf(priv->profFile, "%3d %7.4f %9.0f %7.3f %12f %s\n", iNP, fStep,
fPreset, fTemp, fMF, pBueffel);
} else {
fprintf(priv->profFile, "%3d %7.4f %9.0f %7.3f %12f %s %s\n", iNP, fStep,
fPreset, fTemp, fMF, pBueffel,extra);
}
for (i = 0; i < iNP; i++) {
for (ii = 0; ii < 10 && i < iNP; ii++) {
fprintf(priv->profFile, " %7ld", lCounts[i]);
iLF = 1;
i++;
}
fprintf(priv->profFile, "\n");
i--;
iLF = 0;
}
if (iLF) {
fprintf(priv->profFile, "\n");
}
fflush(priv->profFile);
}
strcpy(pTime, pBueffel);
sprintf(pBueffel, "%3d%8.4f%10.0f%8.3f %s\n", iNP, fStep,
fPreset, fTemp, pTime);
SCWrite(pCon, pBueffel, eLog);
pBueffel[0] = '\0';
for (i = 0; i < iNP; i++) {
for (ii = 0; ii < 10 && i < iNP; ii++) {
sprintf(pNum, " %6ld", lCounts[i]);
strcat(pBueffel, pNum);
iLF = 1;
i++;
}
SCWrite(pCon, pBueffel, eLog);
pBueffel[0] = '\0';
i--;
iLF = 0;
}
if (iLF) {
SCWrite(pCon, pBueffel, eLog);
}
free(lCounts);
return 1;
}
/*----------------------------------------------------------------------------*/
static int FourMessStore(pSICSOBJ self, SConnection * pCon,
pHdb commandNode, pHdb par[], int nPar)
{
double fHkl[3], fPosition[4];
int i;
if (nPar < 7) {
SCWrite(pCon, "ERROR: not enough arguments for store", eError);
return 0;
}
/* load hkl */
for (i = 0; i < 3; i++) {
fHkl[i] = par[i]->value.v.doubleValue;
}
/* load positions */
for (i = 0; i < 4; i++) {
fPosition[i] = par[i + 3]->value.v.doubleValue;
}
return FourMessStoreIntern(self, pCon, fHkl, fPosition, NULL);
}
/*----------------------------------------------------------------------------*/
static int FourMessStoreExtra(pSICSOBJ self, SConnection * pCon,
pHdb commandNode, pHdb par[], int nPar)
{
double fHkl[3], fPosition[4];
int i;
if (nPar < 8) {
SCWrite(pCon, "ERROR: not enough arguments for storeextra", eError);
return 0;
}
/* load hkl */
for (i = 0; i < 3; i++) {
fHkl[i] = par[i]->value.v.doubleValue;
}
/* load positions */
for (i = 0; i < 4; i++) {
fPosition[i] = par[i + 3]->value.v.doubleValue;
}
return FourMessStoreIntern(self, pCon, fHkl, fPosition, par[7]->value.v.text);
}
/*------------------------------------------------------------------*/
static int weakScan(pSICSOBJ self, SConnection * pCon)
{
int i, np;
long low = 99999, high = -99999, *lCounts = NULL;
pFourMess priv = self->pPrivate;
hdbValue v;
/*
the scan is always OK if we do not test for weak conditions or we are in psd mode
*/
SICSHdbGetPar(self, pCon, "weak", &v);
if (v.v.intValue == 0) {
return 0;
}
np = GetScanNP(priv->pScanner);
lCounts = malloc(np * sizeof(long));
if (lCounts == NULL) {
SCWrite(pCon, "ERROR: out of memory in weakScan test", eLogError);
return 0;
}
GetScanCounts(priv->pScanner, lCounts, np);
for (i = 0; i < np; i++) {
if (lCounts[i] < low) {
low = lCounts[i];
}
if (lCounts[i] > high) {
high = lCounts[i];
}
}
/*
I am using the weakest point here as a rough estimate of
the background
*/
SICSHdbGetPar(self, pCon, "weakthreshold", &v);
if (high - 2 * low > v.v.intValue) {
return 0;
} else {
return 1;
}
}
/*----------------------------------------------------------------------------*/
static int FourMessWeak(pSICSOBJ self, SConnection * pCon,
pHdb commandNode, pHdb par[], int nPar)
{
int weak;
weak = weakScan(self, pCon);
SCPrintf(pCon, eLog, "weak = %d", weak);
return 1;
}
/*----------------------------------------------------------------------------
* Usage of the suppress flag:
* 0 do notthing
* 1 suppress symmetriqually equivalent reflections
* 2 opposite: add only reflections NOT allowed by the spacegroup
* ---------------------------------------------------------------------------*/
static int GenIndex(pSICSOBJ self, SConnection * pCon, pHdb commandNode,
pHdb par[], int nPar)
{
double lambda = SXGetLambda(), d, om, hkl[3];
int h, k, l, minh, mink, minl, suppress;
hdbValue hkllim, sttlim;
T_SgInfo *sginfo = SXGetSpaceGroup();
pFourMess priv = self->pPrivate;
int count = 0;
MATRIX B, H, Z1;
const double *cell;
lattice direct;
if (nPar < 1) {
SCWrite(pCon, "ERROR: need a suppression flag for indgen", eError);
return 0;
}
SICSHdbGetPar(self, pCon, "hkllim", &hkllim);
SICSHdbGetPar(self, pCon, "sttlim", &sttlim);
suppress = par[0]->value.v.intValue;
minh = hkllim.v.intArray[3];
mink = hkllim.v.intArray[4];
minl = hkllim.v.intArray[5];
SetListMin_hkl(sginfo, hkllim.v.intArray[1], hkllim.v.intArray[2],
&minh, &mink, &minl);
ClearReflectionList(priv->messList);
cell = SXGetCell();
direct.a = cell[0];
direct.b = cell[1];
direct.c = cell[2];
direct.alpha = cell[3];
direct.beta = cell[4];
direct.gamma = cell[5];
B = mat_creat(3, 3, UNIT_MATRIX);
if (!calculateBMatrix(direct, B)) {
SCWrite(pCon, "ERROR: invalid cell", eError);
return 0;
}
H = mat_creat(3, 1, ZERO_MATRIX);
for (h = hkllim.v.intArray[0]; h <= hkllim.v.intArray[3]; h++) {
for (k = hkllim.v.intArray[1]; k <= hkllim.v.intArray[4]; k++) {
for (l = hkllim.v.intArray[2]; l <= hkllim.v.intArray[5]; l++) {
/* SCPrintf(pCon,eLog, "Testing %d, %d, %d", h,k,l); */
/* first test: extinct */
if(suppress != 2){
if (IsSysAbsent_hkl(sginfo, h, k, l, NULL) != 0) {
/* SCPrintf(pCon,eLog, "%d, %d, %d rejected for sys absent", h, k, l); */
continue;
}
} else {
if (!IsSysAbsent_hkl(sginfo, h, k, l, NULL) != 0) {
/* SCPrintf(pCon,eLog, "%d, %d, %d rejected for not sys absent", h, k, l); */
continue;
}
}
/* second test: a symmetrically equivalent already seen */
if ((suppress > 0) && IsSuppressed_hkl(sginfo, minh, mink, minl,
hkllim.v.intArray[1],
hkllim.v.intArray[2], h, k,
l) != 0) {
/* SCPrintf(pCon,eLog, "%d, %d, %d rejected for symmetrical eqiv", h, k, l); */
continue;
}
/* third test: within stt limits */
H[0][0] = (double) h;
H[1][0] = (double) k;
H[2][0] = (double) l;
Z1 = mat_mul(B, H);
calcTheta(lambda, Z1, &d, &om);
om *= 2.;
mat_free(Z1);
if (om > sttlim.v.floatArray[0] && om < sttlim.v.floatArray[1]) {
hkl[0] = (double) h;
hkl[1] = (double) k;
hkl[2] = (double) l;
AddRefIdx(priv->messList, hkl);
count++;
}
if(SCGetInterrupt(pCon) != eContinue){
SCWrite(pCon,"ERROR: reflection generation aborted", eError);
return 0;
}
}
}
}
mat_free(B);
mat_free(H);
priv->masterCount = count;
SCPrintf(pCon, eValue, "%d reflections generated", count);
return 1;
}
/*-----------------------------------------------------------------------------*/
static int GenInconsumerate(pSICSOBJ self, SConnection * pCon,
pHdb commandNode, pHdb par[], int nPar)
{
double hkl[3], qvec[3];
pFourMess priv = self->pPrivate;
int i, j, iGen = 0, startCount;
if (nPar < 3) {
SCWrite(pCon,
"ERROR: need q displacement vector with three components",
eError);
return 0;
}
qvec[0] = par[0]->value.v.doubleValue;
qvec[1] = par[1]->value.v.doubleValue;
qvec[2] = par[2]->value.v.doubleValue;
startCount = priv->masterCount;
for (i = 0; i < startCount; i++) {
GetRefIndex(priv->messList, i, hkl);
for (j = 0; j < 3; j++) {
hkl[j] += qvec[j];
}
AddRefIdx(priv->messList, hkl);
iGen++;
GetRefIndex(priv->messList, i, hkl);
for (j = 0; j < 3; j++) {
hkl[j] -= qvec[j];
}
if(FindHKL(priv->messList, hkl[0], hkl[1], hkl[2]) == NULL){
AddRefIdx(priv->messList, hkl);
iGen++;
}
if(SCGetInterrupt(pCon) != eContinue){
SCWrite(pCon,"ERROR: generating incommensurate reflections aborted", eError);
return 0;
}
if( (i % 50) == 0 ){
SCPrintf(pCon,eLog, "%d of %d input reflections processed", i, startCount);
}
}
/*
add satellites of 0,0,0,
*/
for(i = 0; i < 3; i++){
hkl[i] = 0. + qvec[i];
}
AddRefIdx(priv->messList, hkl);
iGen++;
for(i = 0; i < 3; i++){
hkl[i] = 0. - qvec[i];
}
if(FindHKL(priv->messList, hkl[0], hkl[1], hkl[2]) == NULL){
AddRefIdx(priv->messList, hkl);
iGen++;
}
SCPrintf(pCon, eValue,
"%d additional inconsumerate reflections generated", iGen);
return 1;
}
/*----------------------------------------------------------------------------*/
static int hklCompare(const void *h1, const void *h2)
{
const double *hkl1 = h1, *hkl2 = h2;
if (hkl1[3] < hkl2[3]) {
return -1;
} else if (hkl1[3] == hkl2[3]) {
return 0;
} else {
return 1;
}
}
/*----------------------------------------------------------------------------
* This sorts on two theta only.
-----------------------------------------------------------------------------*/
static int SortRef(pSICSOBJ self, SConnection * pCon, pHdb commandNode,
pHdb par[], int nPar)
{
double *sortlist, d, lambda, om, hkl[4], ang[4];
const double *cell;
double *hklc;
int nRefl, i, j;
MATRIX B, H, Z1;
lattice direct;
pFourMess priv = self->pPrivate;
lambda = SXGetLambda();
cell = SXGetCell();
direct.a = cell[0];
direct.b = cell[1];
direct.c = cell[2];
direct.alpha = cell[3];
direct.beta = cell[4];
direct.gamma = cell[5];
B = mat_creat(3, 3, UNIT_MATRIX);
if (!calculateBMatrix(direct, B)) {
SCWrite(pCon, "ERROR: invalid cell", eError);
return 0;
}
H = mat_creat(3, 1, ZERO_MATRIX);
nRefl = ReflectionListCount(priv->messList);
sortlist = malloc(nRefl * 4 * sizeof(double));
if (sortlist == NULL) {
SCWrite(pCon, "ERROR: out of memory in SortRef", eError);
return 0;
}
/*
* I am using hkl[3] for storing the theta value to sort for! */
for (i = 0; i < nRefl; i++) {
GetRefIndex(priv->messList, i, hkl);
for (j = 0; j < 3; j++) {
H[j][0] = hkl[j];
}
Z1 = mat_mul(B, H);
calcTheta(lambda, Z1, &d, &om);
mat_free(Z1);
hkl[3] = om;
memcpy(sortlist + i * 4, hkl, 4 * sizeof(double));
}
qsort(sortlist, nRefl, 4 * sizeof(double), hklCompare);
ClearReflectionList(priv->messList);
for (i = 0; i < nRefl; i++) {
ang[1] = sortlist[i * 4 + 3];
ang[0] = 2. * ang[1];
ang[2] = .0;
ang[3] = .0;
/* AddRefIdxAng(priv->messList, sortlist+i*4,ang); */
AddRefIdx(priv->messList, sortlist + i * 4);
}
free(sortlist);
mat_free(B);
mat_free(H);
SCSendOK(pCon);
return 1;
}
/*----------------------------------------------------------------------------
* Sorting reflection is actually a variant of the traveling salesman
* problem. Optimal solutions for this are computationally expensive.
* Optimise for the shortest path between reflections. Not all motors
* are equal though. Om and phi move quick whereas stt and chi go slow.
* As a compromise, what I try here is sorting by a weighted addition of
* stt and chi.
----------------------------------------------------------------------------*/
static int SortRefNew(pSICSOBJ self, SConnection * pCon, pHdb commandNode,
pHdb par[], int nPar)
{
double *sortlist, d, lambda, om, hkl[4], ang[5];
const double *cell;
double *hklc;
int nRefl, i, j, status, count;
SingleXModes mode;
pFourMess priv = self->pPrivate;
pSingleDiff diffi = NULL;
diffi = SXGetDiffractometer();
mode = SXGetMode();
if(diffi == NULL){
SCWrite(pCon,"ERROR: not initialized, no diffractometer", eError);
return 0;
}
nRefl = ReflectionListCount(priv->messList);
sortlist = malloc(nRefl * 4 * sizeof(double));
if (sortlist == NULL) {
SCWrite(pCon, "ERROR: out of memory in SortRef", eError);
return 0;
}
/*
* I am using hkl[3] for storing the sort value to sort for! */
for (i = 0, count = 0; i < nRefl; i++) {
GetRefIndex(priv->messList, i, hkl);
status = diffi->calculateSettings(diffi,hkl,ang);
if(status == 1){
if(mode == Bisecting){
hkl[3] = ang[0] + 0.5* ang[2];
if(ang[2] < 0){
SCWrite(pCon,"WARNING: chi < 0!", eWarning);
}
} else if(mode == NB){
/*
* sort for nu in the first place. In order to cope with negativity, add
* 10.
*/
hkl[3] = ang[2]+10 + 0.1 * ang[0];
} else if(mode == BiNB){
/*
* sort for nu in the first place. In order to cope with negativity, add
* 10.
*/
hkl[3] = ang[4]+10 + 0.1 * ang[0];
} else {
hkl[3] = ang[0];
}
memcpy(sortlist + count * 4, hkl, 4 * sizeof(double));
count++;
}
}
qsort(sortlist, count, 4 * sizeof(double), hklCompare);
ClearReflectionList(priv->messList);
for (i = 0; i < count; i++) {
hklc = sortlist+i*4;
hklc[3] = .0; /* otherwise this will be interpreted as psi! */
status = diffi->calculateSettings(diffi,hklc,ang);
AddRefIdxAng(priv->messList, sortlist+i*4,ang);
/* AddRefIdx(priv->messList, sortlist + i * 4); */
}
free(sortlist);
SCSendOK(pCon);
return 1;
}
/*----------------------------------------------------------------------------
* Write a reflection list in the TSPLIB format for running a traveling
* salesman algorithm against it.
----------------------------------------------------------------------------- */
static int WriteTSP(pSICSOBJ self, SConnection * pCon, pHdb commandNode,
pHdb par[], int nPar)
{
double hkl[4], ang[4];
int nRefl, i, j, status, count;
SingleXModes mode;
pFourMess priv = self->pPrivate;
pSingleDiff diffi = NULL;
FILE *fd = NULL;
if(nPar < 1) {
SCWrite(pCon,"ERROR: need a filename parameter", eError);
return 0;
}
fd = fopen(trim(par[0]->value.v.text),"w");
if(fd == NULL){
SCPrintf(pCon,eError,"ERROR: failed to open %s for writing", par[0]->value.v.text);
return 0;
}
diffi = SXGetDiffractometer();
mode = SXGetMode();
if(diffi == NULL){
SCWrite(pCon,"ERROR: not initialized, no diffractometer", eError);
return 0;
}
nRefl = ReflectionListCount(priv->messList);
/*
* smear header...
*/
fprintf(fd,"NAME : %s\n", par[0]->value.v.text);
fprintf(fd,"TYPE : TSP\n");
fprintf(fd,"DIMENSION : %d\n", nRefl);
fprintf(fd,"EDGE_WEIGHT_TYPE : EUC_3D\n");
fprintf(fd,"NODE_COORD_SECTION\n");
for (i = 0, count = 0; i < nRefl; i++) {
GetRefIndex(priv->messList, i, hkl);
status = diffi->calculateSettings(diffi,hkl,ang);
if(status == 1){
if(mode == Bisecting){
fprintf(fd,"%d %e %e %e\n",i+1,ang[0]*3*100,ang[2]*2*100,ang[3]*100);
} else if(mode == NB){
fprintf(fd,"%d %e %e %e\n",i,ang[2]*4*100,ang[0]*3*100,ang[1]*100);
} else {
fprintf(fd,"%d %e %e %e\n",i,ang[0]*3*100,ang[2]*2*100,ang[3]*2*100);
}
}
}
fclose(fd);
SCSendOK(pCon);
return 1;
}
/*----------------------------------------------------------------------------
* read a tour file as generated by LKH or other TSP solvers and
* reorder the reflection list accordingly.
-----------------------------------------------------------------------------*/
static int ReadTour(pSICSOBJ self, SConnection * pCon, pHdb commandNode,
pHdb par[], int nPar)
{
FILE *fd = NULL;
char buffer[132];
double hkl[4], ang[5];
int idx, tour;
pFourMess priv = self->pPrivate;
pSingleDiff diffi = NULL;
if(nPar < 1) {
SCWrite(pCon,"ERROR: need name of a tour file to read", eError);
return 0;
}
fd = fopen(par[0]->value.v.text,"r");
if(fd == NULL){
SCPrintf(pCon,eError,"ERROR: failed to open tour file %s for reading",
par[0]->value.v.text);
return 0;
}
tour = LLDcreate(sizeof(hkl));
diffi = SXGetDiffractometer();
assert(diffi != NULL);
/* skip header */
while(fgets(buffer,sizeof(buffer),fd) != NULL){
if(strstr(buffer,"TOUR_SECTION") != NULL){
break;
}
}
while(fscanf(fd,"%d",&idx) == 1){
if(idx < 0){
break;
}
GetRefIndex(priv->messList,idx-1,hkl);
LLDnodeAppendFrom(tour,hkl);
}
fclose(fd);
idx = LLDnodePtr2First(tour);
ClearReflectionList(priv->messList);
while(idx == 1){
LLDnodeDataTo(tour,hkl);
hkl[3] = 0;
diffi->calculateSettings(diffi,hkl,ang);
AddRefIdxAng(priv->messList, hkl,ang);
idx = LLDnodePtr2Next(tour);
}
LLDdelete(tour);
SCSendOK(pCon);
return 1;
}
/*----------------------------------------------------------------------------*/
static int FourMessSave(void *data, char *name, FILE * fd)
{
pSICSOBJ self = data;
pFourMess priv = self->pPrivate;
SaveSICSOBJ(data, name, fd);
priv->stepTable->pDes->SaveStatus(priv->stepTable,"fmess table", fd);
return 1;
}
/*----------------------------------------------------------------------------*/
void InstallFourMess(SConnection * pCon, SicsInterp * pSics)
{
pFourMess priv = NULL;
pSICSOBJ pNew = NULL;
pHdb cmd = NULL;
int hkl[] = { -10, -10, 10, 10, 10, 10 };
double sttlim[] = { 5.0, 180 };
pNew = MakeSICSOBJ("fmess", "FourMess");
priv = calloc(1, sizeof(FourMess));
if (pNew == NULL || priv == NULL) {
SCWrite(pCon, "ERROR: out of memory creating fourmess", eError);
}
pNew->pDes->SaveStatus = FourMessSave;
pNew->pPrivate = priv;
pNew->KillPrivate = KillFourMess;
priv->stepTable = MakeFourCircleTable();
AddCommand(pSics, "fmesstable", InterInvokeSICSOBJ, NULL, priv->stepTable);
priv->messList = CreateReflectionList(pCon, pSics, "messref");
if (priv->stepTable < 0 || priv->messList == NULL) {
SCWrite(pCon, "ERROR: out of memory creating fourmess", eError);
}
cmd = AddSICSHdbPar(pNew->objectNode, "weak", usUser, MakeHdbInt(0));
SetHdbProperty(cmd, "__save", "true");
cmd = AddSICSHdbPar(pNew->objectNode, "weakthreshold", usUser,
MakeHdbInt(20));
SetHdbProperty(cmd, "__save", "true");
cmd = AddSICSHdbPar(pNew->objectNode, "fast", usUser, MakeHdbInt(0));
SetHdbProperty(cmd, "__save", "true");
cmd =
AddSICSHdbPar(pNew->objectNode, "mode", usUser,
MakeHdbText("Monitor"));
cmd =
AddSICSHdbPar(pNew->objectNode, "hkllim", usUser,
MakeHdbIntArray(6, hkl));
SetHdbProperty(cmd, "__save", "true");
cmd =
AddSICSHdbPar(pNew->objectNode, "sttlim", usUser,
MakeHdbFloatArray(2, sttlim));
SetHdbProperty(cmd, "__save", "true");
cmd =
AddSICSHdbPar(pNew->objectNode, "start", usUser,
MakeSICSFunc(FourMessStart));
cmd = AddSICSHdbPar(cmd, "filename", usUser, MakeHdbText("Unknown"));
cmd =
AddSICSHdbPar(pNew->objectNode, "close", usUser,
MakeSICSFunc(FourMessClose));
cmd =
AddSICSHdbPar(pNew->objectNode, "reopen", usUser,
MakeSICSFunc(FourMessReopen));
cmd =
AddSICSHdbPar(pNew->objectNode, "prepare", usUser,
MakeSICSFunc(FourMessPrepareCmd));
cmd =
AddSICSHdbPar(pNew->objectNode, "scanpar", usUser,
MakeSICSFunc(FourMessScanPar));
cmd = AddSICSHdbPar(cmd, "twotheta", usUser, MakeHdbFloat(.0));
cmd =
AddSICSHdbPar(pNew->objectNode, "store", usUser,
MakeSICSFunc(FourMessStore));
AddSICSHdbPar(cmd, "h", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "k", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "l", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "stt", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "om", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "chi", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "phi", usUser, MakeHdbFloat(.0));
cmd =
AddSICSHdbPar(pNew->objectNode, "storeextra", usUser,
MakeSICSFunc(FourMessStoreExtra));
AddSICSHdbPar(cmd, "h", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "k", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "l", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "stt", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "om", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "chi", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "phi", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "extra", usUser, MakeHdbText(""));
cmd =
AddSICSHdbPar(pNew->objectNode, "weak", usUser,
MakeSICSFunc(FourMessWeak));
cmd =
AddSICSHdbPar(pNew->objectNode, "indgen", usUser,
MakeSICSFunc(GenIndex));
AddSICSHdbPar(cmd, "sup", usUser, MakeHdbInt(1));
cmd =
AddSICSHdbPar(pNew->objectNode, "genw", usUser,
MakeSICSFunc(GenInconsumerate));
AddSICSHdbPar(cmd, "hw", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "kw", usUser, MakeHdbFloat(.0));
AddSICSHdbPar(cmd, "lw", usUser, MakeHdbFloat(.0));
cmd =
AddSICSHdbPar(pNew->objectNode, "indsort", usUser,
MakeSICSFunc(SortRef));
cmd =
AddSICSHdbPar(pNew->objectNode, "writetsp", usUser,
MakeSICSFunc(WriteTSP));
AddSICSHdbPar(cmd, "filename", usUser, MakeHdbText("Unknown"));
cmd =
AddSICSHdbPar(pNew->objectNode, "readtour", usUser,
MakeSICSFunc(ReadTour));
AddSICSHdbPar(cmd, "filename", usUser, MakeHdbText("Unknown"));
cmd =
AddSICSHdbPar(pNew->objectNode, "template", usMugger,
MakeHdbText("Unknown"));
cmd =
AddSICSHdbPar(pNew->objectNode, "scanobj", usMugger,
MakeHdbText("xxxscan"));
PrependHipadabaCallback(cmd,
MakeHipadabaCallback(SetScannerCB, priv, NULL));
priv->pScanner = FindCommandData(pSics, "xxxscan", "ScanObject");
/* AddHipadabaChild(pNew->objectNode, priv->stepTable->objectNode,pCon); */
AddCommand(pSics, "fmess", FourMessAction, KillSICSOBJ, pNew);
}
Reference in New Issue View Git Blame Copy Permalink