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