linse/sics
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
main
sics/simidx.c
koennecke 9eca96b064 - Fixes to make SL6 work 
- New NeXus libraries
- Added new raw binary transfer mode for mass data
- Added a check script option to configurable virtual motor


SKIPPED:
	psi/dumprot.c
	psi/make_gen
	psi/psi.c
	psi/rebin.c
	psi/sanslirebin.c
2012-03-29 08:41:05 +00:00

763 lines
20 KiB
C
Raw Permalink Blame History

/**
* This is the simple reflection indexer. The algorithm is simple:
* - Three non coplanar reflections at low two theta are selected.
* - Candidate indices are calculated from the lattice constants
* - Permutations of the generated indices are changed for a couple
* of conditions:
* -- Does the angle between the reflections matches the expectaions
* -- Do the reflections form a right handed set
* UB matrics matching these conditions are calculated and stored
* for later retrieval.
*
* The software is organized such that this is a standalone module.
* For reasons of laziness I use statically sized arrays for
* candidate indices and for the reflection list. This is simple and
* this code is for small problems anyway.
*
* Mark Koennecke, August 2008
*/
#include <math.h>
#include <assert.h>
#include <string.h>
#include <stdlib.h>
#include <stdarg.h>
#include "simidx.h"
#include "cell.h"
#include "vector.h"
#include "ubfour.h"
#include "fourlib.h"
/*======================== defines ===========================*/
#define MAXCANDIDATES 20
#define MAXREF 20
#define MAXIDX 20
#define MAXSOLUTION 50
#define ABS(x) (x < 0 ? -(x) : (x))
/*======================== types =============================*/
typedef struct {
int h, k, l;
double diff;
} HKL, *pHKL;
typedef struct {
double uvw[3];
MATRIX UVW;
HKL indices[MAXCANDIDATES];
int nHKL;
int currentIDX;
double twotheta;
int originalID;
} IndexVector, *pIndexVector;
/*================== module parameters ========================*/
static lattice direct;
static double lambda;
static T_SgInfo *spgrp = NULL;
static IndexVector reflections[MAXREF];
static int nReflections = 0;
static double sttlim = .5, anglim = 2.;
static OutFunc outFunc = NULL;
static void *userData;
static int outLevel = 10;
static IndexSolution solutions[MAXSOLUTION];
static int nSolutions;
/*------------------------------------------------------------*/
void SimIdxInit()
{
int i;
for (i = 0; i < MAXREF; i++) {
reflections[i].UVW = mat_creat(3, 1, ZERO_MATRIX);
}
}
/*=============== configuration functions =====================*/
void SimIdxSetCell(double cell[6])
{
direct.a = cell[0];
direct.b = cell[1];
direct.c = cell[2];
direct.alpha = cell[3];
direct.beta = cell[4];
direct.gamma = cell[5];
}
/*-------------------------------------------------------------*/
void SimIdxSetLambda(double lmda)
{
lambda = lmda;
}
/*-------------------------------------------------------------*/
void SimIdxSetSttLim(double lmda)
{
sttlim = lmda;
}
/*-------------------------------------------------------------*/
void SimIdxSetAngLim(double lmda)
{
anglim = lmda;
}
/*-------------------------------------------------------------*/
void SimIdxSetSpacegroup(T_SgInfo * sg)
{
spgrp = sg;
}
/*-------------------------------------------------------------*/
void SimIdxClearReflection()
{
nReflections = 0;
}
/*-------------------------------------------------------------*/
void SimIdxAddReflection(double uvw[3])
{
int i;
if (nReflections < MAXREF) {
memcpy(&reflections[nReflections].uvw, uvw, 3 * sizeof(double));
reflections[nReflections].nHKL = 0;
reflections[nReflections].currentIDX = 0;
reflections[nReflections].originalID = nReflections;
for (i = 0; i < 3; i++) {
reflections[nReflections].UVW[i][0] = uvw[i];
}
nReflections++;
}
}
/*-------------------------------------------------------------*/
void SimIdxOutput(void *data, OutFunc out, int level)
{
userData = data;
outFunc = out;
outLevel = level;
}
/*-------------------------------------------------------------*/
static void SimIdxPrint(int level, char *fmt, ...)
{
va_list ap;
char buf[1024];
int l;
if (level > outLevel) {
return;
}
va_start(ap, fmt);
l = vsnprintf(buf, sizeof buf, fmt, ap);
va_end(ap);
if (outFunc != NULL) {
outFunc(userData, buf);
} else {
printf("%s\n", buf);
}
}
/*=============== The alkoholism ===============================*/
static int thetaCompare(const void *d1, const void *d2)
{
pIndexVector iv1, iv2;
iv1 = (pIndexVector) d1;
iv2 = (pIndexVector) d2;
if (iv1->twotheta == iv2->twotheta) {
return 0;
}
if (iv1->twotheta < iv2->twotheta) {
return -1;
} else {
return 1;
}
}
/*--------------------------------------------------------------*/
static void calcRefTheta()
{
int i;
double theta, d;
for (i = 0; i < nReflections; i++) {
calcTheta(lambda, reflections[i].UVW, &d, &theta);
reflections[i].twotheta = 2. * theta;
}
qsort(reflections, nReflections, sizeof(IndexVector), thetaCompare);
SimIdxPrint(10, "%d Reflections", nReflections);
}
/*---------------------------------------------------------------*/
double calcIdxTwoTheta(int h, int k, int l, MATRIX B)
{
MATRIX H, Z1;
double om, d;
H = mat_creat(3, 1, ZERO_MATRIX);
if (H == NULL) {
SimIdxPrint(1, "ERROR: out of memory calculating H matrix");
return 0.;
}
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);
mat_free(H);
return om;
}
/*---------------------------------------------------------------*/
static void AddCandidate(int n, int h, int k, int l, double diff)
{
int cur = reflections[n].nHKL;
if (cur < MAXCANDIDATES) {
reflections[n].indices[cur].h = h;
reflections[n].indices[cur].k = k;
reflections[n].indices[cur].l = l;
reflections[n].indices[cur].diff = diff;
reflections[n].nHKL++;
}
}
/*---------------------------------------------------------------*/
static int calcIndexes()
{
int h, k, l, i, status;
int minh, mink, minl;
MATRIX B;
double twotheta;
B = mat_creat(3, 3, UNIT_MATRIX);
if (B == NULL) {
SimIdxPrint(1, "ERROR: out of memory calculating B matrix");
return 0;
}
status = calculateBMatrix(direct, B);
if (status < 0) {
SimIdxPrint(1,
"ERROR: invalid cell constants, failed to calculate B matrix");
return 0;
}
minh = -MAXIDX;
mink = -MAXIDX;
minl = -MAXIDX;
SetListMin_hkl(spgrp, MAXIDX, MAXIDX, &minh, &mink, &minl);
for(i = 0; i < nReflections; i++){
reflections[i].nHKL = 0;
}
for (h = MAXIDX; h > -MAXIDX; h--) {
for (k = MAXIDX; k > -MAXIDX; k--) {
for (l = MAXIDX; l > -MAXIDX; l--) {
if (IsSysAbsent_hkl(spgrp, h, k, l, NULL) != 0) {
continue;
}
twotheta = calcIdxTwoTheta(h, k, l, B);
for (i = 0; i < nReflections; i++) {
if (reflections[i].twotheta > twotheta - sttlim &&
reflections[i].twotheta < twotheta + sttlim) {
AddCandidate(i, h, k, l,
ABS(twotheta - reflections[i].twotheta));
}
}
}
}
}
mat_free(B);
status = 1;
for(i = 0; i < nReflections; i++){
if(reflections[i].nHKL < 1) {
SimIdxPrint(10,"Failed to find candidate indices for %4d", reflections[i].originalID);
status = 0;
}
}
return status;
}
/*-------------------------------------------------------------*/
static double angleBetweenScatVec(MATRIX v1, MATRIX v2)
{
double angle;
angle = angleBetween(v1, v2);
return angle;
}
/*---------------------------------------------------------------*/
static void printRefDiagnostic()
{
int i, j;
double angle;
SimIdxPrint(10, "Reflection List and Candidate Indices");
SimIdxPrint(10, " N STT U V W");
for (i = 0; i < nReflections; i++) {
SimIdxPrint(10, "%3.3d %8.4f %8.4f %8.4f %8.4f", i,
reflections[i].twotheta,
reflections[i].uvw[0],
reflections[i].uvw[1], reflections[i].uvw[2]);
for (j = 0; j < reflections[i].nHKL; j++) {
SimIdxPrint(10, "\t%4d %4d %4d",
reflections[i].indices[j].h,
reflections[i].indices[j].k,
reflections[i].indices[j].l);
}
}
SimIdxPrint(10, "Angles between reflections");
SimIdxPrint(10, "IDX1 IDX2 Angle");
for (i = 0; i < nReflections; i++) {
for (j = i; j < nReflections; j++) {
if (i != j) {
angle =
angleBetweenScatVec(reflections[i].UVW, reflections[j].UVW);
SimIdxPrint(10, "%3d %3d %8.2f", i, j, angle);
}
}
}
}
/*-------------------------------------------------------------*/
static double calculateVolume(double v1[3], double v2[3], double v3[3])
{
MATRIX m;
int i;
double vol;
m = mat_creat(3, 3, ZERO_MATRIX);
for (i = 0; i < 3; i++) {
m[i][0] = v1[i];
m[i][1] = v2[i];
m[i][2] = v3[i];
}
vol = mat_det(m);
mat_free(m);
return vol;
}
/*-------------------------------------------------------------*/
static int areCoplanar(MATRIX v1, MATRIX v2, MATRIX v3)
{
MATRIX norm;
double dot;
norm = vectorCrossProduct(v1, v2);
if (norm != NULL) {
dot = vectorDotProduct(norm, v3);
mat_free(norm);
} else {
dot = .0;
}
if (ABS(dot) > .01) {
return 0;
} else {
return 1;
}
}
/*-------------------------------------------------------------*/
static double calcCoplanar(MATRIX v1, MATRIX v2, MATRIX v3)
{
MATRIX norm;
double dot;
norm = vectorCrossProduct(v1, v2);
if (norm != NULL) {
dot = vectorDotProduct(norm, v3);
mat_free(norm);
} else {
dot = .0;
}
return ABS(dot);
}
/*--------------------------------------------------------------
* - We want the shortest vectors
* - We do not want vectors at 180 to each other
* - We do not want the three vectors to be coplanar
*-------------------------------------------------------------*/
static int chooseTriplet(int triplet[3])
{
double angle, vol;
int idx = 1, coIdx = 0;
double coplanarVector[MAXREF], coMax = -9999999.99;
memset(coplanarVector,0,MAXREF*sizeof(double));
triplet[0] = 0;
/*
* test for 180
*/
while (idx < nReflections) {
angle = angleBetweenScatVec(reflections[0].UVW, reflections[idx].UVW);
if (angle < 140 && angle > -140) {
triplet[1] = idx;
break;
}
idx++;
}
if (idx >= nReflections) {
SimIdxPrint(1, "ERROR: no second index found");
return 0;
}
/*
* Try to find a reflection which is most out of the plane build by the first
* two. A good angular separation will give a better UB
*/
for (idx = 1; idx < nReflections; idx++) {
if (idx != triplet[1]) {
coplanarVector[idx] = calcCoplanar(reflections[triplet[0]].UVW,
reflections[triplet[1]].UVW,
reflections[idx].UVW);
}
}
for(idx = 0; idx < nReflections; idx++){
if(coplanarVector[idx] > coMax){
coMax = coplanarVector[idx];
coIdx = idx;
}
}
if(coMax > .01){
triplet[2] = coIdx;
return 1;
}
SimIdxPrint(1, "ERROR: no three non coplanar reflections found");
return 0;
}
/*------------------------------------------------------------*/
static double reflectionsAngle(MATRIX B, int hkl1[3], int hkl2[3])
{
double angle;
reflection r1, r2;
r1.h = hkl1[0];
r1.k = hkl1[1];
r1.l = hkl1[2];
r2.h = hkl2[0];
r2.k = hkl2[1];
r2.l = hkl2[2];
return angleBetweenReflections(B, r1, r2);
}
/*-------------------------------------------------------------*/
static int findAngleMatch(MATRIX B, int idxr1, int r1,
int r2start, int r2, double *diff)
{
double scatAngle, hklAngle;
MATRIX H1, H2;
int i, r, hkl1[3], hkl2[3];
scatAngle = angleBetweenScatVec(reflections[r1].UVW,
reflections[r2].UVW);
hkl1[0] = reflections[r1].indices[idxr1].h;
hkl1[1] = reflections[r1].indices[idxr1].k;
hkl1[2] = reflections[r1].indices[idxr1].l;
for (i = r2start; i < reflections[r2].nHKL; i++) {
hkl2[0] = reflections[r2].indices[i].h;
hkl2[1] = reflections[r2].indices[i].k;
hkl2[2] = reflections[r2].indices[i].l;
hklAngle = reflectionsAngle(B, hkl1, hkl2);
*diff = ABS(scatAngle - hklAngle);
if (*diff < anglim) {
return i;
}
}
return -1;
}
/*-------------------------------------------------------------
* If the system is right handed the determinat of the
* matrix having the indices as columns must be positive
* I now (05/2009) thgink that this test is wrong. The testing
* has to be done separatly, by checking if the determinant of the
* calculated UB is > 0. As I do not want this in this file, which
* ought to be diffractometer independent, this test may need to
* left to a the caller.
-------------------------------------------------------------*/
static int testRightHandedness(int r1, int r1idx,
int r2, int r2idx, int r3, int r3idx)
{
MATRIX T;
double vol;
T = mat_creat(3, 3, ZERO_MATRIX);
if (T == NULL) {
return 0;
}
T[0][0] = reflections[r1].indices[r1idx].h;
T[1][0] = reflections[r1].indices[r1idx].k;
T[2][0] = reflections[r1].indices[r1idx].l;
T[0][1] = reflections[r2].indices[r2idx].h;
T[1][1] = reflections[r2].indices[r2idx].k;
T[2][1] = reflections[r2].indices[r2idx].l;
T[0][2] = reflections[r3].indices[r3idx].h;
T[1][2] = reflections[r3].indices[r3idx].k;
T[2][2] = reflections[r3].indices[r3idx].l;
vol = mat_det(T);
mat_free(T);
if (vol > .0) {
return 0;
}
return 1;
}
/*-------------------------------------------------------------*/
static void storeSolution(int r1, int r1idx,
int r2, int r2idx,
int r3, int r3idx, double diff)
{
IndexSolution is;
is.h[0] = reflections[r1].indices[r1idx].h;
is.k[0] = reflections[r1].indices[r1idx].k;
is.l[0] = reflections[r1].indices[r1idx].l;
is.originalID[0] = reflections[r1].originalID;
is.diff = reflections[r1].indices[r1idx].diff;
is.h[1] = reflections[r2].indices[r2idx].h;
is.k[1] = reflections[r2].indices[r2idx].k;
is.l[1] = reflections[r2].indices[r2idx].l;
is.originalID[1] = reflections[r2].originalID;
is.diff += reflections[r2].indices[r2idx].diff;
if (r3 != 999) {
is.h[2] = reflections[r3].indices[r3idx].h;
is.k[2] = reflections[r3].indices[r3idx].k;
is.l[2] = reflections[r3].indices[r3idx].l;
is.diff += reflections[r3].indices[r3idx].diff;
is.originalID[2] = reflections[r3].originalID;
} else {
is.h[2] = 0;
is.k[2] = 0;
is.l[2] = 0;
is.originalID[2] = 999;
}
is.diff += diff;
if(nSolutions < MAXSOLUTION){
solutions[nSolutions] = is;
nSolutions++;
} else {
SimIdxPrint(10,"WARNING: more solutions then solution space");
}
}
/*----------------------------------------------------------------------*/
static int compareSolution(const void *d1, const void *d2)
{
IndexSolution *iv1, *iv2;
iv1 = (IndexSolution *) d1;
iv2 = (IndexSolution *) d2;
if (iv1->diff == iv2->diff) {
return 0;
}
if (iv1->diff < iv2->diff) {
return -1;
} else {
return 1;
}
}
/*--------------------------------------------------------------*/
static int findSolutionsForTriplet(int triplet[3], int testRight)
{
int r1, r2, r3, i, status;
int match1, match2, r2start, r3start;
double diff1, diff2;
MATRIX B;
r1 = triplet[0];
r2 = triplet[1];
r3 = triplet[2];
B = mat_creat(3, 3, UNIT_MATRIX);
if (B == NULL) {
SimIdxPrint(1, "ERROR: out of memory calculating B matrix");
return 0;
}
status = calculateBMatrix(direct, B);
if (status < 0) {
SimIdxPrint(1,
"ERROR: invalid cell constants, failed to calculate B matrix");
return 0;
}
for (i = 0; i < reflections[r1].nHKL; i++) {
r2start = 0;
while ((match1 = findAngleMatch(B, i, r1, r2start, r2, &diff1)) >= 0) {
r3start = 0;
while ((match2 = findAngleMatch(B, i, r1, r3start, r3, &diff2)) >= 0) {
if (testRight == 1) {
if (testRightHandedness(r1, i, r2, match1, r3, match2)) {
storeSolution(r1, i, r2, match1, r3, match2, diff1 + diff2);
}
} else {
storeSolution(r1, i, r2, match1, r3, match2, diff1 + diff2);
}
r3start = match2 + 1;
}
r2start = match1 + 1;
}
}
qsort(solutions, nSolutions, sizeof(IndexSolution), compareSolution);
return 1;
}
/*--------------------------------------------------------------*/
static int findSolutionsForDuett(int triplet[3])
{
MATRIX B;
int r1, r2, r2start, i, status, match1;
double diff;
if (triplet[1] == 999) {
SimIdxPrint(1, "ERROR: No suitable reflection set found");
return 0;
}
r1 = triplet[0];
r2 = triplet[1];
B = mat_creat(3, 3, UNIT_MATRIX);
if (B == NULL) {
SimIdxPrint(1, "ERROR: out of memory calculating B matrix");
return 0;
}
status = calculateBMatrix(direct, B);
if (status < 0) {
SimIdxPrint(1,
"ERROR: invalid cell constants, failed to calculate B matrix");
return 0;
}
for (i = 0; i < reflections[r1].nHKL; i++) {
r2start = 0;
while ((match1 = findAngleMatch(B, i, r1, r2start, r2, &diff)) >= 0) {
storeSolution(r1, i, r2, match1, 999, 0, diff);
r2start = match1 + 1;
}
}
qsort(solutions, nSolutions, sizeof(IndexSolution), compareSolution);
return 1;
}
/*--------------------------------------------------------------
* This is used if we cannot find a solution for a triplet.
* Then we try to reduce to a duett. So we look for a second
* reflection in the original triplett which is closest
* to 90 degree in angle.
*/
static int secondForDuett(int triplet[3])
{
double diff1, diff2;
diff1 =
ABS(90. -
angleBetweenScatVec(reflections[triplet[0]].UVW,
reflections[triplet[1]].UVW));
diff2 =
ABS(90. -
angleBetweenScatVec(reflections[triplet[0]].UVW,
reflections[triplet[2]].UVW));
if (diff1 < diff2) {
return triplet[1];
} else {
return triplet[2];
}
}
/*--------------------------------------------------------------*/
int SimIdxRun()
{
int triplet[3] = { 999, 999, 999 }, status;
SimIdxPrint(10, "SimIdx calculating with parameters:");
SimIdxPrint(10, "Cell = %f %f %f %f %f %f", direct.a, direct.b, direct.c,
direct.alpha, direct.beta, direct.gamma);
SimIdxPrint(10, "Lambda = %f", lambda);
SimIdxPrint(10, "Sttlim, anglim = %f %f", sttlim, anglim);
nSolutions = 0;
calcRefTheta();
if (!calcIndexes()) {
SimIdxPrint(10,"ERROR: Problems finding candidate indices for reflections\nCannot continue\nCheck limits, cell and spacegroup");
return 0;
}
if (outLevel >= 10) {
printRefDiagnostic();
}
if (nReflections >= 3) {
if (!chooseTriplet(triplet)) {
return findSolutionsForDuett(triplet);
}
} else {
triplet[0] = 0;
triplet[1] = 1;
return findSolutionsForDuett(triplet);
}
SimIdxPrint(10, "Choosen triplet: %d, %d, %d\n", triplet[0],
triplet[1], triplet[2]);
status = findSolutionsForTriplet(triplet, 0);
/*
if (nSolutions == 0) {
SimIdxPrint(1,
"WARNING: found no right handed solution set, trying to find lefthanded");
status = findSolutionsForTriplet(triplet, 0);
}
*/
if (nSolutions == 0) {
SimIdxPrint(10, "Failed to find solution for triplet, trying duett");
status = secondForDuett(triplet);
triplet[1] = status;
status = findSolutionsForDuett(triplet);
}
return status;
}
/*=========================== solution management ===================*/
int SimIdxGetNSolutions()
{
return nSolutions;
}
/*------------------------------------------------------------------*/
void SimIdxRemoveSolution(int id){
if(id >= 0 && id < nSolutions){
solutions[id] = solutions[nSolutions-1];
nSolutions--;
}
}
/*------------------------------------------------------------------*/
IndexSolution SimIdxGetSolution(int id)
{
IndexSolution broken;
if (id >= 0 && id < nSolutions) {
return solutions[id];
}
broken.h[0] = -999;
broken.h[1] = -999;
broken.h[2] = -999;
return broken;
}
Reference in New Issue View Git Blame Copy Permalink