/* 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 #include #include #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