linse/sics
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

- Fixed sign problems with om in tasub

Browse Source 
- Mended tasdrive to drive energy even if Q is askew
- Fixed QM values
- Fixed problems in mesure: om2th, wrong theta selection
- Fixed core dump when driving h,kl, failed
This commit is contained in:
koennecke
2005-06-09 12:04:38 +00:00
parent f33ca7b0d7
commit 3a0b4f293c
16 changed files with 1899 additions and 116 deletions
Download Patch File Download Diff File Expand all files Collapse all files

159
tasublib.c
View File

@ -18,6 +18,8 @@
#define PI 3.141592653589793
#define ECONST 2.072
#define DEGREE_RAD (PI/180.0) /* Radians per degree */
#define VERT 0
#define HOR 1
/*============== monochromator/analyzer stuff =========================*/
double energyToK(double energy){
double K;
@ -33,8 +35,14 @@ double KtoEnergy(double k){
return energy;
}
/*-------------------------------------------------------------------*/
static double calcCurvature(double B1, double B2, double theta){
return B1 + B2/Sind(ABS(theta));
static double calcCurvature(double B1, double B2, double theta,
int ori){
assert(ori == VERT || ori == HOR);
if(ori == VERT){
return B1 + B2/Sind(ABS(theta));
} else {
return B1 + B2*Sind(ABS(theta));
}
}
/*--------------------------------------------------------------------*/
int maCalcTwoTheta(maCrystal data, double k, double *two_theta){
@ -51,11 +59,11 @@ int maCalcTwoTheta(maCrystal data, double k, double *two_theta){
}
/*--------------------------------------------------------------------*/
double maCalcVerticalCurvature(maCrystal data, double two_theta){
return calcCurvature(data.VB1,data.VB2, two_theta/2.);
return calcCurvature(data.VB1,data.VB2, two_theta/2.,VERT);
}
/*-------------------------------------------------------------------*/
double maCalcHorizontalCurvature(maCrystal data, double two_theta){
return calcCurvature(data.HB1,data.HB2, two_theta/2.);
return calcCurvature(data.HB1,data.HB2, two_theta/2.,HOR);
}
/*--------------------------------------------------------------------*/
double maCalcK(maCrystal data, double two_theta){
@ -155,11 +163,22 @@ static MATRIX calcTasUVectorFromAngles(tasReflection r){
/*-------------------------------------------------------------------*/
MATRIX calcPlaneNormal(tasReflection r1, tasReflection r2){
MATRIX u1 = NULL, u2 = NULL, planeNormal = NULL;
int i;
u1 = calcTasUVectorFromAngles(r1);
u2 = calcTasUVectorFromAngles(r2);
if(u1 != NULL && u2 != NULL){
return vectorCrossProduct(u1,u2);
planeNormal = vectorCrossProduct(u1,u2);
/*
The plane normal has to point to the stars and not to the earth
core in order for the algorithm to work.
*/
if(planeNormal[2][0] < .0){
for(i = 0; i < 3; i++){
planeNormal[i][0] = -1.*planeNormal[i][0];
}
}
return planeNormal;
} else {
return NULL;
}
@ -187,11 +206,15 @@ MATRIX calcTasUBFromTwoReflections(lattice cell, tasReflection r1,
h2 = tasReflectionToHC(r2.qe,B);
if(h1 == NULL || h2 == NULL){
*errorCode = UBNOMEMORY;
mat_free(B);
return NULL;
}
HT = matFromTwoVectors(h1,h2);
if(HT == NULL){
*errorCode = UBNOMEMORY;
mat_free(B);
mat_free(h1);
mat_free(h2);
return NULL;
}
@ -202,39 +225,59 @@ MATRIX calcTasUBFromTwoReflections(lattice cell, tasReflection r1,
u2 = calcTasUVectorFromAngles(r2);
if(u1 == NULL || u2 == NULL){
*errorCode = UBNOMEMORY;
mat_free(B);
mat_free(h1);
mat_free(h2);
return NULL;
}
UT = matFromTwoVectors(u1,u2);
if(UT == NULL){
*errorCode = UBNOMEMORY;
mat_free(B);
mat_free(h1);
mat_free(h2);
mat_free(u1);
mat_free(u2);
mat_free(HT);
return NULL;
}
/*
debugging output
printf("B-matrix\n");
mat_dump(B);
printf("HT-matrix\n");
mat_dump(HT);
printf("UT-matrix\n");
mat_dump(UT);
*/
/*
UT = U * HT
*/
HTT = mat_tran(HT);
if(HTT == NULL){
*errorCode = UBNOMEMORY;
mat_free(B);
mat_free(h1);
mat_free(h2);
mat_free(u1);
mat_free(u2);
mat_free(HT);
return NULL;
}
U = mat_mul(UT,HTT);
if(U == NULL){
*errorCode = UBNOMEMORY;
mat_free(B);
mat_free(h1);
mat_free(h2);
mat_free(u1);
mat_free(u2);
mat_free(HT);
mat_free(HTT);
return NULL;
}
UB = mat_mul(U,B);
if(UB == NULL){
mat_free(B);
mat_free(h1);
mat_free(h2);
mat_free(u1);
mat_free(u2);
mat_free(HT);
mat_free(HTT);
mat_free(U);
*errorCode = UBNOMEMORY;
}
@ -296,12 +339,14 @@ static MATRIX tasReflectionToQC(tasQEPosition r, MATRIX UB){
}
/*----------------------------------------------------------------------------*/
static MATRIX buildRMatrix(MATRIX UB, MATRIX planeNormal,
tasQEPosition qe){
tasQEPosition qe, int *errorCode){
MATRIX U1V, U2V, TV, TVINV, M;
*errorCode = 1;
U1V = tasReflectionToQC(qe,UB);
if(U1V == NULL){
*errorCode = UBNOMEMORY;
return NULL;
}
normalizeVector(U1V);
@ -309,6 +354,13 @@ static MATRIX buildRMatrix(MATRIX UB, MATRIX planeNormal,
U2V = vectorCrossProduct(planeNormal,U1V);
if(U2V == NULL){
killVector(U1V);
*errorCode = UBNOMEMORY;
return NULL;
}
if(vectorLength(U2V) < .0001){
*errorCode = BADUBORQ;
killVector(U1V);
killVector(U2V);
return NULL;
}
@ -316,9 +368,14 @@ static MATRIX buildRMatrix(MATRIX UB, MATRIX planeNormal,
if(TV == NULL){
killVector(U1V);
killVector(U2V);
*errorCode = UBNOMEMORY;
return NULL;
}
TVINV = mat_inv(TV);
if(TVINV == NULL){
*errorCode = BADUBORQ;
}
killVector(U1V);
killVector(U2V);
@ -329,30 +386,36 @@ static MATRIX buildRMatrix(MATRIX UB, MATRIX planeNormal,
int calcTasQAngles(MATRIX UB, MATRIX planeNormal, int ss, tasQEPosition qe,
ptasAngles angles){
MATRIX R, QC;
double om, q, theta, cos2t, tmp, sq;
double om, q, theta, cos2t;
int errorCode = 1;
R = buildRMatrix(UB, planeNormal, qe);
R = buildRMatrix(UB, planeNormal, qe, &errorCode);
if(R == NULL){
return UBNOMEMORY;
return errorCode;
}
sq = sqrt(R[0][0]*R[0][0] + R[1][0]*R[1][0]);
if(ABS(sq) < .00001){
return BADRMATRIX;
angles->sgl = Asind(-R[2][0]);
if(ABS(angles->sgl -90.) < .5){
return BADUBORQ;
}
om = Acosd(R[0][0]/sq);
om -= 180.;
tmp = Asind(R[1][0]/sqrt(R[0][0]*R[0][0] + R[1][0]*R[1][0]));
/*
Now, this is slightly different then in the publication by M. Lumsden.
The reason is that the atan2 helps to determine the sign of om
whereas the sin, cos formula given by M. Lumsden yield ambiguous signs
especially for om.
sgu = atan(R[2][1],R[2][2]) where:
R[2][1] = cos(sgl)sin(sgu)
R[2][2] = cos(sgu)cos(sgl)
om = atan(R[1][0],R[0][0]) where:
R[1][0] = sin(om)cos(sgl)
R[0][0] = cos(om)cos(sgl)
The definitions of th R components are taken from M. Lumsden
R-matrix definition.
*/
tmp = Acosd(sqrt(R[0][0]*R[0][0] + R[1][0]*R[1][0]));
angles->sgl = Asind(-R[2][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);
tmp = Acosd(R[2][2]/sqrt(R[0][0]*R[0][0]+R[1][0]*R[1][0]));
om = Atan2d(R[1][0],R[0][0]);
angles->sgu = Atan2d(R[2][1],R[2][2]);
QC = tasReflectionToQC(qe,UB);
if(QC == NULL){
@ -370,6 +433,15 @@ int calcTasQAngles(MATRIX UB, MATRIX planeNormal, int ss, tasQEPosition qe,
theta = calcTheta(qe.ki, qe.kf,angles->sample_two_theta);
angles->a3 = om + theta;
/*
put a3 into -180, 180 properly. We cal always turn by 180 because the
scattering geometry is symmetric in this respect. It is like looking at
the scattering plane from the other side
*/
angles->a3 -= 180.;
if(angles->a3 < -180.){
angles->a3 += 360.;
}
killVector(QC);
mat_free(R);
@ -395,9 +467,10 @@ int calcTasQH(MATRIX UB, tasAngles angles, ptasQEPosition qe){
Thereby take into account the physicists magic fudge
2PI factor
*/
q = sqrt(qe->ki*qe->ki + qe->kf*qe->kf - 2.*qe->ki*qe->kf*Cosd(angles.sample_two_theta));
q /= 2. * PI;
q = sqrt(qe->ki*qe->ki + qe->kf*qe->kf -
2.*qe->ki*qe->kf*Cosd(angles.sample_two_theta));
qe->qm = q;
q /= 2. * PI;
for(i = 0; i < 3; i++){
QV[i][0] *= q;
@ -431,14 +504,14 @@ int calcAllTasAngles(ptasMachine machine, tasQEPosition qe,
return status;
}
status = calcTasQAngles(machine->UB, machine->planeNormal,
machine->ss_sample, qe,angles);
status = maCalcTwoTheta(machine->analyzer,qe.kf,&
angles->analyzer_two_theta);
if(status != 1){
return status;
}
status = maCalcTwoTheta(machine->analyzer,qe.kf,&
angles->analyzer_two_theta);
status = calcTasQAngles(machine->UB, machine->planeNormal,
machine->ss_sample, qe,angles);
if(status != 1){
return status;
}
@ -498,7 +571,6 @@ int calcTasPowderPosition(ptasMachine machine, tasAngles angles,
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 =========================*/
@ -582,6 +654,3 @@ double getTasPar(tasQEPosition qe, int tasVar){
assert(0);
}
}