linse/sics
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
main
sics/fourlib.c
Ferdi Franceschini 10d29d597c Cleaned up ANSTO code to merge with sinqdev.sics 
This is our new RELEASE-4_0 branch which was taken from ansto/93d9a7c
Conflicts:
	.gitignore
	SICSmain.c
	asynnet.c
	confvirtualmot.c
	counter.c
	devexec.c
	drive.c
	event.h
	exebuf.c
	exeman.c
	histmem.c
	interface.h
	motor.c
	motorlist.c
	motorsec.c
	multicounter.c
	napi.c
	napi.h
	napi4.c
	network.c
	nwatch.c
	nxscript.c
	nxxml.c
	nxxml.h
	ofac.c
	reflist.c
	scan.c
	sicshipadaba.c
	sicsobj.c
	site_ansto/docs/Copyright.txt
	site_ansto/instrument/lyrebird/config/tasmad/sicscommon/nxsupport.tcl
	site_ansto/instrument/lyrebird/config/tasmad/taspub_sics/tasscript.tcl
	statusfile.c
	tasdrive.c
	tasub.c
	tasub.h
	tasublib.c
	tasublib.h
2015-04-23 20:49:26 +10:00

944 lines
24 KiB
C
Raw Permalink Blame History

/*
F O U R L I B
This is a library of routines for doing transformations between the
various coordinate systems used on a four circle diffractometer as
used for neutron or X-ray diffraction. The coordinate systems used are
described in Busing, Levy, Acta Cryst (1967),22, 457 ff.
Generally we have:
Z = [OM][CHI][PHI][UB]h
where: Z is a vector in the diffractometer coordinate system
OM CHI PHI are rotation matrices around the respective angles
UB is the UB matrix
h is the reciprocal lattice vector.
The vector Z cannot only be expressed in terms of the angles stt, om, chi,
and phi put also by polar coordinates gamma and nu.
This code is a reimplementation based on a F77 code from Gary McIntyre, ILL
and code extracted from the ILL MAD control program.
Mark Koennecke, November-December 2001
*/
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "fortify.h"
#include "fourlib.h"
#define PI 3.141592653589793
#define RD 57.30
#define ABS(x) (x < 0 ? -(x) : (x))
/*------------------------------------------------------------------------
a safe routine for calculating the atan of y/x
------------------------------------------------------------------------*/
static double myatan(double y, double x)
{
if (ABS(x) < 0.0001) {
if (y > .0) {
return PI / 2;
} else {
return -PI / 2;
}
}
if (x > .0) {
return atan(y / x);
} else {
if (y > .0) {
return PI + atan(y / x);
} else {
return -PI + atan(y / x);
}
}
}
/*-------------------------------------------------------------------------*/
static double rtan(double y, double x)
{
double val;
if ((x == 0.) && (y == 0.)) {
return .0;
}
if (x == 0.) {
if (y < 0.) {
return -PI / 2.;
} else {
return PI / 2.;
}
}
if (ABS(y) < ABS(x)) {
val = atan(ABS(y / x));
if (x < 0.) {
val = PI - val;
}
if (y < 0.) {
val = -val;
}
return val;
} else {
val = PI / 2. - atan(ABS(x / y));
if (x < 0.) {
val = PI - val;
}
if (y < 0.) {
val = -val;
}
}
return val;
}
/*----------------------------------------------------------------------
clear3x3 sets a 3x3 matrix to 0
-----------------------------------------------------------------------*/
static void clear3x3(MATRIX target)
{
int i, j;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
target[i][j] = .0;
}
}
}
/*---------------------------------------------------------------------*/
MATRIX vectorToMatrix(double z[3])
{
int i;
MATRIX res;
res = mat_creat(3, 1, ZERO_MATRIX);
for (i = 0; i < 3; i++) {
res[i][0] = z[i];
}
return res;
}
/*---------------------------------------------------------------------*/
void matrixToVector(MATRIX zm, double z[3])
{
int i;
for (i = 0; i < 3; i++) {
z[i] = zm[i][0];
}
}
/*----------------------------------------------------------------------
A Busing & levy PSI matrix
----------------------------------------------------------------------*/
static void psimat(MATRIX target, double psi)
{
double psird = psi / RD;
clear3x3(target);
target[0][0] = 1.;
target[1][1] = cos(psird);
target[1][2] = sin(psird);
target[2][1] = -target[1][2];
target[2][2] = target[1][1];
}
/*----------------------------------------------------------------------
A Busing & levy CHI matrix
----------------------------------------------------------------------*/
void chimat(MATRIX target, double chi)
{
double chird = chi / RD;
clear3x3(target);
target[0][0] = cos(chird);
target[0][2] = sin(chird);
target[1][1] = 1.;
target[2][0] = -target[0][2];
target[2][2] = target[0][0];
}
/*----------------------------------------------------------------------
A Busing & levy PHI matrix
----------------------------------------------------------------------*/
void phimat(MATRIX target, double phi)
{
double phird = phi / RD;
clear3x3(target);
target[0][0] = cos(phird);
target[0][1] = sin(phird);
target[1][0] = -target[0][1];
target[1][1] = target[0][0];
target[2][2] = 1.;
}
/*------------------- PSD specific code ----------------------------*/
void det2pol(psdDescription * psd, int x, int y, double *gamma, double *nu)
{
double xobs, yobs, b, z, d;
assert(ABS(psd->distance) > .001);
xobs = (x - psd->xZero) * psd->xScale;
yobs = (y - psd->yZero) * psd->yScale;
b = psd->distance * cos(yobs / psd->distance);
z = psd->distance * sin(yobs / psd->distance);
d = sqrt(xobs * xobs + b * b);
*gamma = psd->gamma + myatan(xobs, b) * RD;
*nu = psd->nu + myatan(z, d) * RD;
}
/*----------------------------------------------------------------------*/
void pol2det(psdDescription * psd, double gamma, double nu, int *x, int *y)
{
double delga, delnu, td, tn, e, f, g, zobs, xobs;
delga = gamma - psd->gamma;
delnu = nu - psd->nu;
td = tan(delga / RD);
tn = tan(delnu / RD);
e = sqrt(1. + td * td);
f = tn * e;
g = psd->distance * (1. + tn * tn + tn * tn * td * td);
zobs = psd->distance * (atan(tn * e) + asin(f / g));
xobs = td * (psd->distance * cos(zobs / psd->distance));
*y = (int) rint(psd->yZero + zobs / psd->yScale);
*x = (int) rint(psd->xZero + xobs / psd->xScale);
}
/*--------------- bisecting geometry code -----------------------------*/
/*
turn chi and phi in order to get Z1 into the equatorial plane
*/
static void turnEquatorial(MATRIX z1, double *chi, double *phi)
{
if (ABS(z1[0][0]) < .0001 && ABS(z1[1][0]) < .00001) {
*phi = .0;
*chi = 90.;
if (z1[2][0] < .0) {
*chi = -*chi;
}
} else {
*phi = myatan(z1[1][0], z1[0][0]) * RD;
*chi = myatan(z1[2][0],
sqrt(z1[0][0] * z1[0][0] + z1[1][0] * z1[1][0])) * RD;
}
}
/*----------------------------------------------------------------------
calculate d-spacing and theta from z1
-----------------------------------------------------------------------*/
int calcTheta(double lambda, MATRIX z1, double *d, double *theta)
{
double dstar, sintheta;
dstar =
sqrt(z1[0][0] * z1[0][0] + z1[1][0] * z1[1][0] +
z1[2][0] * z1[2][0]);
if (dstar < .0001) {
*d = .0;
*theta = 0;
return 0;
}
*d = 1. / dstar;
sintheta = lambda * dstar / 2.;
if (ABS(sintheta) > 1.0) {
*d = .0;
*theta = .0;
return 0;
}
*theta = asin(sintheta) * RD;
return 1;
}
/*-------------------------------------------------------------------*/
int z1ToBisecting(double lambda, double z1[3], double *stt, double *om,
double *chi, double *phi)
{
MATRIX zz1;
int status;
zz1 = vectorToMatrix(z1);
status = z1mToBisecting(lambda, zz1, stt, om, chi, phi);
mat_free(zz1);
return status;
}
/*------------------------------------------------------------------------*/
int z1mToBisecting(double lambda, MATRIX z1, double *stt, double *om,
double *chi, double *phi)
{
double d;
int status;
status = calcTheta(lambda, z1, &d, om);
if (!status) {
*stt = .0;
*om = .0;
*chi = .0;
*phi = .0;
return status;
}
turnEquatorial(z1, chi, phi);
*stt = *om * 2.;
*chi = 180. - *chi;
*phi = 180. + *phi;
return 1;
}
/*---------------------------------------------------------------------*/
int z1ToAnglesWithOffset(double lambda, MATRIX z1m, double omOffset,
double *stt, double *om, double *chi, double *phi)
{
int status, i;
double d, delOm, sinchi, q, cosp, sinp, dstar;
MATRIX z1;
status = calcTheta(lambda, z1m, &d, om);
if (!status) {
*stt = .0;
*om = .0;
*chi = .0;
*phi = .0;
return status;
}
*stt = 2. * *om;
*om += omOffset;
z1 = mat_copy(z1m); /* routine messes z1m up! */
dstar =
sqrt(z1[0][0] * z1[0][0] + z1[1][0] * z1[1][0] +
z1[2][0] * z1[2][0]);
for (i = 0; i < 3; i++) {
z1[i][0] /= dstar;
}
delOm = -omOffset / RD;
sinchi = z1[2][0] / cos(delOm);
if (ABS(sinchi) > 1.) {
mat_free(z1);
return 0;
}
*chi = asin(sinchi);
*chi = PI - *chi;
if (*chi > PI) {
*chi = *chi - 2 * PI;
}
q = cos(delOm) * cos(*chi);
cosp = (sin(delOm) * z1[0][0] + q * z1[1][0]) /
(sin(delOm) * sin(delOm) + q * q);
sinp = (z1[0][0] - sin(delOm) * cosp) / q;
*phi = rtan(cosp, sinp);
*phi *= RD;
*chi *= RD;
mat_free(z1);
return 1;
}
/*---------------------------------------------------------------------*/
int psiForOmegaOffset(MATRIX z1m, double omOffset,
double chi, double phi, double *psi)
{
double chibi, sinps, cosps, sinchi, dstar;
MATRIX z1;
z1 = mat_copy(z1m); /* routine messes z1m up! */
dstar =
sqrt(z1[0][0] * z1[0][0] + z1[1][0] * z1[1][0] +
z1[2][0] * z1[2][0]);
z1[2][0] /= dstar;
omOffset /= RD;
if (-z1[2][0] == 1.) {
*psi = phi;
} else {
chibi = PI - asin(-z1[2][0]);
if (chibi > PI) {
chibi -= 2 * PI;
}
sinps = tan(chibi) * tan(-omOffset);
sinchi = -z1[2][0] / cos(-omOffset);
if (ABS(sinchi) > 1.) {
mat_free(z1);
return 0;
}
cosps = (cos(chibi) - sinps * sin(-omOffset) * sinchi) / cos(chi / RD);
*psi = -rtan(sinps, cosps) * RD;
}
mat_free(z1);
return 1;
}
/*----------------------------------------------------------------------*/
void rotatePsi(double om, double chi, double phi, double psi,
double *newom, double *newchi, double *newphi)
{
MATRIX chim, phim, psim, r0, r0psi;
chim = mat_creat(3, 3, ZERO_MATRIX);
chimat(chim, chi);
phim = mat_creat(3, 3, ZERO_MATRIX);
phimat(phim, phi);
r0 = mat_mul(chim, phim);
psim = mat_creat(3, 3, ZERO_MATRIX);
psimat(psim, psi);
r0psi = mat_mul(psim, r0);
*newchi =
rtan(sqrt(r0psi[2][0] * r0psi[2][0] + r0psi[2][1] * r0psi[2][1]),
r0psi[2][2]) * RD;
*newphi = rtan(-r0psi[2][1], -r0psi[2][0]) * RD;
*newom = om + rtan(-r0psi[1][2], r0psi[0][2]) * RD;
mat_free(chim);
mat_free(phim);
mat_free(psim);
mat_free(r0);
mat_free(r0psi);
}
/*--------------------------------------------------------------------
calculate Z1 = [PHI]TZ3
The T means transposed matrixes for the return calculation!
---------------------------------------------------------------------*/
static void z1fromz2(double z1[3], MATRIX z2, double phi)
{
MATRIX phim, phimt, zz1;
phim = mat_creat(3, 3, ZERO_MATRIX);
phimat(phim, phi);
phimt = mat_tran(phim);
zz1 = mat_mul(phimt, z2);
matrixToVector(zz1, z1);
mat_free(phim);
mat_free(phimt);
mat_free(zz1);
}
/*--------------------------------------------------------------------
calculate Z1 = [PHI]T*[CHI]T*Z3
The T means transposed matrixes for the return calculation!
---------------------------------------------------------------------*/
static void z1fromz3(double z1[3], MATRIX z3, double chi, double phi)
{
MATRIX chim, chimt, z2;
chim = mat_creat(3, 3, ZERO_MATRIX);
chimat(chim, chi);
chimt = mat_tran(chim);
z2 = mat_mul(chimt, z3);
z1fromz2(z1, z2, phi);
mat_free(chim);
mat_free(chimt);
mat_free(z2);
}
/*--------------------------------------------------------------------
calculate Z1 = [PHI]T*[CHI]T*[OM]T*Z4
The T means transposed matrixes for the return calculation!
---------------------------------------------------------------------*/
static void z1fromz4(double z1[3], MATRIX z4, double om, double chi,
double phi)
{
MATRIX oma, oma2, z3;
oma = mat_creat(3, 3, ZERO_MATRIX);
phimat(oma, om);
oma2 = mat_tran(oma);
z3 = mat_mul(oma2, z4);
z1fromz3(z1, z3, chi, phi);
mat_free(oma);
mat_free(oma2);
mat_free(z3);
}
/*---------------------------------------------------------------------*/
void z1FromAngles(double lambda, double stt, double om,
double chi, double phi, double z1[3])
{
MATRIX z4;
double th;
z4 = mat_creat(3, 1, ZERO_MATRIX);
th = (stt / 2.) / RD;
z4[0][0] = (2. * sin(th) * cos(th)) / lambda;
z4[1][0] = (-2. * sin(th) * sin(th)) / lambda;
z4[2][0] = .0;
z1fromz4(z1, z4, om, chi, phi);
mat_free(z4);
}
/*----------------------------------------------------------------------
Normal Beam geometry specific calculations. This means the reflection is
expressed through the three angles omega, gamma (two theta detector) and
nu (detector tilt).
bisToNormalBeam was lifted from HKLGEN.
-----------------------------------------------------------------------*/
double sign(double a, double b)
{
if (b >= .0) {
return ABS(a);
} else {
return -ABS(a);
}
}
/*--------------------------------------------------------------------*/
static void makeNull(double *gamma, double *om, double *nu)
{
*gamma = .0;
*om = .0;
*nu = .0;
}
/*---------------------------------------------------------------------*/
int z1mToNormalBeam(double lambda, MATRIX z1m, double *gamma, double *om,
double *nu)
{
MATRIX dum, znew;
double d, a, b, sint, theta, omdeg;
int status;
status = calcTheta(lambda, z1m, &d, &theta);
if (!status) {
makeNull(gamma, om, nu);
return status;
}
/* Everything on omega axis is blind: test for this */
a = sqrt(z1m[0][0] * z1m[0][0] + z1m[1][0] * z1m[1][0]);
if (ABS(a) < .0001) {
makeNull(gamma, om, nu);
return 0;
}
sint = sin(theta / RD);
b = 2. * sint * sint / (lambda * a);
if (b >= 1.) {
makeNull(gamma, om, nu);
return 0;
}
a = -atan2(z1m[1][0], -z1m[0][0]);
b = -asin(b);
*om = a + b;
omdeg = *om * RD;
dum = mat_creat(3, 3, ZERO_MATRIX);
phimat(dum, omdeg);
znew = mat_mul(dum, z1m);
if (znew[0][0] < 0) {
*om = *om - 2. * atan2(-znew[0][0], -znew[1][0]);
omdeg = *om * RD;
}
b = (sign(180., omdeg) + omdeg) / 360.;
b = sign(1, b) * floor(ABS(b));
omdeg = omdeg - 360. * b;
*nu = asin(lambda * z1m[2][0]);
*gamma = acos(cos(2. * (theta / RD)) / cos(*nu));
*om = omdeg;
*nu = *nu * RD;
*gamma = *gamma * RD;
mat_free(dum);
mat_free(znew);
return 1;
}
/*----------------------------------------------------------------------*/
int bisToNormalBeam(double twotheta, double omega, double chi, double phi,
double *omeganb, double *gamma, double *nu)
{
double tth, fom, fchi, fphi, sint, sinnu, del;
tth = twotheta / RD;
fom = omega / RD;
fchi = chi / RD;
fphi = phi / RD;
/* nu calculation */
sint = sin(tth / 2.);
*nu = 2. * sin(fchi) * sint;
if (*nu > 1.) {
return 0;
}
*nu = asin(*nu);
*nu = ABS(*nu) * sign(1., fchi);
sinnu = sin(*nu);
/* gamma calculation */
*gamma = cos(tth) / cos(*nu);
if (ABS(*gamma) > 1.0) {
return 0;
}
*gamma = acos(*gamma);
*gamma = *gamma * sign(1., tth);
/* omega normal beam */
del = asin(sint / cos(fchi));
*omeganb = del + fphi;
*omeganb *= RD;
*gamma *= RD;
*nu *= RD;
if (ABS(*omeganb) > 180.) {
if (*omeganb < -180.) {
*omeganb += 360.;
}
if (*omeganb > 180.) {
*omeganb -= 360.;
}
}
return 1;
}
/* --------------------------------------------------------------------*/
static void z4FromNormalBeam(MATRIX z4, double lambda, double gamma,
double nu)
{
double gar, nur;
gar = gamma / RD;
nur = nu / RD;
/*
this is basically a conversion from polar coordinates to
x,y,z coordinates. However, sin and cos are exchanged in
comparison to textbook formulas for this conversion. But this
is only a shift of origin
*/
z4[0][0] = (sin(gar) * cos(nur)) / lambda;
z4[1][0] = (cos(gar) * cos(nur) - 1.) / lambda;
z4[2][0] = (sin(nur)) / lambda;
}
/*----------------------------------------------------------------------*/
void z1FromNormalBeam(double lambda, double omega, double gamma,
double nu, double z1[3])
{
MATRIX z4, omm, omt, z1m;
z4 = mat_creat(3, 1, ZERO_MATRIX);
z4FromNormalBeam(z4, lambda, gamma, nu);
omm = mat_creat(3, 3, ZERO_MATRIX);
phimat(omm, omega);
omt = mat_tran(omm);
z1m = mat_mul(omt, z4);
matrixToVector(z1m, z1);
mat_free(z4);
mat_free(omm);
mat_free(omt);
mat_free(z1m);
}
/*--------------------------------------------------------------------*/
double circlify(double val)
{
while (val >= 359.8) {
val -= 360.;
}
while (val < 0.) {
val += 360.;
}
return val;
}
/*----------------------------------------------------------------------*/
void z1FromAllAngles(double lambda, double omega, double gamma,
double nu, double chi, double phi, double z1[3])
{
MATRIX z3;
z1FromNormalBeam(lambda, omega, gamma, nu, z1);
z3 = vectorToMatrix(z1);
z1fromz3(z1, z3, chi, phi);
mat_free(z3);
}
/*-----------------------------------------------------------------------*/
int findAllowedBisecting(double lambda, MATRIX z1, double *fSet,
inRange testFunc, void *userData)
{
int status, i, mask[4];
double stt, om, chi, phi, psi, ompsi, chipsi, phipsi;
double fTest[4];
status = z1mToBisecting(lambda, z1, &stt, &om, &chi, &phi);
chi = circlify(chi);
phi = circlify(phi);
fSet[0] = stt;
fSet[1] = om;
fSet[2] = chi;
fSet[3] = phi;
if (status != 1) {
return 0;
}
if (testFunc(userData, fSet, mask) == 1) {
return 1;
}
for (psi = .0; psi < 360.; psi += .5) {
rotatePsi(om, chi, phi, psi, &ompsi, &chipsi, &phipsi);
fTest[0] = stt;
fTest[1] = ompsi;
fTest[2] = circlify(chipsi);
fTest[3] = circlify(phipsi);
status = testFunc(userData, fTest, mask);
if (status == 1) {
for (i = 0; i < 4; i++) {
fSet[i] = fTest[i];
}
return 1;
}
/*
* if chi close to 0, or 180, try to wrap phi onto om
*/
if (ABS(fTest[2] - .0) < .1 || ABS(fTest[2] - 180.) < .1) {
fTest[1] -= fTest[3];
fTest[3] = .0;
if (fTest[1] < 0.) {
fTest[1] += 360.;
}
if (fTest[1] > 360.0) {
fTest[1] -= 360.;
}
status = testFunc(userData, fTest, mask);
if (status == 1) {
for (i = 0; i < 4; i++) {
fSet[i] = fTest[i];
}
return 1;
}
}
if (mask[0] == 0) {
/*
* useless: when two theta problem there is no solution
*/
return 0;
}
}
return 0;
}
/*------------------- a test program ------------------------------------*/
#ifdef TESTCODE
int main(int argc, char *argv[])
{
psdDescription psd;
double gamma, nu, lambda, stt, om, chi, phi, z1[3], psi;
double psiom, psichi, psiphi;
double sttex, omex, chiex, phiex;
int x, y, status, i;
MATRIX UB, UBinv, zz1, zz1b, hkl, hklback;
double gaex, omnbex, nuex, omeganb;
/* initializations */
UB = mat_creat(3, 3, ZERO_MATRIX);
UB[0][0] = 0.016169;
UB[0][1] = 0.011969;
UB[0][2] = 0.063195;
UB[1][0] = -.000545;
UB[1][1] = 0.083377;
UB[1][2] = -0.009117;
UB[2][0] = -0.162051;
UB[2][1] = 0.000945;
UB[2][2] = 0.006321;
UBinv = mat_inv(UB);
lambda = 1.179;
/* test PSD code */
psd.xScale = .7;
psd.yScale = 1.4;
psd.xZero = 128;
psd.yZero = 64;
psd.distance = 450.;
psd.gamma = 43.;
psd.nu = 12.3;
det2pol(&psd, 200, 111, &gamma, &nu);
pol2det(&psd, gamma, nu, &x, &y);
if (x == 200 && y == 111) {
printf("PSD test: OK\n");
}
/* test bisecting calculations. The tests were taken from a system known
to work
*/
status = 1;
sttex = 30.58704;
omex = 15.293520;
chiex = 129.053604;
phiex = 134.191132;
hkl = mat_creat(3, 1, ZERO_MATRIX);
hkl[0][0] = -2.;
hkl[1][0] = -2.;
hkl[2][0] = 4.;
zz1 = mat_mul(UB, hkl);
z1mToBisecting(lambda, zz1, &stt, &om, &chi, &phi);
if (ABS(stt - sttex) > .01 ||
ABS(om - omex) > .01 ||
ABS(chi - chiex) > .01 || ABS(phi - phiex) > .01) {
printf("Bisecting Angles calculation FAILED!!!!!!!!!\n");
printf("Expected: %12.4f %12.4f %12.4f %12.4f\n", sttex, omex, chiex,
phiex);
printf("Got: %12.4f %12.4f %12.4f %12.4f\n", stt, om, chi, phi);
status = 0;
}
z1FromAngles(lambda, sttex, omex, chiex, phiex, z1);
for (i = 0; i < 3; i++) {
if (ABS(zz1[i][0] - z1[i]) > .0001) {
printf("Calculation of Z1 from Angles FAILED!!!!\n");
printf("Expected: %8.4f %8.4f %8.4f\n", zz1[0][0], zz1[1][0],
zz1[2][0]);
printf("GOT: %8.4f %8.4f %8.4f\n", z1[0], z1[1], z1[2]);
status = 0;
}
}
hklback = mat_mul(UBinv, zz1);
for (i = 0; i < 3; i++) {
if (ABS(hklback[i][0] - hkl[i][0]) > .0001) {
printf("Back Calculation of HKL FAILED!!!!!!!!!\n");
printf("Expected: %8.4f %8.4f %8.4f\n", hkl[0][0], hkl[1][0],
hkl[2][0]);
printf("GOT: %8.4f %8.4f %8.4f\n", hklback[0][0],
hklback[1][0], hklback[2][0]);
status = 0;
}
}
if (status == 1) {
printf("Bisecting test: OK\n");
}
/* try turning in psi */
rotatePsi(omex, chiex, phiex, 30., &psiom, &psichi, &psiphi);
omex = 37.374298;
chiex = 123.068192;
phiex = 170.8209099;
if (ABS(psiom - omex) > .01 ||
ABS(psichi - chiex) > .01 || ABS(psiphi - phiex) > .01) {
printf("Psi Rotation test FAILED!!!!!!!!!\n");
printf("Expected: %12.4f %12.4f %12.4f\n", omex, chiex, phiex);
printf("Got: %12.4f %12.4f %12.4f\n", psiom, psichi, psiphi);
status = 0;
} else {
printf("Psi rotation test: OK\n");
}
/*
try to calculate a psi for a given omega offset
*/
omex = 15.293520;
chiex = 129.053604;
phiex = 134.191132;
if (psiForOmegaOffset(zz1, -5., chiex, phiex, &psi) != 1) {
printf("Failed to calculate PSI for a Omega Offset\n");
} else {
rotatePsi(omex, chiex, phiex, psi, &psiom, &psichi, &psiphi);
if (ABS(psiom + 5. - omex) > .5) {
printf("Failed to reach desired omega with calculated psi\n");
printf("Expected: %12.4f Got: %12.4f PSI = %12.4f\n", omex - 10.,
psiom, psi);
} else {
printf("Psi for Omega Offset: OK\n");
}
}
/*
try to calculate angles with a omega offset
*/
omex = 10.29;
chiex = 128.78;
phiex = 126.24;
sttex = 30.59;
z1ToAnglesWithOffset(lambda, zz1, -5., &stt, &om, &chi, &phi);
if (ABS(stt - sttex) > .01 ||
ABS(om - omex) > .01 ||
ABS(chi - chiex) > .01 || ABS(phi - phiex) > .01) {
printf("Angles calculation with Omega Offset FAILED!!!!!!!!!\n");
printf("Expected: %12.4f %12.4f %12.4f %12.4f\n", sttex, omex, chiex,
phiex);
printf("Got: %12.4f %12.4f %12.4f %12.4f\n", stt, om, chi, phi);
status = 0;
} else {
printf("Angles with Omega Offset: OK\n");
}
/*
test normal beam codes
*/
omex = 37.374298;
gaex = 19.322;
omnbex = -21.057;
nuex = 24.186;
z1mToBisecting(lambda, zz1, &stt, &om, &chi, &phi);
chi = 50.949;
phi = -45.8088;
if (!bisToNormalBeam(stt, om, chi, phi, &omeganb, &gamma, &nu)) {
printf("Error on normal beam calculation!!\n");
}
if (ABS(omeganb - omnbex) > .01 ||
ABS(gamma - gaex) > .01 || ABS(nu - nuex) > .01) {
printf("Normal Beam Calculation test FAILED!!!!!!!!!\n");
printf("Expected: %12.4f %12.4f %12.4f\n", omnbex, gaex, nuex);
printf("Got: %12.4f %12.4f %12.4f\n", omeganb, gamma, nu);
} else {
printf("Normal Beam Angle Calculation: OK\n");
}
z1FromNormalBeam(lambda, omnbex, gaex, nuex, z1);
status = 1;
for (i = 0; i < 3; i++) {
if (ABS(zz1[i][0] - z1[i]) > .0001) {
printf("Calculation of Z1 from Normal Beam Angles FAILED!!!!\n");
printf("Expected: %8.4f %8.4f %8.4f\n", zz1[0][0], zz1[1][0],
zz1[2][0]);
printf("GOT: %8.4f %8.4f %8.4f\n", z1[0], z1[1], z1[2]);
break;
status = 0;
}
}
if (status) {
printf("Normal Beam Backwards Calculation: OK\n");
}
/*
for interest: try to see how normal beam angles go if we rotate psi
*/
/*
sttex = 30.58704;
om = 15.293520;
chi = 50.949;
phi = -45.8088;
printf(" PSI Omega Gamma Nu\n");
for(i = 0; i < 36; i++){
psi = i*10.;
rotatePsi(om,chi,phi,psi,&psiom,&psichi,&psiphi);
status = bisToNormalBeam(sttex,psiom,psichi,psiphi,&omeganb, &gamma, &nu);
if(status){
printf("%12.4f%12.4f%12.4f%12.4f\n",psi,omeganb,gamma,nu);
}
}
*/
}
#endif
Reference in New Issue View Git Blame Copy Permalink