/** * 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 20 #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(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); return 1; } /*-------------------------------------------------------------*/ 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) > .00001){ return 0; } else { return 1; } } /*-------------------------------------------------------------- * - 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; triplet[0] = 0; /* * test for 180 */ while(idx < nReflections){ angle = angleBetweenScatVec(reflections[0].UVW, reflections[idx].UVW); if(angle < 160 && angle > -160){ triplet[1] = idx; break; } idx++; } if(idx >= nReflections){ SimIdxPrint(1,"ERROR: no second index found"); return 0; } for(idx = 1; idx < nReflections; idx++){ if(idx != triplet[1]) { if(!areCoplanar(reflections[triplet[0]].UVW, reflections[triplet[1]].UVW, reflections[idx].UVW)){ triplet[2] = idx; 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 -------------------------------------------------------------*/ 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 1; } else { return 0; } } /*-------------------------------------------------------------*/ 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; solutions[nSolutions] = is; nSolutions++; } /*----------------------------------------------------------------------*/ 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; } /*------------------------------------------------------------- * If the system is right handed the determinat of the * matrix having the indices as columns must be positive. * As I have only two vectors, I simulate the third by * using the nromal on the other two. -------------------------------------------------------------*/ static int testDuoRightHandedness(int r1, int r1idx, int r2, int r2idx){ MATRIX T; double vol; int r3, r3idx; 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 1; } else { return 0; } } /*--------------------------------------------------------------*/ 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()){ 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,1); 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 retrieval ===================*/ int SimIdxGetNSolutions(){ return 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; }