- First working version of the triple axis UB matrix code
This commit is contained in:
koennecke
335
tasublib.c
335
tasublib.c
@ -1,13 +1,15 @@
|
||||
/**
|
||||
* This is a library of functions and data structures for performing
|
||||
* triple axis spectrometer angle calculations using the UB-matrix
|
||||
* formalism as described by Mark Lumsden.
|
||||
* formalism as described by Mark Lumsden, to appear in Acta Cryst.
|
||||
*
|
||||
* copyright: see file COPYRIGHT
|
||||
*
|
||||
* Mark Koennecke, April 2005
|
||||
*/
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <assert.h>
|
||||
#include "trigd.h"
|
||||
#include "vector.h"
|
||||
#include "tasublib.h"
|
||||
@ -30,42 +32,53 @@ double KtoEnergy(double k){
|
||||
energy = ECONST*k*k;
|
||||
return energy;
|
||||
}
|
||||
/*-------------------------------------------------------------------*/
|
||||
static double calcCurvature(double B1, double B2, double theta){
|
||||
return B1 + B2/Sind(ABS(theta));
|
||||
}
|
||||
/*--------------------------------------------------------------------*/
|
||||
int maCalcAngles(maCrystal data, pmaAngles angles, double k){
|
||||
double fd;
|
||||
int maCalcTwoTheta(maCrystal data, double k, double *two_theta){
|
||||
double fd, theta;
|
||||
|
||||
/* fd = k/(2.*data.dd); */
|
||||
fd = PI/(data.dd*k);
|
||||
if(fd > 1.0) {
|
||||
return ENERGYTOBIG;
|
||||
}
|
||||
angles->theta = Asind(fd)*data.ss;
|
||||
angles->two_theta = 2.*angles->theta;
|
||||
angles->horizontal_curvature = data.HB1 + data.HB2/Sind(angles->theta);
|
||||
angles->vertical_curvature = data.VB1 + data.VB2/Sind(angles->theta);
|
||||
theta = Asind(fd)*data.ss;
|
||||
*two_theta = 2.*theta;
|
||||
return 1;
|
||||
}
|
||||
/*--------------------------------------------------------------------*/
|
||||
int maCalcK(maCrystal data, maAngles angles, double *k){
|
||||
|
||||
*k = ABS(data.dd * Sind(angles.two_theta/2));
|
||||
*k = PI / *k;
|
||||
if(ABS(angles.two_theta/2. - angles.theta) > .1) {
|
||||
return BADSYNC;
|
||||
double maCalcVerticalCurvature(maCrystal data, double two_theta){
|
||||
return calcCurvature(data.VB1,data.VB2, two_theta/2.);
|
||||
}
|
||||
/*-------------------------------------------------------------------*/
|
||||
double maCalcHorizontalCurvature(maCrystal data, double two_theta){
|
||||
return calcCurvature(data.HB1,data.HB2, two_theta/2.);
|
||||
}
|
||||
/*--------------------------------------------------------------------*/
|
||||
double maCalcK(maCrystal data, double two_theta){
|
||||
double k;
|
||||
k = ABS(data.dd * Sind(two_theta/2));
|
||||
if(ABS(k) > .001){
|
||||
k = PI / k;
|
||||
} else {
|
||||
k = .0;
|
||||
}
|
||||
return 1;
|
||||
return k;
|
||||
}
|
||||
/*==================== reciprocal space ==============================*/
|
||||
static MATRIX tasReflectionToHC(tasReflection r, MATRIX B){
|
||||
static MATRIX tasReflectionToHC(tasQEPosition r, MATRIX B){
|
||||
MATRIX h = NULL, hc = NULL;
|
||||
|
||||
h = makeVector();
|
||||
if(h == NULL){
|
||||
return NULL;
|
||||
}
|
||||
vectorSet(h,0,r.h);
|
||||
vectorSet(h,1,r.k);
|
||||
vectorSet(h,2,r.l);
|
||||
vectorSet(h,0,r.qh);
|
||||
vectorSet(h,1,r.qk);
|
||||
vectorSet(h,2,r.ql);
|
||||
|
||||
hc = mat_mul(B,h);
|
||||
killVector(h);
|
||||
@ -135,9 +148,9 @@ static MATRIX uFromAngles(double om, double sgu, double sgl){
|
||||
static MATRIX calcTasUVectorFromAngles(tasReflection r){
|
||||
double theta, om;
|
||||
|
||||
theta = calcTheta(r.ki,r.kf,r.two_theta);
|
||||
om = r.a3 - theta;
|
||||
return uFromAngles(om,r.sgu, r.sgl);
|
||||
theta = calcTheta(r.qe.ki,r.qe.kf,r.angles.sample_two_theta);
|
||||
om = r.angles.a3 - theta;
|
||||
return uFromAngles(om,r.angles.sgu, r.angles.sgl);
|
||||
}
|
||||
/*-------------------------------------------------------------------*/
|
||||
MATRIX calcPlaneNormal(tasReflection r1, tasReflection r2){
|
||||
@ -170,8 +183,8 @@ MATRIX calcTasUBFromTwoReflections(lattice cell, tasReflection r1,
|
||||
*errorCode = status;
|
||||
return NULL;
|
||||
}
|
||||
h1 = tasReflectionToHC(r1,B);
|
||||
h2 = tasReflectionToHC(r2,B);
|
||||
h1 = tasReflectionToHC(r1.qe,B);
|
||||
h2 = tasReflectionToHC(r2.qe,B);
|
||||
if(h1 == NULL || h2 == NULL){
|
||||
*errorCode = UBNOMEMORY;
|
||||
return NULL;
|
||||
@ -265,34 +278,33 @@ static MATRIX buildTVMatrix(MATRIX U1V, MATRIX U2V){
|
||||
return T;
|
||||
}
|
||||
/*-----------------------------------------------------------------------------*/
|
||||
static MATRIX tasReflectionToQC(ptasReflection r, MATRIX UB){
|
||||
static MATRIX tasReflectionToQC(tasQEPosition r, MATRIX UB){
|
||||
MATRIX Q, QC;
|
||||
|
||||
Q = makeVector();
|
||||
if(Q == NULL){
|
||||
return NULL;
|
||||
}
|
||||
vectorSet(Q,0,r->h);
|
||||
vectorSet(Q,1,r->k);
|
||||
vectorSet(Q,2,r->l);
|
||||
vectorSet(Q,0,r.qh);
|
||||
vectorSet(Q,1,r.qk);
|
||||
vectorSet(Q,2,r.ql);
|
||||
QC = mat_mul(UB,Q);
|
||||
killVector(Q);
|
||||
return QC;
|
||||
}
|
||||
/*----------------------------------------------------------------------------*/
|
||||
static MATRIX buildRMatrix(MATRIX UB, MATRIX planeNormal,
|
||||
ptasReflection r){
|
||||
MATRIX U3V, U1V, U2V, TV, TVINV;
|
||||
tasQEPosition qe){
|
||||
MATRIX U1V, U2V, TV, TVINV, M;
|
||||
|
||||
U3V = planeNormal;
|
||||
|
||||
U1V = tasReflectionToQC(r,UB);
|
||||
U1V = tasReflectionToQC(qe,UB);
|
||||
if(U1V == NULL){
|
||||
return NULL;
|
||||
}
|
||||
normalizeVector(U1V);
|
||||
|
||||
U2V = vectorCrossProduct(U3V,U1V);
|
||||
U2V = vectorCrossProduct(planeNormal,U1V);
|
||||
if(U2V == NULL){
|
||||
killVector(U1V);
|
||||
return NULL;
|
||||
@ -305,73 +317,69 @@ static MATRIX buildRMatrix(MATRIX UB, MATRIX planeNormal,
|
||||
return NULL;
|
||||
}
|
||||
TVINV = mat_inv(TV);
|
||||
if(TVINV == NULL){
|
||||
killVector(U1V);
|
||||
killVector(U2V);
|
||||
mat_free(TV);
|
||||
return NULL;
|
||||
}
|
||||
|
||||
|
||||
killVector(U1V);
|
||||
killVector(U2V);
|
||||
mat_free(TVINV);
|
||||
|
||||
mat_free(TV);
|
||||
return TVINV;
|
||||
}
|
||||
/*-------------------------------------------------------------------------------*/
|
||||
int calcTasQAngles(MATRIX UB, MATRIX planeNormal, int ss, ptasReflection r){
|
||||
int calcTasQAngles(MATRIX UB, MATRIX planeNormal, int ss, tasQEPosition qe,
|
||||
ptasAngles angles){
|
||||
MATRIX R, QC;
|
||||
double om, q, theta, cos2t, tmp;
|
||||
double om, q, theta, cos2t, tmp, sq;
|
||||
|
||||
R = buildRMatrix(UB, planeNormal, r);
|
||||
R = buildRMatrix(UB, planeNormal, qe);
|
||||
if(R == NULL){
|
||||
return UBNOMEMORY;
|
||||
}
|
||||
|
||||
om = Acosd(R[0][0]/sqrt(R[0][0]*R[0][0] + R[1][0]*R[1][0]));
|
||||
r->sgl = Acosd(sqrt(R[0][0]*R[0][0] + R[1][0]*R[1][0]));
|
||||
r->sgu = Asind(R[2][1]/sqrt(R[0][0]*R[0][0] + R[1][0]*R[1][0]));
|
||||
|
||||
/*
|
||||
r->sgl = Asind(-R[2][0]);
|
||||
tmp = sqrt(R[0][0]*R[0][0] + R[1][0]*R[1][0]);
|
||||
if(tmp > .001){
|
||||
r->sgu = Acosd(R[2][2]/tmp);
|
||||
} else {
|
||||
sq = sqrt(R[0][0]*R[0][0] + R[1][0]*R[1][0]);
|
||||
if(ABS(sq) < .00001){
|
||||
return BADRMATRIX;
|
||||
}
|
||||
*/
|
||||
|
||||
|
||||
QC = tasReflectionToQC(r,UB);
|
||||
om = Acosd(R[0][0]/sq);
|
||||
om -= 180.;
|
||||
angles->sgl = Acosd(sqrt(R[0][0]*R[0][0] + R[1][0]*R[1][0]));
|
||||
sq = sqrt(R[0][0]*R[0][0] + R[1][0]*R[1][0]);
|
||||
if(ABS(sq) < .00001){
|
||||
return BADRMATRIX;
|
||||
}
|
||||
angles->sgu = Asind(R[2][1]/sq);
|
||||
|
||||
QC = tasReflectionToQC(qe,UB);
|
||||
if(QC == NULL){
|
||||
return UBNOMEMORY;
|
||||
}
|
||||
|
||||
q = vectorLength(QC);
|
||||
q = 2.*PI*vectorLength(QC);
|
||||
cos2t = (r->ki*r->ki + r->kf*r->kf - q*q)/(2. * ABS(r->ki) * ABS(r->kf));
|
||||
cos2t = (qe.ki*qe.ki + qe.kf*qe.kf - q*q)/(2. * ABS(qe.ki) * ABS(qe.kf));
|
||||
if(cos2t > 1.){
|
||||
return TRIANGLENOTCLOSED;
|
||||
}
|
||||
r->two_theta = ss*Acosd(cos2t);
|
||||
angles->sample_two_theta = ss*Acosd(cos2t);
|
||||
|
||||
theta = calcTheta(r->ki, r->kf,r->two_theta);
|
||||
theta = calcTheta(qe.ki, qe.kf,angles->sample_two_theta);
|
||||
|
||||
r->a3 = om + theta;
|
||||
angles->a3 = om + theta;
|
||||
|
||||
killVector(QC);
|
||||
mat_free(R);
|
||||
|
||||
return 1;
|
||||
}
|
||||
/*------------------------------------------------------------------------*/
|
||||
int calcTasQH(MATRIX UB, ptasReflection r){
|
||||
int calcTasQH(MATRIX UB, tasAngles angles, ptasQEPosition qe){
|
||||
MATRIX UBINV = NULL, QV = NULL, Q = NULL;
|
||||
double q;
|
||||
tasReflection r;
|
||||
int i;
|
||||
|
||||
UBINV = mat_inv(UB);
|
||||
QV = calcTasUVectorFromAngles(*r);
|
||||
r.angles = angles;
|
||||
r.qe = *qe;
|
||||
QV = calcTasUVectorFromAngles(r);
|
||||
if(UBINV == NULL || QV == NULL){
|
||||
return UBNOMEMORY;
|
||||
}
|
||||
@ -380,26 +388,23 @@ int calcTasQH(MATRIX UB, ptasReflection r){
|
||||
Thereby take into account the physicists magic fudge
|
||||
2PI factor
|
||||
*/
|
||||
q = sqrt(r->ki*r->ki + r->kf*r->kf - 2.*r->ki*r->kf*Cosd(r->two_theta));
|
||||
q = sqrt(qe->ki*qe->ki + qe->kf*qe->kf - 2.*qe->ki*qe->kf*Cosd(angles.sample_two_theta));
|
||||
q /= 2. * PI;
|
||||
qe->qm = q;
|
||||
|
||||
for(i = 0; i < 3; i++){
|
||||
QV[i][0] *= q;
|
||||
}
|
||||
|
||||
/*
|
||||
mat_dump(UB);
|
||||
mat_dump(UBINV);
|
||||
mat_dump(QV);
|
||||
*/
|
||||
Q = mat_mul(UBINV,QV);
|
||||
if(Q == NULL){
|
||||
mat_free(UBINV);
|
||||
killVector(QV);
|
||||
return UBNOMEMORY;
|
||||
}
|
||||
r->h = Q[0][0];
|
||||
r->k = Q[1][0];
|
||||
r->l = Q[2][0];
|
||||
qe->qh = Q[0][0];
|
||||
qe->qk = Q[1][0];
|
||||
qe->ql = Q[2][0];
|
||||
|
||||
killVector(QV);
|
||||
killVector(Q);
|
||||
@ -413,29 +418,20 @@ int calcAllTasAngles(ptasMachine machine, tasQEPosition qe,
|
||||
int status;
|
||||
tasReflection r;
|
||||
|
||||
status = maCalcAngles(machine->monochromator,&angles->monochromator,
|
||||
qe.ki);
|
||||
status = maCalcTwoTheta(machine->monochromator,qe.ki,
|
||||
&angles->monochromator_two_theta);
|
||||
if(status != 1){
|
||||
return status;
|
||||
}
|
||||
|
||||
r.h = qe.qh;
|
||||
r.k = qe.qk;
|
||||
r.l = qe.ql;
|
||||
r.ki = qe.ki;
|
||||
r.kf = qe.kf;
|
||||
status = calcTasQAngles(machine->UB, machine->planeNormal,
|
||||
machine->ss_sample, &r);
|
||||
machine->ss_sample, qe,angles);
|
||||
if(status != 1){
|
||||
return status;
|
||||
}
|
||||
angles->a3 = r.a3;
|
||||
angles->sample_two_theta = r.two_theta;
|
||||
angles->sgu = r.sgu;
|
||||
angles->sgl = r.sgl;
|
||||
|
||||
status = maCalcAngles(machine->analyzer,&angles->analyzer,
|
||||
qe.kf);
|
||||
status = maCalcTwoTheta(machine->analyzer,qe.kf,&
|
||||
angles->analyzer_two_theta);
|
||||
if(status != 1){
|
||||
return status;
|
||||
}
|
||||
@ -445,47 +441,140 @@ int calcAllTasAngles(ptasMachine machine, tasQEPosition qe,
|
||||
/*----------------------------------------------------------------------*/
|
||||
int calcTasQEPosition(ptasMachine machine, tasAngles angles,
|
||||
ptasQEPosition qe){
|
||||
int status, retVal = 1;
|
||||
tasReflection r;
|
||||
double k;
|
||||
int status;
|
||||
|
||||
status = maCalcK(machine->monochromator,angles.monochromator,&k);
|
||||
if(status != 1){
|
||||
if(status != BADSYNC){
|
||||
retVal = BADSYNC;
|
||||
} else {
|
||||
return status;
|
||||
}
|
||||
}
|
||||
qe->ki = k;
|
||||
qe->ki = maCalcK(machine->monochromator,angles.monochromator_two_theta);
|
||||
qe->kf = maCalcK(machine->analyzer,angles.analyzer_two_theta);
|
||||
|
||||
status = maCalcK(machine->analyzer,angles.analyzer,&k);
|
||||
if(status != 1){
|
||||
if(status != BADSYNC){
|
||||
retVal = BADSYNC;
|
||||
} else {
|
||||
return status;
|
||||
}
|
||||
}
|
||||
qe->kf = k;
|
||||
|
||||
|
||||
r.sgu = angles.sgu;
|
||||
r.sgl = angles.sgl;
|
||||
r.a3 = angles.a3;
|
||||
r.two_theta = angles.sample_two_theta;
|
||||
r.ki = qe->ki;
|
||||
r.kf = qe->kf;
|
||||
status = calcTasQH(machine->UB,&r);
|
||||
status = calcTasQH(machine->UB,angles,qe);
|
||||
if(status != 1){
|
||||
return status;
|
||||
}
|
||||
qe->qh = r.h;
|
||||
qe->qk = r.k;
|
||||
qe->ql = r.l;
|
||||
|
||||
return retVal;
|
||||
return 1;
|
||||
}
|
||||
|
||||
/*================== POWDER Implementation ===========================*/
|
||||
int calcTasPowderAngles(ptasMachine machine, tasQEPosition qe,
|
||||
ptasAngles angles){
|
||||
double cos2t;
|
||||
int status;
|
||||
tasReflection r;
|
||||
|
||||
status = maCalcTwoTheta(machine->monochromator,qe.ki,
|
||||
&angles->monochromator_two_theta);
|
||||
if(status != 1){
|
||||
return status;
|
||||
}
|
||||
|
||||
cos2t = (qe.ki*qe.ki + qe.kf*qe.kf - qe.qm*qe.qm)/(2. * ABS(qe.ki) * ABS(qe.kf));
|
||||
if(cos2t > 1.){
|
||||
return TRIANGLENOTCLOSED;
|
||||
}
|
||||
angles->sample_two_theta = machine->ss_sample*Acosd(cos2t);
|
||||
|
||||
|
||||
status = maCalcTwoTheta(machine->analyzer,qe.kf,&
|
||||
angles->analyzer_two_theta);
|
||||
if(status != 1){
|
||||
return status;
|
||||
}
|
||||
|
||||
return 1;
|
||||
}
|
||||
/*---------------------------------------------------------------------*/
|
||||
int calcTasPowderPosition(ptasMachine machine, tasAngles angles,
|
||||
ptasQEPosition qe){
|
||||
|
||||
int status;
|
||||
|
||||
qe->ki = maCalcK(machine->monochromator,angles.monochromator_two_theta);
|
||||
qe->kf = maCalcK(machine->analyzer,angles.analyzer_two_theta);
|
||||
|
||||
qe->qm = sqrt(qe->ki*qe->ki + qe->kf*qe->kf -
|
||||
2.*qe->ki*qe->kf*Cosd(angles.sample_two_theta));
|
||||
|
||||
return 1;
|
||||
}
|
||||
/*====================== Logic implementation =========================*/
|
||||
void setTasPar(ptasQEPosition qe, int tasMode, int tasVar, double value){
|
||||
double et;
|
||||
|
||||
assert(tasMode == KICONST || tasMode == KFCONST);
|
||||
|
||||
switch(tasVar){
|
||||
case KF:
|
||||
qe->kf = value;
|
||||
break;
|
||||
case EF:
|
||||
qe->kf = energyToK(value);
|
||||
break;
|
||||
case KI:
|
||||
qe->ki = value;
|
||||
break;
|
||||
case EI:
|
||||
qe->ki = energyToK(value);
|
||||
break;
|
||||
case QH:
|
||||
qe->qh = value;
|
||||
break;
|
||||
case QK:
|
||||
qe->qk = value;
|
||||
break;
|
||||
case QL:
|
||||
qe->ql = value;
|
||||
break;
|
||||
case EN:
|
||||
if(tasMode == KICONST){
|
||||
et = KtoEnergy(qe->ki) - value;
|
||||
qe->kf = energyToK(et);
|
||||
} else if(tasMode == KFCONST){
|
||||
et = KtoEnergy(qe->kf) + value;
|
||||
qe->ki = energyToK(et);
|
||||
} else {
|
||||
assert(0);
|
||||
}
|
||||
break;
|
||||
case QM:
|
||||
qe->qm = value;
|
||||
break;
|
||||
default:
|
||||
assert(0);
|
||||
break;
|
||||
}
|
||||
}
|
||||
/*-------------------------------------------------------------------------*/
|
||||
double getTasPar(tasQEPosition qe, int tasVar){
|
||||
switch(tasVar){
|
||||
case EI:
|
||||
return KtoEnergy(qe.ki);
|
||||
break;
|
||||
case KI:
|
||||
return qe.ki;
|
||||
break;
|
||||
case EF:
|
||||
return KtoEnergy(qe.kf);
|
||||
break;
|
||||
case KF:
|
||||
return qe.kf;
|
||||
break;
|
||||
case QH:
|
||||
return qe.qh;
|
||||
break;
|
||||
case QK:
|
||||
return qe.qk;
|
||||
break;
|
||||
case QL:
|
||||
return qe.ql;
|
||||
break;
|
||||
case EN:
|
||||
return KtoEnergy(qe.ki) - KtoEnergy(qe.kf);
|
||||
break;
|
||||
case QM:
|
||||
return qe.qm;
|
||||
break;
|
||||
default:
|
||||
assert(0);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user