- Removed SCStart/EndBuffering as far as possible and fixed an issue with
the capture command in that it not put resluts into the Tcl interpreter. This broke scriptcontext scripts in complicated situations. - Resolved some issues with the TAS calculation and negative scattering sense. - Fixed a bug which did not reset the state to idle after checking reachability in confvirtualmot.c SKIPPED: psi/autowin.c psi/eigera2.c psi/jvlprot.c psi/makefile_linux psi/sinqhttpopt.c psi/tasscan.c
This commit is contained in:
koennecke
252
tasublib.c
252
tasublib.c
@ -202,22 +202,45 @@ static MATRIX uFromAngles(double om, double sgu, double sgl)
|
||||
if (u == NULL) {
|
||||
return NULL;
|
||||
}
|
||||
vectorSet(u, 0, -Cosd(sgl) * Cosd(om));
|
||||
vectorSet(u, 1, Cosd(sgu) * Sind(om) - Sind(sgu) * Sind(sgl) * Cosd(om));
|
||||
vectorSet(u, 0, Cosd(om) * Cosd(sgl));
|
||||
vectorSet(u, 1, -Sind(om) * Cosd(sgu) + Cosd(om)* Sind(sgl)*Sind(sgu));
|
||||
vectorSet(u, 2,
|
||||
-Sind(sgu) * Sind(om) - Cosd(sgu) * Sind(sgl) * Cosd(om));
|
||||
Sind(om) * Sind(sgu) + Cosd(om) * Sind(sgl) * Cosd(sgu));
|
||||
|
||||
return u;
|
||||
}
|
||||
|
||||
/*---------------------------------------------------------------*/
|
||||
MATRIX calcTasUVectorFromAngles(tasReflection r)
|
||||
MATRIX calcTasUVectorFromAngles(tasReflection rr)
|
||||
{
|
||||
double theta, om;
|
||||
double theta, om, ss;
|
||||
MATRIX m;
|
||||
tasReflection r;
|
||||
|
||||
|
||||
r = rr;
|
||||
if(r.angles.sample_two_theta < .0){
|
||||
ss = -1.;
|
||||
} else {
|
||||
ss = 1.;
|
||||
}
|
||||
r.angles.sample_two_theta = ABS(r.angles.sample_two_theta);
|
||||
|
||||
|
||||
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);
|
||||
om = r.angles.a3 - ss*theta;
|
||||
/*
|
||||
This here may have to do with scattering sense.....
|
||||
if (om < -180.){
|
||||
om += 180;
|
||||
}
|
||||
*/
|
||||
m = uFromAngles(om, r.angles.sgu, ss*r.angles.sgl);
|
||||
/*
|
||||
printf("Q = %f, %f, %f inwards Uv = %f %f %f, om = %f\n", r.qe.qh, r.qe.qk, r.qe.ql,
|
||||
m[0][0], m[1][0],m[2][0], om);
|
||||
*/
|
||||
return m ;
|
||||
}
|
||||
|
||||
/*-----------------------------------------------------------------------------*/
|
||||
@ -316,19 +339,179 @@ MATRIX calcPlaneNormal(tasReflection r1, tasReflection r2)
|
||||
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];
|
||||
}
|
||||
}
|
||||
*/
|
||||
mat_free(u1);
|
||||
mat_free(u2);
|
||||
normalizeVector(planeNormal);
|
||||
normalizeVector(planeNormal);
|
||||
return planeNormal;
|
||||
} else {
|
||||
return NULL;
|
||||
}
|
||||
}
|
||||
/*-------------------------------------------------------------------*/
|
||||
MATRIX calcPlaneNormalQOld(MATRIX UB, tasReflection r1, tasReflection r2)
|
||||
{
|
||||
MATRIX u1 = NULL, u2 = NULL, planeNormal = NULL;
|
||||
int i;
|
||||
|
||||
u1 = tasReflectionToQC(r1.qe, UB);
|
||||
u2 = tasReflectionToQC(r2.qe,UB);
|
||||
if (u1 != NULL && u2 != NULL) {
|
||||
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];
|
||||
}
|
||||
}
|
||||
*/
|
||||
mat_free(u1);
|
||||
mat_free(u2);
|
||||
normalizeVector(planeNormal);
|
||||
return planeNormal;
|
||||
} else {
|
||||
return NULL;
|
||||
}
|
||||
}
|
||||
/*-------------------------------------------------------------------*/
|
||||
MATRIX calcPlaneNormalQ(MATRIX UB, tasReflection r1, tasReflection r2)
|
||||
{
|
||||
MATRIX q1 = NULL, q2 = NULL, q3 = NULL, planeNormal = NULL;
|
||||
int i;
|
||||
|
||||
q1 = makeVector();
|
||||
q2 = makeVector();
|
||||
q1[0][0] = r1.qe.qh;
|
||||
q1[1][0] = r1.qe.qk;
|
||||
q1[2][0] = r1.qe.ql;
|
||||
q2[0][0] = r2.qe.qh;
|
||||
q2[1][0] = r2.qe.qk;
|
||||
q2[2][0] = r2.qe.ql;
|
||||
|
||||
q3 = vectorCrossProduct(q1,q2);
|
||||
planeNormal = mat_mul(UB,q3);
|
||||
normalizeVector(planeNormal);
|
||||
if (planeNormal[2][0] < .0) {
|
||||
for (i = 0; i < 3; i++) {
|
||||
planeNormal[i][0] = -1. * planeNormal[i][0];
|
||||
}
|
||||
}
|
||||
mat_free(q1);
|
||||
mat_free(q2);
|
||||
mat_free(q3);
|
||||
|
||||
planeNormal[0][0] *= -1.;
|
||||
planeNormal[1][0] *= -1.;
|
||||
|
||||
return planeNormal;
|
||||
}
|
||||
/*--------------------------------------------------------------------*/
|
||||
MATRIX calcTestUB(lattice cell, double om, double sgu, double sgl)
|
||||
{
|
||||
MATRIX B, M, N, OM, UB;
|
||||
int status;
|
||||
|
||||
/*
|
||||
* create matrices
|
||||
*/
|
||||
B = mat_creat(3, 3, ZERO_MATRIX);
|
||||
|
||||
status = calculateBMatrix(cell, B);
|
||||
if (status < 0) {
|
||||
return NULL;
|
||||
}
|
||||
|
||||
M = mat_creat(3, 3, ZERO_MATRIX);
|
||||
N = mat_creat(3, 3, ZERO_MATRIX);
|
||||
OM = mat_creat(3, 3, ZERO_MATRIX);
|
||||
|
||||
if(M == NULL || N == NULL || OM == NULL){
|
||||
mat_free(B);
|
||||
return NULL;
|
||||
}
|
||||
|
||||
N[0][0] = 1.;
|
||||
N[1][1] = Cosd(sgu);
|
||||
N[1][2] = -Sind(sgu);
|
||||
N[2][1] = Sind(sgu);
|
||||
N[2][2] = Cosd(sgu);
|
||||
|
||||
M[0][0] = Cosd(sgl);
|
||||
M[0][2] = Sind(sgl);
|
||||
M[1][1] = 1.;
|
||||
M[2][0] = -Sind(sgl);
|
||||
M[2][2] = Cosd(sgl);
|
||||
|
||||
OM[0][0] = Cosd(om);
|
||||
OM[0][1] = -Sind(om);
|
||||
OM[1][0] = Sind(om);
|
||||
OM[1][1] = Cosd(om);
|
||||
OM[2][2] = 1.;
|
||||
|
||||
/*
|
||||
Multiply....
|
||||
*/
|
||||
UB = mat_mul(OM,M);
|
||||
UB = mat_mul(UB,N);
|
||||
UB = mat_mul(UB,B);
|
||||
|
||||
mat_free(OM);
|
||||
mat_free(N);
|
||||
mat_free(M);
|
||||
mat_free(B);
|
||||
|
||||
return UB;
|
||||
}
|
||||
/*--------------------------------------------------------------------*/
|
||||
MATRIX calcTestNormal(double sgu, double sgl)
|
||||
{
|
||||
MATRIX M, N, Z, NORM;
|
||||
int status;
|
||||
|
||||
|
||||
M = mat_creat(3, 3, ZERO_MATRIX);
|
||||
N = mat_creat(3, 3, ZERO_MATRIX);
|
||||
Z = mat_creat(3, 1, ZERO_MATRIX);
|
||||
|
||||
|
||||
if(M == NULL || N == NULL || Z == NULL){
|
||||
return NULL;
|
||||
}
|
||||
|
||||
N[0][0] = 1.;
|
||||
N[1][1] = Cosd(sgu);
|
||||
N[1][2] = -Sind(sgu);
|
||||
N[2][1] = Sind(sgu);
|
||||
N[2][2] = Cosd(sgu);
|
||||
|
||||
M[0][0] = Cosd(sgl);
|
||||
M[0][2] = Sind(sgl);
|
||||
M[1][1] = 1.;
|
||||
M[2][0] = -Sind(sgl);
|
||||
M[2][2] = Cosd(sgl);
|
||||
|
||||
Z[2][0] = 1.;
|
||||
|
||||
NORM = mat_inv(N);
|
||||
NORM = mat_mul(NORM,mat_inv(M));
|
||||
NORM = mat_mul(NORM,Z);
|
||||
|
||||
mat_free(N);
|
||||
mat_free(M);
|
||||
mat_free(Z);
|
||||
|
||||
return NORM;
|
||||
}
|
||||
|
||||
/*--------------------------------------------------------------------*/
|
||||
MATRIX calcTasUBFromTwoReflections(lattice cell, tasReflection r1,
|
||||
@ -345,11 +528,18 @@ MATRIX calcTasUBFromTwoReflections(lattice cell, tasReflection r1,
|
||||
calculate the B matrix and the HT matrix
|
||||
*/
|
||||
B = mat_creat(3, 3, ZERO_MATRIX);
|
||||
|
||||
status = calculateBMatrix(cell, B);
|
||||
if (status < 0) {
|
||||
*errorCode = status;
|
||||
return NULL;
|
||||
}
|
||||
/*
|
||||
printf("B-matrix\n");
|
||||
mat_dump(B);
|
||||
printf("\n");
|
||||
*/
|
||||
|
||||
h1 = tasReflectionToHC(r1.qe, B);
|
||||
h2 = tasReflectionToHC(r2.qe, B);
|
||||
if (h1 == NULL || h2 == NULL) {
|
||||
@ -469,6 +659,19 @@ static MATRIX buildTVMatrix(MATRIX U1V, MATRIX U2V)
|
||||
T[i][1] = U2V[i][0];
|
||||
T[i][2] = T3V[i][0];
|
||||
}
|
||||
/*
|
||||
printf("\n");
|
||||
printf("U1V\n");
|
||||
mat_dump(U1V);
|
||||
printf("U2V\n");
|
||||
mat_dump(U2V);
|
||||
printf("U3V\n");
|
||||
mat_dump(T3V);
|
||||
printf("TV-matrix\n");
|
||||
mat_dump(T);
|
||||
printf("\n");
|
||||
*/
|
||||
|
||||
killVector(T3V);
|
||||
return T;
|
||||
}
|
||||
@ -487,8 +690,12 @@ static MATRIX buildRMatrix(MATRIX UB, MATRIX planeNormal,
|
||||
return NULL;
|
||||
}
|
||||
normalizeVector(U1V);
|
||||
/*
|
||||
printf("Q = %f, %f, %f outwards Uv = %f %f %f\n", qe.qh, qe.qk, qe.ql,
|
||||
U1V[0][0], U1V[1][0],U1V[2][0]);
|
||||
*/
|
||||
|
||||
U2V = vectorCrossProduct(planeNormal, U1V);
|
||||
U2V = vectorCrossProduct(planeNormal, U1V);
|
||||
if (U2V == NULL) {
|
||||
killVector(U1V);
|
||||
*errorCode = UBNOMEMORY;
|
||||
@ -514,6 +721,12 @@ static MATRIX buildRMatrix(MATRIX UB, MATRIX planeNormal,
|
||||
*errorCode = BADUBORQ;
|
||||
}
|
||||
|
||||
/*
|
||||
printf("TV after inversion\n");
|
||||
mat_dump(TVINV);
|
||||
printf("\n");
|
||||
*/
|
||||
|
||||
killVector(U1V);
|
||||
killVector(U2V);
|
||||
mat_free(TV);
|
||||
@ -544,7 +757,7 @@ int calcTasQAngles(MATRIX UB, MATRIX planeNormal, int ss, double a3offset,
|
||||
tasQEPosition qe, ptasAngles angles)
|
||||
{
|
||||
MATRIX R, QC;
|
||||
double om, q, theta, cos2t;
|
||||
double om, q, theta, cos2t, cossgl;
|
||||
int errorCode = 1;
|
||||
|
||||
R = buildRMatrix(UB, planeNormal, qe, &errorCode);
|
||||
@ -552,8 +765,10 @@ int calcTasQAngles(MATRIX UB, MATRIX planeNormal, int ss, double a3offset,
|
||||
return errorCode;
|
||||
}
|
||||
|
||||
/* mat_dump(R); */
|
||||
|
||||
angles->sgl = Asind(-R[2][0]);
|
||||
cossgl = sqrt(R[0][0]*R[0][0]+R[1][0]*R[1][0]);
|
||||
angles->sgl = ss*Atan2d(-R[2][0],cossgl);
|
||||
if (ABS(angles->sgl - 90.) < .5) {
|
||||
mat_free(R);
|
||||
return BADUBORQ;
|
||||
@ -573,14 +788,18 @@ int calcTasQAngles(MATRIX UB, MATRIX planeNormal, int ss, double a3offset,
|
||||
R-matrix definition.
|
||||
*/
|
||||
|
||||
om = Atan2d(R[1][0], R[0][0]);
|
||||
angles->sgu = Atan2d(R[2][1], R[2][2]);
|
||||
om = Atan2d(R[1][0]/cossgl, R[0][0]/cossgl);
|
||||
angles->sgu = Atan2d(R[2][1]/cossgl, R[2][2]/cossgl);
|
||||
|
||||
QC = tasReflectionToQC(qe, UB);
|
||||
if (QC == NULL) {
|
||||
mat_free(R);
|
||||
return UBNOMEMORY;
|
||||
}
|
||||
|
||||
/*
|
||||
printf("Q = %f, %f, %f calculated om = %f\n", qe.qh, qe.qk, qe.ql, om);
|
||||
*/
|
||||
|
||||
q = vectorLength(QC);
|
||||
q = 2. * PI * vectorLength(QC);
|
||||
@ -592,20 +811,23 @@ int calcTasQAngles(MATRIX UB, MATRIX planeNormal, int ss, double a3offset,
|
||||
killVector(QC);
|
||||
return TRIANGLENOTCLOSED;
|
||||
}
|
||||
theta = calcTheta(qe.ki, qe.kf, Acosd(cos2t));
|
||||
angles->sample_two_theta = ss * Acosd(cos2t);
|
||||
|
||||
theta = calcTheta(qe.ki, qe.kf, angles->sample_two_theta);
|
||||
|
||||
angles->a3 = om + theta + a3offset;
|
||||
angles->a3 = om + ss*theta + a3offset;
|
||||
/*
|
||||
put a3 into -180, 180 properly. We can 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 - a3offset) {
|
||||
angles->a3 += 360.;
|
||||
}
|
||||
*/
|
||||
angles->a3 = fmod(angles->a3 + ss*180.,360.) - ss*180.;
|
||||
|
||||
killVector(QC);
|
||||
mat_free(R);
|
||||
|
||||
Reference in New Issue
Block a user