- Reworked the connection object and the IO system
- Reworked the support for TRICS - Added a second generation motor
This commit is contained in:
koennecke
679
simidx.c
Normal file
679
simidx.c
Normal file
@@ -0,0 +1,679 @@
|
||||
/**
|
||||
* 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 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;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user