680 lines
18 KiB
C
680 lines
18 KiB
C
/**
|
|
* 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;
|
|
}
|
|
|