From 0011ff02a297441f60716c7ccc20615deae42dba Mon Sep 17 00:00:00 2001 From: Zack Singer Date: Fri, 13 Feb 2026 12:17:42 -0500 Subject: [PATCH 1/2] Apply transforms cleanup changes --- hexrd/core/extensions/__init__.py | 3 - hexrd/core/extensions/_transforms_CAPI.pyi | 174 -- ...ansforms_capi.pyi => transforms_c_api.pyi} | 0 hexrd/core/transforms/Makefile | 47 - hexrd/core/transforms/__init__.py | 28 - hexrd/core/transforms/debug_helpers.h | 61 - hexrd/core/transforms/new_capi/README.md | 9 - .../core/transforms/new_capi/angles_to_dvec.c | 124 -- .../core/transforms/new_capi/angles_to_gvec.c | 124 -- hexrd/core/transforms/new_capi/gvec_to_xy.c | 32 +- .../core/transforms/new_capi/make_beam_rmat.c | 27 +- .../transforms/new_capi/make_binary_rmat.c | 65 - .../transforms/new_capi/make_detector_rmat.c | 25 +- .../transforms/new_capi/make_rmat_of_expmap.c | 79 - .../transforms/new_capi/make_sample_rmat.c | 30 +- hexrd/core/transforms/new_capi/module.c | 90 +- .../core/transforms/new_capi/ndargs_helper.c | 3 +- .../core/transforms/new_capi/ndargs_helper.h | 2 +- hexrd/core/transforms/new_capi/new_func.c | 22 - .../new_capi/oscill_angles_of_HKLs.c | 26 +- .../core/transforms/new_capi/quat_distance.c | 31 +- hexrd/core/transforms/new_capi/reference.py | 57 - .../new_capi/rotate_vecs_about_axis.c | 27 +- .../new_capi/transforms_prototypes.h | 36 +- .../transforms/new_capi/unit_row_vector.c | 33 +- .../new_capi/validate_angle_ranges.c | 26 +- hexrd/core/transforms/new_capi/xf_new_capi.py | 691 -------- hexrd/core/transforms/new_capi/xy_to_gvec.c | 26 +- hexrd/core/transforms/old_xfcapi.py | 548 ------ hexrd/core/transforms/stdbool.h | 6 - hexrd/core/transforms/transforms_CAPI.c | 1221 ------------- hexrd/core/transforms/transforms_CAPI.h | 79 - hexrd/core/transforms/transforms_CFUNC.c | 1521 ----------------- hexrd/core/transforms/transforms_CFUNC.h | 95 - hexrd/core/transforms/xf.py | 1191 ------------- hexrd/core/transforms/xfcapi.py | 734 +++++++- hexrd/powder/wppf/texture.py | 4 +- setup.py | 17 +- tests/transforms/common.py | 2 +- .../test_angles_to_dvec_from_file.py | 2 +- .../test_angles_to_gvec_from_file.py | 2 +- tests/transforms/test_gvec_to_xy.py | 2 +- tests/transforms/test_gvec_to_xy_from_file.py | 2 +- .../test_make_beam_rmat_from_file.py | 2 +- tests/transforms/test_make_binary_rmat.py | 2 +- .../test_make_detector_rmat_from_file.py | 2 +- .../test_make_rmat_of_expmap_from_file.py | 2 +- .../test_make_sample_rmat_from_file.py | 2 +- .../test_quat_distance_from_file.py | 2 +- .../transforms/test_rotate_vecs_about_axis.py | 2 +- .../test_validate_angle_ranges_from_file.py | 2 +- tests/transforms/test_xy_to_gvec.py | 2 +- 52 files changed, 773 insertions(+), 6569 deletions(-) delete mode 100644 hexrd/core/extensions/__init__.py delete mode 100644 hexrd/core/extensions/_transforms_CAPI.pyi rename hexrd/core/extensions/{_new_transforms_capi.pyi => transforms_c_api.pyi} (100%) delete mode 100644 hexrd/core/transforms/Makefile delete mode 100644 hexrd/core/transforms/__init__.py delete mode 100644 hexrd/core/transforms/debug_helpers.h delete mode 100644 hexrd/core/transforms/new_capi/README.md delete mode 100644 hexrd/core/transforms/new_capi/angles_to_dvec.c delete mode 100644 hexrd/core/transforms/new_capi/angles_to_gvec.c delete mode 100644 hexrd/core/transforms/new_capi/make_binary_rmat.c delete mode 100644 hexrd/core/transforms/new_capi/make_rmat_of_expmap.c delete mode 100644 hexrd/core/transforms/new_capi/new_func.c delete mode 100644 hexrd/core/transforms/new_capi/reference.py delete mode 100644 hexrd/core/transforms/new_capi/xf_new_capi.py delete mode 100644 hexrd/core/transforms/old_xfcapi.py delete mode 100644 hexrd/core/transforms/stdbool.h delete mode 100644 hexrd/core/transforms/transforms_CAPI.c delete mode 100644 hexrd/core/transforms/transforms_CAPI.h delete mode 100644 hexrd/core/transforms/transforms_CFUNC.c delete mode 100644 hexrd/core/transforms/transforms_CFUNC.h delete mode 100644 hexrd/core/transforms/xf.py diff --git a/hexrd/core/extensions/__init__.py b/hexrd/core/extensions/__init__.py deleted file mode 100644 index 424934712..000000000 --- a/hexrd/core/extensions/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from . import _new_transforms_capi -from . import _transforms_CAPI -from . import inverse_distortion diff --git a/hexrd/core/extensions/_transforms_CAPI.pyi b/hexrd/core/extensions/_transforms_CAPI.pyi deleted file mode 100644 index e45b3c3d2..000000000 --- a/hexrd/core/extensions/_transforms_CAPI.pyi +++ /dev/null @@ -1,174 +0,0 @@ -from __future__ import annotations - -import numpy as np -from numpy.typing import NDArray - -def anglesToGVec( - angs: NDArray[np.float64], # (n, 3) - bHat_l: NDArray[np.float64], # (3,) - eHat_l: NDArray[np.float64], # (3,) - chi: float, - rMat_c: NDArray[np.float64], # (3, 3) -) -> NDArray[np.float64]: # (n, 3) - ... - -def anglesToDVec( - angs: NDArray[np.float64], # (n, 3) - bHat_l: NDArray[np.float64], # (3,) - eHat_l: NDArray[np.float64], # (3,) - chi: float, - rMat_c: NDArray[np.float64], # (3, 3) -) -> NDArray[np.float64]: # (n, 3) - ... - -def gvecToDetectorXY( - gVec_c: NDArray[np.float64], # (n, 3) - rMat_d: NDArray[np.float64], # (3, 3) - rMat_s: NDArray[np.float64], # (3, 3) - rMat_c: NDArray[np.float64], # (3, 3) - tVec_d: NDArray[np.float64], # (3,) - tVec_s: NDArray[np.float64], # (3,) - tVec_c: NDArray[np.float64], # (3,) - beamVec: NDArray[np.float64], # (3,) -) -> NDArray[np.float64]: # (n, 2) - ... - -def gvecToDetectorXYArray( - gVec_c: NDArray[np.float64], # (n, 3) - rMat_d: NDArray[np.float64], # (3, 3) - rMat_s: NDArray[np.float64], # (n, 3, 3) - rMat_c: NDArray[np.float64], # (3, 3) - tVec_d: NDArray[np.float64], # (3,) - tVec_s: NDArray[np.float64], # (3,) - tVec_c: NDArray[np.float64], # (3,) - beamVec: NDArray[np.float64], # (3,) -) -> NDArray[np.float64]: # (n, 2) - ... - -def detectorXYToGvec( - xy_det: NDArray[np.float64], # (n, 2) - rMat_d: NDArray[np.float64], # (3, 3) - rMat_s: NDArray[np.float64], # (3, 3) - tVec_d: NDArray[np.float64], # (3,) - tVec_s: NDArray[np.float64], # (3,) - tVec_c: NDArray[np.float64], # (3,) - beamVec: NDArray[np.float64], # (3,) - etaVec: NDArray[np.float64], # (3,) -) -> tuple[tuple[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]]: - """ - Returns: - ((tTh, eta), gVec_l) - where: - tTh: (n,) float64 - eta: (n,) float64 - gVec_l: (n, 3) float64 - """ - ... - -def detectorXYToGvecArray( - xy_det: NDArray[np.float64], # (n, 2) - rMat_d: NDArray[np.float64], # (3, 3) - rMat_s: NDArray[np.float64], # (n, 3, 3) - tVec_d: NDArray[np.float64], # (3,) - tVec_s: NDArray[np.float64], # (3,) - tVec_c: NDArray[np.float64], # (3,) - beamVec: NDArray[np.float64], # (3,) - etaVec: NDArray[np.float64], # (3,) -) -> tuple[tuple[NDArray[np.float64], NDArray[np.float64]], NDArray[np.float64]]: - """ - Returns: - ((tTh, eta), gVec_l) - where: - tTh: (n,) float64 - eta: (n,) float64 - gVec_l: (n, 3) float64 - """ - ... - -def oscillAnglesOfHKLs( - hkls: NDArray[np.float64], # (n, 3) - chi: float, - rMat_c: NDArray[np.float64], # (3, 3) - bMat: NDArray[np.float64], # (3, 3) - wavelength: float, - vInv_s: NDArray[np.float64], # (6,) - beamVec: NDArray[np.float64], # (3,) - etaVec: NDArray[np.float64], # (3,) -) -> tuple[NDArray[np.float64], NDArray[np.float64]]: - """ - Returns: - (oangs0, oangs1) each (n, 3) float64 - """ - ... - -def unitRowVector( - vecIn: NDArray[np.float64], # (n,) -) -> NDArray[np.float64]: # (n,) - ... - -def unitRowVectors( - vecIn: NDArray[np.float64], # (m, n) -) -> NDArray[np.float64]: # (m, n) - ... - -def makeDetectorRotMat( - tiltAngles: NDArray[np.float64], # (3,) -) -> NDArray[np.float64]: # (3, 3) - ... - -def makeOscillRotMat( - oscillAngles: NDArray[np.float64], # (2,) -) -> NDArray[np.float64]: # (3, 3) - ... - -def makeOscillRotMatArray( - chi: float, - omeArray: NDArray[np.float64], # (n,) -) -> NDArray[np.float64]: # (n, 3, 3) - ... - -def makeRotMatOfExpMap( - expMap: NDArray[np.float64], # (3,) -) -> NDArray[np.float64]: # (3, 3) - ... - -def makeRotMatOfQuat( - quat: NDArray[np.float64], # (n, 4) -) -> NDArray[np.float64]: # (n, 3, 3) - ... - -def makeBinaryRotMat( - axis: NDArray[np.float64], # (3,) -) -> NDArray[np.float64]: # (3, 3) - ... - -def makeEtaFrameRotMat( - bHat: NDArray[np.float64], # (3,) - eHat: NDArray[np.float64], # (3,) -) -> NDArray[np.float64]: # (3, 3) - ... - -def validateAngleRanges( - angList: NDArray[np.float64], # (na,) - angMin: NDArray[np.float64], # (nmin,) - angMax: NDArray[np.float64], # (nmin,) - ccw: object, # expects bool-ish -) -> NDArray[np.bool_]: # (na,) - ... - -def rotate_vecs_about_axis( - angles: NDArray[np.float64], # (1,) or (n,) - axes: NDArray[np.float64], # (1,3) or (n,3) - vecs: NDArray[np.float64], # (n,3) -) -> NDArray[np.float64]: # (n,3) - ... - -def quat_distance( - q1: NDArray[np.float64], # (4,) - q2: NDArray[np.float64], # (4,) - qsym: NDArray[np.float64], # (4, nsym) -) -> float: ... -def homochoricOfQuat( - quat: NDArray[np.float64], # (n, 4) -) -> NDArray[np.float64]: # (n, 3) - ... diff --git a/hexrd/core/extensions/_new_transforms_capi.pyi b/hexrd/core/extensions/transforms_c_api.pyi similarity index 100% rename from hexrd/core/extensions/_new_transforms_capi.pyi rename to hexrd/core/extensions/transforms_c_api.pyi diff --git a/hexrd/core/transforms/Makefile b/hexrd/core/transforms/Makefile deleted file mode 100644 index 54d99d813..000000000 --- a/hexrd/core/transforms/Makefile +++ /dev/null @@ -1,47 +0,0 @@ -#============================================================================ -# sources - -SRCS = transforms_CAPI.c transforms_CFUNC.c -OBJS = $(SRCS:.c=.o) - -#============================================================================ -# library name - -LIBBASENAME = _transforms_CAPI - -LIBNAME = $(LIBBASENAME).so - -#============================================================================ -# compiler options - -CC = gcc -CFLAGS = -O3 -fPIC -CDEFINES = - -#============================================================================ -# header location - -PYTHON_INCLUDE_DIR = ${HOME}/opt/include/python2.7 -NUMPY_INCLUDE_DIR = ${HOME}/opt/lib/python2.7/site-packages/numpy/core/include/numpy - -INCPATH = -I$(PYTHON_INCLUDE_DIR) -I$(NUMPY_INCLUDE_DIR) - -#============================================================================ -# targets - -default: lib - -lib: $(LIBNAME) - -$(LIBBASENAME).so: $(OBJS) - $(CC) -shared $(CFLAGS) -flat_namespace -o $(LIBBASENAME).so $(OBJS) -lm - -clean: - $(RM) -rf *.o $(LIBBASENAME).so - -#============================================================================ -# suffix rules - -.SUFFIXES: .c -.c.o: - $(CC) $(CFLAGS) $(CDEFINES) $(INCPATH) -c $< diff --git a/hexrd/core/transforms/__init__.py b/hexrd/core/transforms/__init__.py deleted file mode 100644 index 12ff60c85..000000000 --- a/hexrd/core/transforms/__init__.py +++ /dev/null @@ -1,28 +0,0 @@ -# ============================================================ -# Copyright (c) 2012, Lawrence Livermore National Security, LLC. -# Produced at the Lawrence Livermore National Laboratory. -# Written by Joel Bernier and others. -# LLNL-CODE-529294. -# All rights reserved. -# -# This file is part of HEXRD. For details on dowloading the source, -# see the file COPYING. -# -# Please also see the file LICENSE. -# -# This program is free software; you can redistribute it and/or modify it under the -# terms of the GNU Lesser General Public License (as published by the Free Software -# Foundation) version 2.1 dated February 1999. -# -# This program is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the IMPLIED WARRANTY OF MERCHANTABILITY -# or FITNESS FOR A PARTICULAR PURPOSE. See the terms and conditions of the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU Lesser General Public -# License along with this program (see file LICENSE); if not, write to -# the Free Software Foundation, Inc., 59 Temple Place, Suite 330, -# Boston, MA 02111-1307 USA or visit . -# ============================================================ -"""Tools or X-ray diffraction analysis""" -from . import xfcapi diff --git a/hexrd/core/transforms/debug_helpers.h b/hexrd/core/transforms/debug_helpers.h deleted file mode 100644 index aacce7d5a..000000000 --- a/hexrd/core/transforms/debug_helpers.h +++ /dev/null @@ -1,61 +0,0 @@ -#include - - - -static void debug_dump_val(const char *name, double val) -{ - printf("%s: %10.6f\n", name, val); -} - -static void debug_dump_m33(const char *name, double *array) -{ - printf("%s:\n", name); - printf("\t%10.6f %10.6f %10.6f\n", array[0], array[1], array[2]); - printf("\t%10.6f %10.6f %10.6f\n", array[3], array[4], array[5]); - printf("\t%10.6f %10.6f %10.6f\n", array[6], array[7], array[8]); -} - -static void debug_dump_v3(const char *name, double *vec) -{ - printf("%s: %10.6f %10.6f %10.6f\n", name, vec[0], vec[1], vec[2]); -} - -/* ============================================================================ - * These can be used to initialize and check buffers if we suspect the - * code may leave it uninitialized - * ============================================================================ - */ - -#define SNAN_HI 0x7ff700a0 -#define SNAN_LO 0xbad0feed -void fill_signaling_nans(double *arr, int count) -{ - int i; - npy_uint32 *arr_32 = (npy_uint32 *)arr; - /* Fills an array with signaling nans to detect errors - * Use the 0x7ff700a0bad0feed as the pattern - */ - for (i = 0; i < count; ++i) - { - arr_32[2*i+0] = SNAN_LO; - arr_32[2*i+1] = SNAN_HI; - } -} - -int detect_signaling_nans(double *arr, int count) -{ - int i; - npy_uint32 *arr_32 = (npy_uint32 *) arr; - for (i = 0; i < count; ++i) - { - if (arr_32[2*i+0] == SNAN_LO && - arr_32[2*i+1] == SNAN_HI) - { - return 1; - } - } - - return 0; -} -#undef SNAN_HI -#undef SNAN_LO diff --git a/hexrd/core/transforms/new_capi/README.md b/hexrd/core/transforms/new_capi/README.md deleted file mode 100644 index c1641d18f..000000000 --- a/hexrd/core/transforms/new_capi/README.md +++ /dev/null @@ -1,9 +0,0 @@ -These files were copied from Oscar Villellas' xrd-transforms library -located here: https://github.com/ovillellas/xrd-transforms - -At commit: b94f8b2d7839d883829d00a2adc5bec9c80e0116 - -They have been modified a little since then (see commit history). - -Our goal is to eventually migrate all of our xfcapi function calls to these, -but we need to test them thoroughly, and do them one at a time. diff --git a/hexrd/core/transforms/new_capi/angles_to_dvec.c b/hexrd/core/transforms/new_capi/angles_to_dvec.c deleted file mode 100644 index f1618c389..000000000 --- a/hexrd/core/transforms/new_capi/angles_to_dvec.c +++ /dev/null @@ -1,124 +0,0 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -# include "ndargs_helper.h" -#endif - -XRD_CFUNCTION void -angles_to_dvec(size_t nvecs, double * angs, - double * bHat_l, double * eHat_l, - double chi, double * rMat_c, - double * gVec_c) -{ - /* - * takes an angle spec (2*theta, eta, omega) for nvecs g-vectors and - * returns the unit d-vector components in the crystal frame - * - * For unit d-vector in the lab frame, spec rMat_c = Identity and - * overwrite the omega values with zeros - */ - size_t i, j, k, l; - double rMat_e[9], rMat_s[9], rMat_ctst[9]; - double gVec_e[3], gVec_l[3], gVec_c_tmp[3]; - - /* Need eta frame cob matrix (could omit for standard setting) */ - make_beam_rmat(bHat_l, eHat_l, rMat_e); - - /* make vector array */ - for (i=0; i -# include -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_anglesToDVec = - "c module implementation of angles_to_dvec.\n" - "Please use the Python wrapper.\n"; - -XRD_PYTHON_WRAPPER PyObject * -python_anglesToDVec(PyObject * self, PyObject * args) -{ - /* - API interface in Python is: - angs, beam_vec=None, eta_vec=None, chi=None, rmat_c=None - - right now, defaults are handled by the wrappers. This C module function - expects angs to be already 2d, and all optional arguments to have the - defaults filled by the calling wrapper. - - TODO: handle defaults here? - */ - nah_array angs = { NULL, "angs", NAH_TYPE_DP_FP, { 3, NAH_DIM_ANY }}; - nah_array beam_vec = { NULL, "beam_vec", NAH_TYPE_DP_FP, { 3 }}; - nah_array eta_vec = { NULL, "eta_vec", NAH_TYPE_DP_FP, { 3 }}; - nah_array rmat_c = { NULL, "rmat_c", NAH_TYPE_DP_FP, { 3, 3 }}; - PyArrayObject *result = NULL; - double chi; - - if ( !PyArg_ParseTuple(args,"O&O&O&dO&", - nah_array_converter, &angs, - nah_array_converter, &beam_vec, - nah_array_converter, &eta_vec, - &chi, - nah_array_converter, &rmat_c)) - return NULL; - - /* Allocate C-style array for return data */ - result = (PyArrayObject*)PyArray_EMPTY(PyArray_NDIM(angs.pyarray), - PyArray_SHAPE(angs.pyarray), - NPY_DOUBLE, - 0); - - /* Call the actual function */ - angles_to_dvec(PyArray_DIM(angs.pyarray, 0), - (double *)PyArray_DATA(angs.pyarray), - (double *)PyArray_DATA(beam_vec.pyarray), - (double *)PyArray_DATA(eta_vec.pyarray), - chi, - (double *)PyArray_DATA(rmat_c.pyarray), - (double *)PyArray_DATA(result)); - - /* Build and return the nested data structure */ - return((PyObject*)result); -} - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/angles_to_gvec.c b/hexrd/core/transforms/new_capi/angles_to_gvec.c deleted file mode 100644 index cdd2edec8..000000000 --- a/hexrd/core/transforms/new_capi/angles_to_gvec.c +++ /dev/null @@ -1,124 +0,0 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -# include "ndargs_helper.h" -#endif - - -XRD_CFUNCTION void -angles_to_gvec(size_t nvecs, double * angs, - double * bHat_l, double * eHat_l, - double chi, double * rMat_c, - double * gVec_c) -{ - /* - * takes an angle spec (2*theta, eta, omega) for nvecs g-vectors and - * returns the unit g-vector components in the crystal frame - * - * For unit g-vector in the lab frame, spec rMat_c = Identity and - * overwrite the omega values with zeros - */ - size_t i, j, k, l; - double rMat_e[9], rMat_s[9], rMat_ctst[9]; - double gVec_e[3], gVec_l[3], gVec_c_tmp[3]; - - /* Need eta frame cob matrix (could omit for standard setting) */ - make_beam_rmat(bHat_l, eHat_l, rMat_e); - - /* make vector array */ - for (i=0; i -# include -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_anglesToGVec = - "c module implementation of angles_to_gvec.\n" - "Please use the Python wrapper.\n"; - -XRD_PYTHON_WRAPPER PyObject * -python_anglesToGVec(PyObject * self, PyObject * args) -{ - /* - API interface in Python is: - angs, beam_vec=None, eta_vec=None, chi=None, rmat_c=None - - currently, defaults are handled by the wrapper, so the C-module - is always called with all arguments. Wrapper always guarantees - that the input array is 2d. - TODO: handle defaults here? - */ - nah_array angs = { NULL, "angs", NAH_TYPE_DP_FP, { 3, NAH_DIM_ANY }}; - nah_array beam_vec = { NULL, "beam_vec", NAH_TYPE_DP_FP, { 3 }}; - nah_array eta_vec = { NULL, "eta_vec", NAH_TYPE_DP_FP, { 3 }}; - nah_array rmat_c = { NULL, "rmat_c", NAH_TYPE_DP_FP, { 3, 3 }}; - PyArrayObject *result = NULL; - double chi; - - if ( !PyArg_ParseTuple(args,"O&O&O&dO&", - nah_array_converter, &angs, - nah_array_converter, &beam_vec, - nah_array_converter, &eta_vec, - &chi, - nah_array_converter, &rmat_c)) - return NULL; - - /* Allocate C-style array for return data */ - result = (PyArrayObject*)PyArray_EMPTY(PyArray_NDIM(angs.pyarray), - PyArray_SHAPE(angs.pyarray), - NPY_DOUBLE, - 0); - - /* Call the actual function */ - angles_to_gvec(PyArray_DIM(angs.pyarray, 0), - (double *)PyArray_DATA(angs.pyarray), - (double *)PyArray_DATA(beam_vec.pyarray), - (double *)PyArray_DATA(eta_vec.pyarray), - chi, - (double *)PyArray_DATA(rmat_c.pyarray), - (double *)PyArray_DATA(result)); - - /* Build and return the nested data structure */ - return((PyObject*)result); -} - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/gvec_to_xy.c b/hexrd/core/transforms/new_capi/gvec_to_xy.c index 3c1e9d7af..7662798d7 100644 --- a/hexrd/core/transforms/new_capi/gvec_to_xy.c +++ b/hexrd/core/transforms/new_capi/gvec_to_xy.c @@ -1,10 +1,4 @@ -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -# include "ndargs_helper.h" -#endif - #define GV2XY_GROUP_SIZE 128 #define GVEC_TO_XY_FUNC gvec_to_xy_vect @@ -326,7 +320,7 @@ ray_plane_intersect(const double *origin, const double *vect, } -XRD_CFUNCTION void +static void gvec_to_xy(size_t npts, const double *gVec_cs, const double *rMat_d, const double *rMat_ss, const double *rMat_c, const double *tVec_d, const double *tVec_s, const double *tVec_c, @@ -459,7 +453,7 @@ typedef void (*rays_to_detector_func) (const void *detector_data, double * restrict projected_xy, size_t projected_xy_count); /* experimental: */ -XRD_CFUNCTION void +static void gvec_to_xy_detector(size_t npts, const double *gVec_cs, const double *rMat_ss, const double *rMat_c, const double *tVec_s, const double *tVec_c, @@ -571,7 +565,7 @@ gvec_to_xy_detector(size_t npts, const double *gVec_cs, } } -XRD_CFUNCTION void +static void gvec_to_xy_vect(size_t npts, const double *gVec_cs, const double *rMat_d, const double *rMat_ss, const double *rMat_c, const double *tVec_d, const double *tVec_s, const double *tVec_c, @@ -592,20 +586,11 @@ gvec_to_xy_vect(size_t npts, const double *gVec_cs, result, flags, rays_to_planar_detector, &planar_detector_data); } -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_gvecToDetectorXY = +static const char *docstring_gvecToDetectorXY = "c module implementation of gvec_to_xy (single sample).\n" "Please use the Python wrapper.\n"; -XRD_PYTHON_WRAPPER const char *docstring_gvecToDetectorXYArray = +static const char *docstring_gvecToDetectorXYArray = "c module implementation of gvec_to_xy (multi sample).\n" "Please use the Python wrapper.\n"; @@ -671,7 +656,7 @@ normalize_beam(double *in, double *work) (m, 2) ndarray containing the intersections of m <= n diffracted beams associated with gVecs */ -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_gvecToDetectorXY(PyObject * self, PyObject * args) { nah_array gVec_c = { NULL, "gvec_c", NAH_TYPE_DP_FP, { 3, NAH_DIM_ANY }}; @@ -756,7 +741,7 @@ python_gvecToDetectorXY(PyObject * self, PyObject * args) (m, 2) ndarray containing the intersections of m <= n diffracted beams associated with gVecs */ -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_gvecToDetectorXYArray(PyObject * self, PyObject * args) { nah_array gVec_c = { NULL, "gVec_c", NAH_TYPE_DP_FP, { 3, NAH_DIM_ANY }}; @@ -828,6 +813,3 @@ python_gvecToDetectorXYArray(PyObject * self, PyObject * args) return PyErr_NoMemory(); } - - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/make_beam_rmat.c b/hexrd/core/transforms/new_capi/make_beam_rmat.c index 1ea02870c..edd70d23a 100644 --- a/hexrd/core/transforms/new_capi/make_beam_rmat.c +++ b/hexrd/core/transforms/new_capi/make_beam_rmat.c @@ -1,13 +1,4 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -# include "ndargs_helper.h" -#endif - - -XRD_CFUNCTION int -make_beam_rmat(double * bPtr, double * ePtr, double * rPtr) +static int make_beam_rmat(double * bPtr, double * ePtr, double * rPtr) { /* * This function generates a COB matrix that takes components in the @@ -57,20 +48,11 @@ make_beam_rmat(double * bPtr, double * ePtr, double * rPtr) } -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_makeEtaFrameRotMat = +static const char *docstring_makeEtaFrameRotMat = "c module implementation of make_beam_mat.\n" "Please use the Python wrapper.\n"; -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_makeEtaFrameRotMat(PyObject * self, PyObject * args) { PyArrayObject *rMat=NULL; @@ -127,6 +109,3 @@ python_makeEtaFrameRotMat(PyObject * self, PyObject * args) done: return (PyObject *)rMat; } - - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/make_binary_rmat.c b/hexrd/core/transforms/new_capi/make_binary_rmat.c deleted file mode 100644 index cc50430fb..000000000 --- a/hexrd/core/transforms/new_capi/make_binary_rmat.c +++ /dev/null @@ -1,65 +0,0 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -#endif - - -XRD_CFUNCTION void -make_binary_rmat(const double *aPtr, double * restrict rPtr) -{ - int i, j; - - for (i=0; i<3; i++) { - for (j=0; j<3; j++) { - rPtr[3*i+j] = 2.0*aPtr[i]*aPtr[j]; - } - rPtr[3*i+i] -= 1.0; - } -} - - -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# include "ndargs_helper.h" -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_makeBinaryRotMat = - "c module implementation of make_binary_rmat.\n" - "Please use the Python wrapper.\n"; - -XRD_PYTHON_WRAPPER PyObject * -python_makeBinaryRotMat(PyObject * self, PyObject * args) -{ - nah_array axis = { NULL, "axis", NAH_TYPE_DP_FP, { 3 }}; - PyArrayObject *result = NULL; - npy_intp dims[2] = { 3, 3 }; - - /* Parse arguments */ - if (!PyArg_ParseTuple(args, "O&", - nah_array_converter, &axis)) - return NULL; - - /* Allocate the result matrix with appropriate dimensions and type */ - result = (PyArrayObject*)PyArray_EMPTY(2, dims, NPY_DOUBLE, 0); - if (NULL == result) - goto fail_alloc; - - /* Call the actual function */ - make_binary_rmat((double *)PyArray_DATA(axis.pyarray), - (double *)PyArray_DATA(result)); - - return (PyObject*)result; - - fail_alloc: - Py_XDECREF(result); - - return PyErr_NoMemory(); -} - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/make_detector_rmat.c b/hexrd/core/transforms/new_capi/make_detector_rmat.c index 9be3bf096..2712005a6 100644 --- a/hexrd/core/transforms/new_capi/make_detector_rmat.c +++ b/hexrd/core/transforms/new_capi/make_detector_rmat.c @@ -1,12 +1,4 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -#endif - - -XRD_CFUNCTION void -make_detector_rmat(double * tPtr, double * rPtr) +static void make_detector_rmat(double * tPtr, double * rPtr) { int i; double c[3],s[3]; @@ -28,20 +20,11 @@ make_detector_rmat(double * tPtr, double * rPtr) } -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_makeDetectorRotMat = +static const char *docstring_makeDetectorRotMat = "c module implementation of makeDetectorRotMat.\n" "Please use the Python wrapper.\n"; -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_makeDetectorRotMat(PyObject * self, PyObject * args) { PyArrayObject *tiltAngles, *rMat; @@ -74,5 +57,3 @@ python_makeDetectorRotMat(PyObject * self, PyObject * args) return((PyObject*)rMat); } - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/make_rmat_of_expmap.c b/hexrd/core/transforms/new_capi/make_rmat_of_expmap.c deleted file mode 100644 index e21209d27..000000000 --- a/hexrd/core/transforms/new_capi/make_rmat_of_expmap.c +++ /dev/null @@ -1,79 +0,0 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -# include "ndargs_helper.h" -#endif - - -XRD_CFUNCTION void -make_rmat_of_expmap(double *ePtr, double *rPtr) -{ - double e0 = ePtr[0]; - double e1 = ePtr[1]; - double e2 = ePtr[2]; - double sqr_phi = e0*e0 + e1*e1 + e2*e2; - - if ( sqr_phi > epsf ) { - double phi = sqrt(sqr_phi); - double s = sin(phi)/phi; - double c = (1.0-cos(phi))/sqr_phi; - - rPtr[0] = 1.0 - c*(e1*e1 + e2*e2); - rPtr[1] = -s*e2 + c*e0*e1; - rPtr[2] = s*e1 + c*e0*e2; - rPtr[3] = s*e2 + c*e1*e0; - rPtr[4] = 1.0 - c*(e2*e2 + e0*e0); - rPtr[5] = -s*e0 + c*e1*e2; - rPtr[6] = -s*e1 + c*e2*e0; - rPtr[7] = s*e0 + c*e2*e1; - rPtr[8] = 1.0 - c*(e0*e0 + e1*e1); - } else { - m33_set_identity(rPtr); - } -} - - -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_makeRotMatOfExpMap = - "c module implementation of make_rmat_of_expmap.\n" - "Please use the Python wrapper.\n"; - -XRD_PYTHON_WRAPPER PyObject * -python_makeRotMatOfExpMap(PyObject * self, PyObject * args) -{ - nah_array exp_map = { NULL, "exp_map", NAH_TYPE_DP_FP, { 3 }}; - PyArrayObject *result = NULL; - npy_intp dims[2] = { 3, 3 }; - - /* Parse arguments */ - if (!PyArg_ParseTuple(args, "O&", - nah_array_converter, &exp_map)) - return NULL; - - /* Allocate the result matrix with appropriate dimensions and type */ - result = (PyArrayObject*)PyArray_EMPTY(2, dims, NPY_DOUBLE, 0); - if (NULL == result) - goto fail_alloc; - - /* Call the actual function */ - make_rmat_of_expmap(PyArray_DATA(exp_map.pyarray), - PyArray_DATA(result)); - - return (PyObject *)result; - - fail_alloc: - Py_XDECREF(result); - - return PyErr_NoMemory(); -} - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/make_sample_rmat.c b/hexrd/core/transforms/new_capi/make_sample_rmat.c index 3c7c2f501..a0989bce1 100644 --- a/hexrd/core/transforms/new_capi/make_sample_rmat.c +++ b/hexrd/core/transforms/new_capi/make_sample_rmat.c @@ -1,12 +1,4 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -# include "ndargs_helper.h" -#endif - -XRD_CFUNCTION void -make_sample_rmat(double chi, double ome, double *result_rmat) +static void make_sample_rmat(double chi, double ome, double *result_rmat) { double cchi, schi, come, some; double * restrict rmat = result_rmat; @@ -28,7 +20,7 @@ make_sample_rmat(double chi, double ome, double *result_rmat) } -XRD_CFUNCTION void +static void make_sample_rmat_array(double chi, const double *ome_ptr, size_t ome_count, double *result_ptr) { double cchi,schi; @@ -57,23 +49,11 @@ make_sample_rmat_array(double chi, const double *ome_ptr, size_t ome_count, doub } } - - -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_makeOscillRotMat = +static const char *docstring_makeOscillRotMat = "c module implementation of make_sample_rotmat.\n" "Please use the Python wrapper.\n"; - -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_makeOscillRotMat(PyObject * self, PyObject * args) { nah_array ome = { NULL, "ome", NAH_TYPE_DP_FP, { NAH_DIM_ANY }}; @@ -105,5 +85,3 @@ python_makeOscillRotMat(PyObject * self, PyObject * args) return PyErr_NoMemory(); } - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/module.c b/hexrd/core/transforms/new_capi/module.c index bbff03a30..4516fc0ea 100644 --- a/hexrd/core/transforms/new_capi/module.c +++ b/hexrd/core/transforms/new_capi/module.c @@ -15,7 +15,7 @@ * Note the supporting CONCAT and STR macros... * ============================================================================= */ -#define THIS_MODULE_NAME _new_transforms_capi +#define THIS_MODULE_NAME transforms_c_api #define _CONCAT(a,b) a ## b #define CONCAT(a,b) _CONCAT(a,b) @@ -50,27 +50,18 @@ * statics. * ============================================================================= */ -#define XRD_SINGLE_COMPILE_UNIT 1 -#define XRD_INCLUDE_PYTHON_WRAPPERS 1 -#define XRD_CFUNCTION static -#define XRD_PYTHON_WRAPPER static - #include "transforms_types.h" #include "transforms_utils.h" #include "transforms_prototypes.h" #include "ndargs_helper.h" #include "ndargs_helper.c" -#include "angles_to_gvec.c" -#include "angles_to_dvec.c" #include "gvec_to_xy.c" #include "xy_to_gvec.c" #include "oscill_angles_of_HKLs.c" #include "unit_row_vector.c" #include "make_detector_rmat.c" #include "make_sample_rmat.c" -#include "make_rmat_of_expmap.c" -#include "make_binary_rmat.c" #include "make_beam_rmat.c" #include "validate_angle_ranges.c" #include "rotate_vecs_about_axis.c" @@ -83,15 +74,10 @@ * ============================================================================= */ -/*#define EXPORT_METHOD(name) \ - { STR(name), CONCAT(python_, name), METH_VARARGS, CONCAT(docstring_, name) } -*/ #define EXPORT_METHOD(name) \ { STR(name), CONCAT(python_, name), METH_VARARGS, "" } static PyMethodDef _module_methods[] = { - EXPORT_METHOD(anglesToGVec), /* angles_to_gvec */ - EXPORT_METHOD(anglesToDVec), /* angles_to_dvec */ EXPORT_METHOD(gvecToDetectorXY), /* gvec_to_xy */ EXPORT_METHOD(gvecToDetectorXYArray), /* gvec_to_xy */ EXPORT_METHOD(detectorXYToGvec), /* xy_to_gvec */ @@ -99,57 +85,14 @@ static PyMethodDef _module_methods[] = { EXPORT_METHOD(unitRowVector), /* unit_vector */ EXPORT_METHOD(unitRowVectors), /* unit_vector */ EXPORT_METHOD(makeOscillRotMat), /* make_sample_rmat */ - EXPORT_METHOD(makeRotMatOfExpMap), /* make_rmat_of_expmap */ EXPORT_METHOD(makeDetectorRotMat), /* make_detector_rmat */ - EXPORT_METHOD(makeBinaryRotMat), /* make_binary_rmat */ EXPORT_METHOD(makeEtaFrameRotMat), /* make_beam_rmat */ EXPORT_METHOD(validateAngleRanges), /* validate_angle_ranges */ EXPORT_METHOD(rotate_vecs_about_axis), /* rotate_vecs_about_axis */ EXPORT_METHOD(quat_distance), /* quat_distance */ - /* adapted */ - /* {"anglesToGVec",anglesToGVec,METH_VARARGS,"take angle tuples to G-vectors"},*/ - /* {"anglesToDVec",anglesToDVec,METH_VARARGS,"take angle tuples to unit diffraction vectors"},*/ - /* {"gvecToDetectorXY",gvecToDetectorXY,METH_VARARGS,""},*/ - /* {"gvecToDetectorXYArray",gvecToDetectorXYArray,METH_VARARGS,""},*/ - /* {"detectorXYToGvec",detectorXYToGvec,METH_VARARGS,"take cartesian coordinates to G-vectors"},*/ - /* {"oscillAnglesOfHKLs",oscillAnglesOfHKLs,METH_VARARGS,"solve angle specs for G-vectors"},*/ - /* {"unitRowVector",unitRowVector,METH_VARARGS,"Normalize a single row vector"},*/ - /* {"unitRowVectors",unitRowVectors,METH_VARARGS,"Normalize a collection of row vectors"},*/ - /* {"makeOscillRotMat",makeOscillRotMat,METH_VARARGS,""},*/ - /* {"makeOscillRotMatArray",makeOscillRotMatArray,METH_VARARGS,""},*/ - /* {"makeRotMatOfExpMap",makeRotMatOfExpMap,METH_VARARGS,""},*/ - /* {"makeBinaryRotMat",makeBinaryRotMat,METH_VARARGS,""},*/ - /* {"makeEtaFrameRotMat",makeEtaFrameRotMat,METH_VARARGS,"Make eta basis COB matrix"},*/ - /* {"validateAngleRanges",validateAngleRanges,METH_VARARGS,""}, */ - /* {"rotate_vecs_about_axis",rotate_vecs_about_axis,METH_VARARGS,"Rotate vectors about an axis"},*/ - /* {"quat_distance",quat_distance,METH_VARARGS,"Compute distance between two unit quaternions"},*/ - - /* adapted but not part of the API, so not exported */ - /* {"makeDetectorRotMat",makeDetectorRotMat,METH_VARARGS,""},*/ - - /* removed... */ - /* {"makeGVector",makeGVector,METH_VARARGS,"Make G-vectors from hkls and B-matrix"},*/ - /* {"arccosSafe",arccosSafe,METH_VARARGS,""},*/ - /* {"angularDifference",angularDifference,METH_VARARGS,"difference for cyclical angles"},*/ - /* {"mapAngle",mapAngle,METH_VARARGS,"map cyclical angles to specified range"},*/ - /* {"columnNorm",columnNorm,METH_VARARGS,""},*/ - /* {"rowNorm",rowNorm,METH_VARARGS,""},*/ - /* {"makeRotMatOfQuat",makeRotMatOfQuat,METH_VARARGS,""},*/ - /* {"homochoricOfQuat",homochoricOfQuat,METH_VARARGS,"Compute homochoric parameterization of list of unit quaternions"},*/ {NULL,NULL,0,NULL} /* sentinel */ }; -/* - * In Python 3 the entry point for a C module changes slightly, but both can - * be supported with little effort with some conditionals and macro magic - */ - -#if PY_VERSION_HEX >= 0x03000000 -# define IS_PY3K -#endif - -#if defined(IS_PY3K) -/* a module definition structure is required in Python 3 */ static struct PyModuleDef moduledef = { PyModuleDef_HEAD_INIT, STR(THIS_MODULE_NAME), @@ -161,30 +104,15 @@ static struct PyModuleDef moduledef = { NULL, NULL }; -#endif - -#if defined(IS_PY3K) -# define MODULE_INIT_FUNC_NAME CONCAT(PyInit_, THIS_MODULE_NAME) -# define MODULE_RETURN(module) return module -#else -# define MODULE_INIT_FUNC_NAME CONCAT(init, THIS_MODULE_NAME) -# define MODULE_RETURN(module) return -#endif -PyMODINIT_FUNC -MODULE_INIT_FUNC_NAME(void) -{ - PyObject *module; -#if defined(IS_PY3K) - module = PyModule_Create(&moduledef); -#else - module = Py_InitModule(STR(THIS_MODULE_NAME),_module_methods); -#endif - if (NULL == module) - MODULE_RETURN(module); +PyMODINIT_FUNC CONCAT(PyInit_, THIS_MODULE_NAME)(void) +{ + PyObject *module = PyModule_Create(&moduledef); + if (module == NULL) { + return NULL; + } import_array(); - - MODULE_RETURN(module); -} + return module; +} \ No newline at end of file diff --git a/hexrd/core/transforms/new_capi/ndargs_helper.c b/hexrd/core/transforms/new_capi/ndargs_helper.c index bad63929b..5270d14dd 100644 --- a/hexrd/core/transforms/new_capi/ndargs_helper.c +++ b/hexrd/core/transforms/new_capi/ndargs_helper.c @@ -1,7 +1,6 @@ #include "ndargs_helper.h" -XRD_CFUNCTION int -nah_array_converter(PyObject *op, void *result) +static int nah_array_converter(PyObject *op, void *result) { nah_array *na = (nah_array*)result; diff --git a/hexrd/core/transforms/new_capi/ndargs_helper.h b/hexrd/core/transforms/new_capi/ndargs_helper.h index 6c7d867f3..4e43bcdb7 100644 --- a/hexrd/core/transforms/new_capi/ndargs_helper.h +++ b/hexrd/core/transforms/new_capi/ndargs_helper.h @@ -33,6 +33,6 @@ typedef struct nah_array_struct * returns 1 on success, 0 otherwise. When returning 0, a python * exception has been setup. */ -XRD_CFUNCTION int nah_array_converter(PyObject *op, void *result); +static int nah_array_converter(PyObject *op, void *result); #endif /* NDARRAY_ARGS_HELPER_H */ diff --git a/hexrd/core/transforms/new_capi/new_func.c b/hexrd/core/transforms/new_capi/new_func.c deleted file mode 100644 index 3f425b5f0..000000000 --- a/hexrd/core/transforms/new_capi/new_func.c +++ /dev/null @@ -1,22 +0,0 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -#endif - - -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *<>_docstring = -"c module implementation of <>. Please use the Python wrapper.\n" - - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ - diff --git a/hexrd/core/transforms/new_capi/oscill_angles_of_HKLs.c b/hexrd/core/transforms/new_capi/oscill_angles_of_HKLs.c index 22fb34539..e97a51881 100644 --- a/hexrd/core/transforms/new_capi/oscill_angles_of_HKLs.c +++ b/hexrd/core/transforms/new_capi/oscill_angles_of_HKLs.c @@ -1,13 +1,5 @@ -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -# include "ndargs_helper.h" -#endif - - -XRD_CFUNCTION void -oscill_angles_of_HKLs(size_t npts, double * hkls, double chi, +static void oscill_angles_of_HKLs(size_t npts, double * hkls, double chi, double * rMat_c, double * bMat, double wavelength, double * vInv_s, double * beamVec, double * etaVec, double * oangs0, double * oangs1) @@ -175,21 +167,11 @@ oscill_angles_of_HKLs(size_t npts, double * hkls, double chi, } } - -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_oscillAnglesOfHKLs = +static const char *docstring_oscillAnglesOfHKLs = "c module implementation of solve_omega.\n" "Please use the Python wrapper.\n"; -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_oscillAnglesOfHKLs(PyObject * self, PyObject * args) { nah_array hkls = { NULL, "hkls", NAH_TYPE_DP_FP, { 3, NAH_DIM_ANY }}; @@ -259,5 +241,3 @@ python_oscillAnglesOfHKLs(PyObject * self, PyObject * args) return PyErr_NoMemory(); } - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/quat_distance.c b/hexrd/core/transforms/new_capi/quat_distance.c index 2d1b17753..f385e050d 100644 --- a/hexrd/core/transforms/new_capi/quat_distance.c +++ b/hexrd/core/transforms/new_capi/quat_distance.c @@ -1,17 +1,5 @@ -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -/* TODO: No printf should go here. Also, it should be possible to avoid malloc - */ -# include -# include - -# include "transforms_utils.h" -# include "transforms_prototypes.h" -#endif - - -XRD_CFUNCTION double -quat_distance(size_t nsym, double * q1, double * q2, double * qsym) +static double quat_distance(size_t nsym, double * q1, double * q2, double * qsym) { size_t i; double q0_max = 0.0, dist = 0.0; @@ -47,21 +35,11 @@ quat_distance(size_t nsym, double * q1, double * q2, double * qsym) } -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# include "ndargs_helper.h" -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_quat_distance = +static const char *docstring_quat_distance = "c module implementation of quat_distance.\n" "Please use the Python wrapper.\n"; -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_quat_distance(PyObject * self, PyObject * args) { nah_array q1 = { NULL, "q1", NAH_TYPE_DP_FP, { 4 }}; @@ -84,6 +62,3 @@ python_quat_distance(PyObject * self, PyObject * args) return PyFloat_FromDouble(dist); } - - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/reference.py b/hexrd/core/transforms/new_capi/reference.py deleted file mode 100644 index 36dfe681e..000000000 --- a/hexrd/core/transforms/new_capi/reference.py +++ /dev/null @@ -1,57 +0,0 @@ -import numpy as np - -# Functions to use for reference in testing and the like. These are made to be -# easily testable and clear. They will be use to test against in unit tests. -# They may be slow and not vectorized. - - -def intersect_ray_plane(ro, rv, p): - ''' - ray-plane intersection - - Parameters - ---------- - ro: array[3] - point of origin of the ray (ray origin). - - rv: array[3] - direction vector of the ray (ray vector). - - p: array[4] - (A, B, C, D) coeficients for the plane formula (Ax+By+Cz+D=0). - - Returns - ------- - t : scalar - t where ray intersects with plane, such as (ro + t*rv) lies within p. - - Notes - ----- - If t is negative, the intersection happens 'behind' the origin point. - - In the case where (A,B,C) -the plane normal- is orthogonal to rv, no - intersection will happen (or the ray is fully on the plane). In that case, - t will be non-finite. - - Vector rv needs not to be an unit vector, result will be scaled accordingly - based on the modulo of rv. However, a 0 vector will result in non-finite - result. - - The code is based on the segment-plane intersection in [1]_ - - .. [1] Ericson. (2005). Real-Time Collision Detection, pp 176. - Morgan Kaufmann. - ''' - assert ro.shape == (3,) - assert rv.shape == (3,) - assert p.shape == (4,) - - normal = p[:3] - D = p[3] - # disable divide by 0 and invalid as the division in t can cause both. - # the behavior of the function actually relies in IEEE754 with a division - # by 0 generating the appropriate infinity, or a NAN if it is a 0/0. - with np.errstate(divide='ignore', invalid='ignore'): - t = (D - normal @ ro) / (normal @ rv) - - return t diff --git a/hexrd/core/transforms/new_capi/rotate_vecs_about_axis.c b/hexrd/core/transforms/new_capi/rotate_vecs_about_axis.c index 8b8af3955..5ebed5ec0 100644 --- a/hexrd/core/transforms/new_capi/rotate_vecs_about_axis.c +++ b/hexrd/core/transforms/new_capi/rotate_vecs_about_axis.c @@ -1,12 +1,4 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -#endif - - -XRD_CFUNCTION void -rotate_vecs_about_axis(size_t na, double *angles, +static void rotate_vecs_about_axis(size_t na, double *angles, size_t nax, double *axes, size_t nv, double *vecs, double *rVecs) @@ -54,22 +46,11 @@ rotate_vecs_about_axis(size_t na, double *angles, } } - -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# include "ndargs_helper.h" -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_rotate_vecs_about_axis = +static const char *docstring_rotate_vecs_about_axis = "c module implementation of rotate_vecs_about_axis.\n" "Please use the Python wrapper.\n"; -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_rotate_vecs_about_axis(PyObject *self, PyObject *args) { /* API interface in Python is: @@ -139,5 +120,3 @@ python_rotate_vecs_about_axis(PyObject *self, PyObject *args) return 0; } - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/transforms_prototypes.h b/hexrd/core/transforms/new_capi/transforms_prototypes.h index 82b6877e9..3c6e1931f 100644 --- a/hexrd/core/transforms/new_capi/transforms_prototypes.h +++ b/hexrd/core/transforms/new_capi/transforms_prototypes.h @@ -1,77 +1,73 @@ #ifndef XRD_TRANSFORMS_PROTOTYPES_H #define XRD_TRANSFORMS_PROTOTYPES_H -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_types.h" -#endif - -XRD_CFUNCTION void +static void angles_to_gvec(size_t nvecs, double *angs, double *bHat_l, double *eHat_l, double chi, double *rMat_c, double *gVec_c); -XRD_CFUNCTION void +static void angles_to_dvec(size_t nvecs, double *angs, double *bHat_l, double *eHat_l, double chi, double *rMat_c, double * gVec_c); #define GV2XY_SINGLE_RMAT_S 1 -XRD_CFUNCTION void +static void gvec_to_xy(size_t npts, const double *gVec_c_Ptr, const double *rMat_d_Ptr, const double *rMat_s_Ptr, const double *rMat_c_Ptr, const double *tVec_d_Ptr, const double *tVec_s_Ptr, const double *tVec_c_Ptr, const double *beamVec_Ptr, double * restrict result_Ptr, unsigned int flags); -XRD_CFUNCTION void +static void xy_to_gvec(size_t npts, double *xy, double *rMat_d, double *rMat_s, double *tVec_d, double *tVec_s, double *tVec_c, double *beamVec, double *etaVec, double *tTh, double *eta, double *gVec_l); -XRD_CFUNCTION void +static void oscill_angles_of_HKLs(size_t npts, double *hkls, double chi, double *rMat_c, double *bMat, double wavelength, double *vInv_s, double *beamVec, double *etaVec, double *oangs0, double *oangs1); /* this should probably be just a utility function in util... */ -XRD_CFUNCTION int +static int unit_row_vector(size_t n, double *cIn, double *cOut); -XRD_CFUNCTION void +static void unit_row_vectors(size_t m, size_t n, double *cIn, double *cOut); /* -XRD_CFUNCTION void +static void make_detector_rmat(double *tPtr, double *rPtr); */ -XRD_CFUNCTION void +static void make_sample_rmat(double chi, const double ome, double *result_rmat); -XRD_CFUNCTION void +static void make_sample_rmat_array(double chi, const double *ome_ptr, size_t ome_count, double *result_ptr); -XRD_CFUNCTION void +static void make_rmat_of_expmap(double *ePtr, double *rPtr); -XRD_CFUNCTION void +static void make_binary_rmat(const double *aPtr, double * restrict rPtr); #define TF_MAKE_BEAM_RMAT_ERR_BEAM_ZERO 1 #define TF_MAKE_BEAM_RMAT_ERR_COLLINEAR 2 -XRD_CFUNCTION int +static int make_beam_rmat(double *bPtr, double *ePtr, double *rPtr); /* aka make_eta_frame_rotmat */ -XRD_CFUNCTION void +static void validate_angle_ranges(size_t na, double *aPtr, size_t nr, double *minPtr, double *maxPtr, bool *rPtr, int ccw); -XRD_CFUNCTION void +static void rotate_vecs_about_axis(size_t na, double *angles, size_t nax, double *axes, size_t nv, double * vecs, double * rVecs); -XRD_CFUNCTION double +static double quat_distance(size_t nsym, double *q1, double *q2, double *qsym); diff --git a/hexrd/core/transforms/new_capi/unit_row_vector.c b/hexrd/core/transforms/new_capi/unit_row_vector.c index e402114be..666b90fa5 100644 --- a/hexrd/core/transforms/new_capi/unit_row_vector.c +++ b/hexrd/core/transforms/new_capi/unit_row_vector.c @@ -1,12 +1,4 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -# include "ndargs_helper.h" -#endif - - -XRD_CFUNCTION int +static int unit_row_vector(size_t n, double * cIn, double * cOut) { size_t j; @@ -31,7 +23,7 @@ unit_row_vector(size_t n, double * cIn, double * cOut) } } -XRD_CFUNCTION void +static void unit_row_vectors(size_t m, size_t n, double *cIn, double *cOut) { size_t i, j; @@ -55,25 +47,15 @@ unit_row_vectors(size_t m, size_t n, double *cIn, double *cOut) } } - -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_unitRowVector = +static const char *docstring_unitRowVector = "c module implementation of unit_row_vector (one row).\n" "Please use the Python wrapper.\n"; -XRD_PYTHON_WRAPPER const char *docstring_unitRowVectors = +static const char *docstring_unitRowVectors = "c module implementation of unit_row_vector (multiple rows).\n" "Please use the Python wrapper.\n"; -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_unitRowVector(PyObject * self, PyObject * args) { PyArrayObject *aop_out = NULL; @@ -110,7 +92,7 @@ python_unitRowVector(PyObject * self, PyObject * args) } -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_unitRowVectors(PyObject *self, PyObject *args) { PyArrayObject *aop_out = NULL; @@ -140,6 +122,3 @@ python_unitRowVectors(PyObject *self, PyObject *args) return (PyObject*)aop_out; } - - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/new_capi/validate_angle_ranges.c b/hexrd/core/transforms/new_capi/validate_angle_ranges.c index c87ffd9ba..d4896723a 100644 --- a/hexrd/core/transforms/new_capi/validate_angle_ranges.c +++ b/hexrd/core/transforms/new_capi/validate_angle_ranges.c @@ -1,10 +1,4 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -#endif - -XRD_CFUNCTION void +static void validate_angle_ranges(size_t na, double *aPtr, size_t nr, double *minPtr, double *maxPtr, bool *rPtr, int ccw) { @@ -77,22 +71,11 @@ validate_angle_ranges(size_t na, double *aPtr, size_t nr, double *minPtr, } } - -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# include "ndargs_helper.h" -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_validateAngleRanges = +static const char *docstring_validateAngleRanges = "c module implementation of validate_angle_ranges.\n" "Please use the Python wrapper.\n"; -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_validateAngleRanges(PyObject * self, PyObject * args) { nah_array ang_list = { NULL, "ang_list", NAH_TYPE_DP_FP, { NAH_DIM_ANY }}; @@ -145,6 +128,3 @@ python_validateAngleRanges(PyObject * self, PyObject * args) Py_XDECREF(result); return NULL; } - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ - diff --git a/hexrd/core/transforms/new_capi/xf_new_capi.py b/hexrd/core/transforms/new_capi/xf_new_capi.py deleted file mode 100644 index 4e41ed299..000000000 --- a/hexrd/core/transforms/new_capi/xf_new_capi.py +++ /dev/null @@ -1,691 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Transforms module implementation using a support C extension module. - -Currently, this implementation contains code for the following functions: - - - angles_to_gvec - - angles_to_dvec - - gvec_to_xy - - xy_to_gvec (partial) - - - unit_vector - - make_rmat_of_exp_map - - make_binary_rmat - - make_beam_rmat - - validate_angle_ranges - - rotate_vecs_about_axis - - quat_distance - -There are also some functions that maybe would be needed in the transforms -module: - - makeGVector - - makeRotMatOfQuat - - homochoricOfQuat -""" -from typing import Optional, Tuple, Sequence -import numpy as np -from numpy.typing import NDArray - -from hexrd.core.extensions import _new_transforms_capi as _impl -from hexrd.core.extensions import transforms as cpp_transforms -from hexrd.core.distortion.distortionabc import DistortionABC -from hexrd.core import constants as cnst - - -def angles_to_gvec( - angs: NDArray[np.float64], - beam_vec: NDArray[np.float64] = cnst.beam_vec, - eta_vec: NDArray[np.float64] = cnst.eta_vec, - chi: float = 0.0, - rmat_c: NDArray[np.float64] = cnst.identity_3x3, -) -> NDArray[np.float64]: - """ - - Takes triplets of angles in the beam frame (2*theta, eta[, omega]) - to components of unit G-vectors in the LAB frame. If the omega - values are not trivial (i.e. angs[:, 2] = 0.), then the components - are in the SAMPLE frame. If the crystal rmat is specified and - is not the identity, then the components are in the CRYSTAL frame. - - G vectors here referes to the reciprocal lattice vectors. - - Parameters - ---------- - angs : ndarray - The euler angles of diffraction. The last dimension must be 2 or 3. In - (2*theta, eta[, omega]) format. - beam_vec : ndarray, optional - Unit vector pointing towards the X-ray source in the lab frame. - Defaults to [0,0,-1] - eta_vec : ndarray, optional - Vector defining eta=0 in the lab frame. Defaults to [1,0,0] - chi : float, optional - The inclination angle of the sample frame about the lab frame X. - rmat_c : ndarray, optional - The change of basis matrix from the reciprocal frame to the crystal - frame. Defaults to the identity. - - Returns - ------- - ndarray - (3,n) array of unit reciprocal lattice vectors, frame depends on the - input parameters - """ - - orig_ndim = angs.ndim - - # if only a pair is provided... convert to a triplet with omegas == 0 - # so that behavior is preserved. - if angs.shape[-1] == 2: - angs = np.hstack((angs, np.zeros(angs.shape[:-1] + (1,)))) - - angs = np.ascontiguousarray(np.atleast_2d(angs)) - beam_vec = np.ascontiguousarray(beam_vec.flatten()) - eta_vec = np.ascontiguousarray(eta_vec.flatten()) - rmat_c = np.ascontiguousarray(rmat_c) - - result = cpp_transforms.anglesToGVec(angs, beam_vec, eta_vec, chi, rmat_c) - - return result[0] if orig_ndim == 1 else result - - -def angles_to_dvec( - angs: NDArray[np.float64], - beam_vec: NDArray[np.float64] = cnst.beam_vec, - eta_vec: NDArray[np.float64] = cnst.eta_vec, - chi: float = 0.0, - rmat_c: NDArray[np.float64] = cnst.identity_3x3, -) -> NDArray[np.float64]: - """Calculate diffraction vectors from beam frame angles. - - Takes triplets of angles in the beam frame (2*theta, eta[, omega]) - to components of unit diffraction vectors in the LAB frame. If the omega - values are not trivial (i.e. angs[:, 2] = 0.), then the components - are in the SAMPLE frame. If the crystal rmat is specified and - is not the identity, then the components are in the CRYSTAL frame. - - - Parameters - ---------- - angs : ndarray - The euler angles of diffraction. The last dimension must be 2 or 3. In - (2*theta, eta[, omega]) format. - beam_vec : ndarray, optional - Unit vector pointing towards the X-ray source in the lab frame. - Defaults to [0,0,-1] - eta_vec : ndarray, optional - Vector defining eta=0 in the lab frame. Defaults to [1,0,0] - chi : float, optional - The inclination angle of the sample frame about the lab frame X. - rmat_c : ndarray, optional - The change of basis matrix from the reciprocal frame to the crystal - frame. Defaults to the identity. - - Returns - ------- - ndarray - (3,n) array of unit diffraction vectors, frame depends on the input - parameters - """ - # TODO: Improve capi to avoid multiplications when rmat_c is None - - # if only a pair is provided... convert to a triplet with omegas == 0 - # so that behavior is preserved. - if angs.shape[-1] == 2: - angs = np.hstack((angs, np.zeros(angs.shape[:-1] + (1,)))) - - angs = np.ascontiguousarray(np.atleast_2d(angs)) - beam_vec = np.ascontiguousarray(beam_vec.flatten()) - eta_vec = np.ascontiguousarray(eta_vec.flatten()) - rmat_c = np.ascontiguousarray(rmat_c) - - return cpp_transforms.anglesToDVec(angs, beam_vec, eta_vec, chi, rmat_c) - - -def makeGVector( - hkl: NDArray[np.float64], bMat: NDArray[np.float64] -) -> NDArray[np.float64]: - """ - Take a crystal relative b matrix onto a list of hkls to output unit - reciprocal latice vectors (a.k.a. lattice plane normals) - - - Parameters - ---------- - hkl : ndarray - (3,n) ndarray of n hstacked reciprocal lattice vector component - triplets - bMat : ndarray - (3,3) ndarray of the change of basis matrix from the reciprocal - lattice to the crystal reference frame - - Returns - ------- - ndarray - (3,n) ndarray of n unit reciprocal lattice vectors (a.k.a. lattice - plane normals) - - """ - if hkl.shape[0] != 3: - raise ValueError('hkl input must be (3, n)') - return unit_vector(np.dot(bMat, hkl)) - - -def gvec_to_xy( - gvec_c: NDArray[np.float64], - rmat_d: NDArray[np.float64], - rmat_s: NDArray[np.float64], - rmat_c: NDArray[np.float64], - tvec_d: NDArray[np.float64], - tvec_s: NDArray[np.float64], - tvec_c: NDArray[np.float64], - beam_vec: Optional[NDArray[np.float64]] = None, - vmat_inv: Optional[NDArray[np.float64]] = None, - bmat: Optional[NDArray[np.float64]] = None, -) -> NDArray[np.float64]: - """ - Takes a concatenated list of reciprocal lattice vectors components in the - CRYSTAL FRAME to the specified detector-relative frame, subject to the - following: - - 1) it must be able to satisfy a bragg condition - 2) the associated diffracted beam must intersect the detector plane - - Parameters - ---------- - gvec_c : array_like - ([N,] 3) G-vector components in the CRYSTAL FRAME. - rmat_d : array_like - The (3, 3) COB matrix taking components in the - DETECTOR FRAME to the LAB FRAME - rmat_s : array_like - The ([N,] 3, 3) COB matrix taking components in the SAMPLE FRAME to the - LAB FRAME. It may be a single (3, 3) rotation matrix to use for all - gvec_c, or just one rotation matrix per gvec. - rmat_c : array_like - The (3, 3) COB matrix taking components in the - CRYSTAL FRAME to the SAMPLE FRAME - tvec_d : array_like - The (3, ) translation vector connecting LAB FRAME to DETECTOR FRAME - tvec_s : array_like - The (3, ) translation vector connecting LAB FRAME to SAMPLE FRAME - tvec_c : array_like - The ([M,] 3, ) translation vector(s) connecting SAMPLE FRAME to - CRYSTAL FRAME - beam_vec : array_like, optional - The (3, ) incident beam propagation vector components in the LAB FRAME; - the default is [0, 0, -1], which is the standard setting. - vmat_inv : array_like, optional - The (3, 3) matrix of inverse stretch tensor components in the - SAMPLE FRAME. The default is None, which implies a strain-free state - (i.e. V = I). - bmat : array_like, optional - The (3, 3) COB matrix taking components in the - RECIPROCAL LATTICE FRAME to the CRYSTAL FRAME; if supplied, it is - assumed that the input `gvecs` are G-vector components in the - RECIPROCL LATTICE FRAME (the default is None, which implies components - in the CRYSTAL FRAME) - - Returns - ------- - array_like - The ([M, ][N, ] 2) array of [x, y] diffracted beam intersections for - each of the N input G-vectors in the DETECTOR FRAME (all Z_d - coordinates are 0 and excluded) and for each of the M candidate - positions. For each input G-vector that cannot satisfy a Bragg - condition or intersect the detector plane, [NaN, Nan] is returned. - - Notes - ----- - Previously only a single candidate position was allowed. This is in - fact a vectored version of the previous API function. It is backwards - compatible, as passing single tvec_c is supported and has the same - result. - - """ - beam_vec = beam_vec if beam_vec is not None else cnst.beam_vec - - orig_ndim = gvec_c.ndim - gvec_c = np.ascontiguousarray(np.atleast_2d(gvec_c)) - rmat_d = np.ascontiguousarray(rmat_d) - rmat_s = np.ascontiguousarray(rmat_s) - rmat_c = np.ascontiguousarray(rmat_c) - tvec_d = np.ascontiguousarray(tvec_d.flatten()) - tvec_s = np.ascontiguousarray(tvec_s.flatten()) - tvec_c = np.ascontiguousarray(tvec_c.flatten()) - beam_vec = np.ascontiguousarray(beam_vec.flatten()) - - # depending on the number of dimensions of rmat_s use either the array - # version or the "scalar" (over rmat_s) version. Note that rmat_s is either - # a 3x3 matrix (ndim 2) or an nx3x4 array of matrices (ndim 3) - if rmat_s.ndim > 2: - result = _impl.gvecToDetectorXYArray( - gvec_c, rmat_d, rmat_s, rmat_c, tvec_d, tvec_s, tvec_c, beam_vec - ) - else: - result = _impl.gvecToDetectorXY( - gvec_c, rmat_d, rmat_s, rmat_c, tvec_d, tvec_s, tvec_c, beam_vec - ) - return result[0] if orig_ndim == 1 else result - - -def xy_to_gvec( - xy_d: NDArray[np.float64], - rmat_d: NDArray[np.float64], - rmat_s: NDArray[np.float64], - tvec_d: NDArray[np.float64], - tvec_s: NDArray[np.float64], - tvec_c: NDArray[np.float64], - rmat_b: Optional[NDArray[np.float64]] = None, - distortion: Optional[DistortionABC] = None, - output_ref: Optional[bool] = False, -) -> tuple[NDArray[np.float64], ...]: - """ - Takes a list cartesian (x, y) pairs in the DETECTOR FRAME and calculates - the associated reciprocal lattice (G) vectors and (bragg angle, azimuth) - pairs with respect to the specified beam and azimth (eta) reference - directions. - - Parameters - ---------- - xy_d : array_like - (n, 2) array of n (x, y) coordinates in DETECTOR FRAME - rmat_d : array_like - (3, 3) COB matrix taking components in the - DETECTOR FRAME to the LAB FRAME - rmat_s : array_like - (3, 3) COB matrix taking components in the - SAMPLE FRAME to the LAB FRAME - tvec_d : array_like - (3, ) translation vector connecting LAB FRAME to DETECTOR FRAME - tvec_s : array_like - (3, ) translation vector connecting LAB FRAME to SAMPLE FRAME - tvec_c : array_like - (3, ) translation vector connecting SAMPLE FRAME to CRYSTAL FRAME - rmat_b : array_like, optional - (3, 3) COB matrix taking components in the BEAM FRAME to the LAB FRAME; - defaults to None, which implies the standard setting of identity. - distortion : distortion class, optional - Default is None - output_ref : bool, optional - If True, prepends the apparent bragg angle and azimuth with respect to - the SAMPLE FRAME (ignoring effect of non-zero tvec_c) - - Returns - ------- - array_like - (n, 2) ndarray containing the (tth, eta) pairs associated with each - (x, y) associated with gVecs - array_like - (n, 3) ndarray containing the associated G vector directions in the - LAB FRAME - array_like, optional - if output_ref is True - """ - # It seems that the output_ref version is not present as the argument - # gets ignored - - rmat_b = rmat_b if rmat_b is not None else cnst.identity_3x3 - - # the code seems to ignore this argument, assume output_ref == True not - # implemented - assert not output_ref, 'output_ref not implemented' - - if distortion is not None: - xy_d = distortion.apply_inverse(xy_d) - - xy_d = np.ascontiguousarray(np.atleast_2d(xy_d)) - rmat_d = np.ascontiguousarray(rmat_d) - rmat_s = np.ascontiguousarray(rmat_s) - tvec_d = np.ascontiguousarray(tvec_d.flatten()) - tvec_s = np.ascontiguousarray(tvec_s.flatten()) - tvec_c = np.ascontiguousarray(tvec_c.flatten()) - beam_vec = np.ascontiguousarray((-rmat_b[:, 2]).flatten()) - eta_vec = np.ascontiguousarray(rmat_b[:, 0].flatten()) - return _impl.detectorXYToGvec( - xy_d, rmat_d, rmat_s, tvec_d, tvec_s, tvec_c, beam_vec, eta_vec - ) - - -def oscill_angles_of_hkls( - hkls: NDArray[np.float64], - chi: float, - rmat_c: NDArray[np.float64], - bmat: NDArray[np.float64], - wavelength: float, - v_inv: Optional[NDArray[np.float64]] = None, - beam_vec: NDArray[np.float64] = cnst.beam_vec, - eta_vec: NDArray[np.float64] = cnst.eta_vec, -) -> Tuple[NDArray[np.float64], NDArray[np.float64]]: - """ - Takes a list of unit reciprocal lattice vectors in crystal frame to the - specified detector-relative frame, subject to the conditions: - - 1) the reciprocal lattice vector must be able to satisfy a bragg condition - 2) the associated diffracted beam must intersect the detector plane - - Parameters - ---------- - hkls -- (n, 3) ndarray of n reciprocal lattice vectors - in the CRYSTAL FRAME - chi -- float representing the inclination angle of the - oscillation axis (std coords) - rmat_c -- (3, 3) ndarray, the COB taking CRYSTAL FRAME components - to SAMPLE FRAME - bmat -- (3, 3) ndarray, the COB taking RECIPROCAL LATTICE components - to CRYSTAL FRAME - wavelength -- float representing the x-ray wavelength in Angstroms - - Optional Keyword Arguments: - v_inv -- (6, ) ndarray, vec representation inverse stretch tensor - components in the SAMPLE FRAME. The default is None, which - implies a strain-free state (i.e. V = I). - beam_vec -- (3, 1) mdarray containing the incident beam direction - components in the LAB FRAME - eta_vec -- (3, 1) mdarray containing the reference azimuth direction - components in the LAB FRAME - - Outputs: - ome0 -- (n, 3) ndarray containing the feasible (tTh, eta, ome) triplets - for each input hkl (first solution) - ome1 -- (n, 3) ndarray containing the feasible (tTh, eta, ome) triplets - for each input hkl (second solution) - - Notes: - ------------------------------------------------------------------------ - The reciprocal lattice vector, G, will satisfy the the Bragg condition - when: - - b.T * G / ||G|| = -sin(theta) - - where b is the incident beam direction (k_i) and theta is the Bragg - angle consistent with G and the specified wavelength. The components of - G in the lab frame in this case are obtained using the crystal - orientation, Rc, and the single-parameter oscillation matrix, Rs(ome): - - Rs(ome) * Rc * G / ||G|| - - The equation above can be rearranged to yield an expression of the form: - - a*sin(ome) + b*cos(ome) = c - - which is solved using the relation: - - a*sin(x) + b*cos(x) = sqrt(a**2 + b**2) * sin(x + alpha) - - --> sin(x + alpha) = c / sqrt(a**2 + b**2) - - where: - - alpha = atan2(b, a) - - The solutions are: - - / - | arcsin(c / sqrt(a**2 + b**2)) - alpha - x = < - | pi - arcsin(c / sqrt(a**2 + b**2)) - alpha - \ - - There is a double root in the case the reflection is tangent to the - Debye-Scherrer cone (c**2 = a**2 + b**2), and no solution if the - Laue condition cannot be satisfied (filled with NaNs in the results - array here) - """ - # this was oscillAnglesOfHKLs - hkls = np.array(hkls, dtype=float, order='C') - if v_inv is None: - v_inv = np.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]) - else: - v_inv = np.ascontiguousarray(v_inv.flatten()) - rmat_c = np.ascontiguousarray(rmat_c) - beam_vec = np.ascontiguousarray(beam_vec.flatten()) - eta_vec = np.ascontiguousarray(eta_vec.flatten()) - bmat = np.ascontiguousarray(bmat) - return _impl.oscillAnglesOfHKLs( - hkls, chi, rmat_c, bmat, wavelength, v_inv, beam_vec, eta_vec - ) - - -def unit_vector(vec_in: NDArray[np.float64]) -> NDArray[np.float64]: - """ - Normalize the input vector(s) to unit length. - - Parameters - ---------- - vec_in : ndarray - The input vector(s) (n,3) to normalize. - - Returns - ------- - ndarray - The normalized vector(s) of the same shape as the input. - - Raises - ------ - ValueError - If the input vector(s) do not have the correct shape. - """ - vec_in = np.ascontiguousarray(vec_in) - if vec_in.ndim == 1: - return _impl.unitRowVector(vec_in) - elif vec_in.ndim == 2: - return _impl.unitRowVectors(vec_in) - else: - raise ValueError( - "incorrect arg shape; must be 1-d or 2-d, yours is %d-d" % (vec_in.ndim) - ) - - -def make_detector_rmat(tilt_angles: NDArray[np.float64]) -> NDArray[np.float64]: - """ - Form the (3, 3) tilt rotations from the tilt angle list: - - tilt_angles = [gamma_Xl, gamma_Yl, gamma_Zl] in radians - Returns the (3, 3) rotation matrix from the detector frame to the lab frame - """ - t_angles = np.ascontiguousarray(np.r_[tilt_angles].flatten()) - return _impl.makeDetectorRotMat(t_angles) - - -# make_sample_rmat in CAPI is split between makeOscillRotMat -# and makeOscillRotMatArray... - - -def make_oscill_rmat(oscillAngles): - chi, ome = oscillAngles - ome = np.atleast_1d(ome) - result = _impl.makeOscillRotMat(chi, ome) - return result.reshape((3, 3)) - - -def make_oscill_rmat_array(chi, omeArray): - arg = np.ascontiguousarray(omeArray) - return _impl.makeOscillRotMat(chi, arg) - - -def make_sample_rmat( - chi: float, ome: float | NDArray[np.float64] -) -> NDArray[np.float64]: - # TODO: Check this docstring - """ - Make a rotation matrix representing the COB from sample frame to the lab - frame. - - Parameters - ---------- - chi : float - The inclination angle of the sample frame about the lab frame Y. - ome : float or ndarray - The oscillation angle(s) about the sample frame Y. - - Returns - ------- - ndarray - A 3x3 rotation matrix representing the input oscillation angles. - - """ - ome_array = np.atleast_1d(ome) - if ome is ome_array: - ome_array = np.ascontiguousarray(ome_array) - result = cpp_transforms.makeOscillRotMat(chi, ome_array).reshape( - len(ome_array), 3, 3 - ) - else: - # converted to 1d array of 1 element, no need - # to call ascontiguousarray, but need to remove - # the outer dimension from the result - result = cpp_transforms.makeOscillRotMat(chi, ome_array) - - return result - - -def make_rmat_of_expmap( - exp_map: Sequence[float] | NDArray[np.float64], -) -> NDArray[np.float64]: - """ - Calculate the rotation matrix of an exponential map. - - Parameters - ---------- - exp_map : array_like - A 3-element sequence representing the exponential map n*phi. - - Returns - ------- - NDArray[np.float64] - A 3x3 rotation matrix representing the input exponential map - """ - arg = np.ascontiguousarray(exp_map).flatten() - return cpp_transforms.make_rot_mat_of_exp_map(arg) - - -def make_binary_rmat(axis: NDArray[np.float64]) -> NDArray[np.float64]: - """ - Create a 180 degree rotation matrix about the input axis. - - Parameters - ---------- - axis : ndarray - 3-element sequence representing the rotation axis, in radians. - - Returns - ------- - ndarray - A 3x3 rotation matrix representing a 180 degree rotation about the - input axis. - """ - - arg = np.ascontiguousarray(axis.flatten()) - return cpp_transforms.make_binary_rot_mat(arg) - - -def make_beam_rmat( - bvec_l: NDArray[np.float64], evec_l: NDArray[np.float64] -) -> NDArray[np.float64]: - """ - Creates a COB matrix from the beam frame to the lab frame - Note: beam and eta vectors must not be colinear - - Parameters - ---------- - bvec_l : ndarray - The beam vector in the lab frame, The (3, ) incident beam propagation - vector components in the lab frame. - the default is [0, 0, -1], which is the standard setting. - evec_l : ndarray - Vector defining eta=0 in the lab frame. Defaults to [1,0,0] - """ - arg1 = np.ascontiguousarray(bvec_l.flatten()) - arg2 = np.ascontiguousarray(evec_l.flatten()) - return _impl.makeEtaFrameRotMat(arg1, arg2) - - -def validate_angle_ranges( - ang_list: float | NDArray[np.float64], - start_angs: float | NDArray[np.float64], - stop_angs: float | NDArray[np.float64], - ccw: bool = True, -) -> NDArray[np.bool_]: - """ - Find out if angles are in the CCW or CW range from start to stop - - Parameters - ---------- - ang_list : ndarray - The list of angles to validate - start_angs : ndarray - The starting angles - stop_angs : ndarray - The stopping angles - ccw : bool, optional - If True, the angles are in the CCW range from start to stop. - Defaults to True. - - Returns - ------- - ndarray - List of bools indicating if the angles are in the correct range - - """ - ang_list = np.atleast_1d(ang_list).astype(np.double, order='C') - start_angs = np.atleast_1d(start_angs).astype(np.double, order='C') - stop_angs = np.atleast_1d(stop_angs).astype(np.double, order='C') - - return _impl.validateAngleRanges(ang_list, start_angs, stop_angs, ccw) - - -def rotate_vecs_about_axis( - angle: float | NDArray[np.float64], - axis: NDArray[np.float64], - vecs: NDArray[np.float64], -) -> NDArray[np.float64]: - """ - Rotate vectors about an axis - - Parameters - ---------- - angle : array_like - Array of angles (len==1 or n) - axis : ndarray - Array of unit vectors (shape == (3,1) or (3,n)) - vecs : ndarray - Array of vectors to rotate (shape == (3,n)) - - Returns - ------- - ndarray - Array of rotated vectors (shape == (3,n)) - """ - angle = np.asarray(angle) - axis = np.ascontiguousarray(axis.T) - vecs = np.ascontiguousarray(vecs.T) - result = _impl.rotate_vecs_about_axis(angle, axis, vecs) - return result.T - - -def quat_distance( - q1: NDArray[np.float64], q2: NDArray[np.float64], qsym: NDArray[np.float64] -) -> NDArray[np.float64]: - """Distance between two quaternions, taking quaternions of symmetry into account. - - Parameters - ---------- - q1 : arary_like - First quaternion. - q2 : arary_like - Second quaternion. - qsym : ndarray - List of symmetry quaternions. - - Returns - ------- - float - The distance between the two quaternions. - """ - q1 = np.ascontiguousarray(q1.flatten()) - q2 = np.ascontiguousarray(q2.flatten()) - # C module expects quaternions in row major, numpy code in column major. - qsym = np.ascontiguousarray(qsym.T) - return _impl.quat_distance(q1, q2, qsym) diff --git a/hexrd/core/transforms/new_capi/xy_to_gvec.c b/hexrd/core/transforms/new_capi/xy_to_gvec.c index 9b7d546f2..611a83985 100644 --- a/hexrd/core/transforms/new_capi/xy_to_gvec.c +++ b/hexrd/core/transforms/new_capi/xy_to_gvec.c @@ -1,12 +1,4 @@ - -#if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# include "transforms_utils.h" -# include "transforms_prototypes.h" -# include "ndargs_helper.h" -#endif - - -XRD_CFUNCTION void +static void xy_to_gvec(size_t npts, double *xy, double *rMat_d, double *rMat_s, double *tVec_d, double *tVec_s, double *tVec_c, double *beamVec, double *etaVec, @@ -90,17 +82,7 @@ xy_to_gvec(size_t npts, double *xy, double *rMat_d, double *rMat_s, } } - -#if defined(XRD_INCLUDE_PYTHON_WRAPPERS) && XRD_INCLUDE_PYTHON_WRAPPERS - -# if !defined(XRD_SINGLE_COMPILE_UNIT) || !XRD_SINGLE_COMPILE_UNIT -# define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -# include -# include -# endif /* XRD_SINGLE_COMPILE_UNIT */ - -XRD_PYTHON_WRAPPER const char *docstring_detectorXYToGvec = +static const char *docstring_detectorXYToGvec = "c module implementation of xy_to_gvec.\n" "Please use the Python wrapper.\n"; @@ -127,7 +109,7 @@ XRD_PYTHON_WRAPPER const char *docstring_detectorXYToGvec = gvec_l -- (n, 3) ndarray containing the associated G vector directions in the LAB FRAME associated with gVecs */ -XRD_PYTHON_WRAPPER PyObject * +static PyObject * python_detectorXYToGvec(PyObject * self, PyObject * args) { /* Right now, the Python wrapper guarantees that: @@ -221,5 +203,3 @@ python_detectorXYToGvec(PyObject * self, PyObject * args) return PyErr_NoMemory(); } - -#endif /* XRD_INCLUDE_PYTHON_WRAPPERS */ diff --git a/hexrd/core/transforms/old_xfcapi.py b/hexrd/core/transforms/old_xfcapi.py deleted file mode 100644 index 0bbaa6ebf..000000000 --- a/hexrd/core/transforms/old_xfcapi.py +++ /dev/null @@ -1,548 +0,0 @@ -#! /usr/bin/env python -# ============================================================ -# Copyright (c) 2012, Lawrence Livermore National Security, LLC. -# Produced at the Lawrence Livermore National Laboratory. -# Written by Joel Bernier and others. -# LLNL-CODE-529294. -# All rights reserved. -# -# This file is part of HEXRD. For details on dowloading the source, -# see the file COPYING. -# -# Please also see the file LICENSE. -# -# This program is free software; you can redistribute it and/or modify it under -# the terms of the GNU Lesser General Public License (as published by the Free -# Software Foundation) version 2.1 dated February 1999. -# -# This program is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the IMPLIED WARRANTY OF MERCHANTABILITY -# or FITNESS FOR A PARTICULAR PURPOSE. See the terms and conditions of the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU Lesser General Public -# License along with this program (see file LICENSE); if not, write to -# the Free Software Foundation, Inc., 59 Temple Place, Suite 330, -# Boston, MA 02111-1307 USA or visit . -# ============================================================ - -import numpy as np -import sys - -from hexrd.core.extensions import _transforms_CAPI - -# Imports so that others can import from this module -from hexrd.core.rotations import mapAngle -from hexrd.core.matrixutil import columnNorm, rowNorm - -# ###################################################################### -# Module Data -epsf = np.finfo(float).eps # ~2.2e-16 -ten_epsf = 10 * epsf # ~2.2e-15 -sqrt_epsf = np.sqrt(epsf) # ~1.5e-8 - -periodDict = {'degrees': 360.0, 'radians': 2 * np.pi} -angularUnits = 'radians' # module-level angle units - -# basis vectors -I3 = np.eye(3) # (3, 3) identity -Xl = np.ascontiguousarray(I3[:, 0].reshape(3, 1)) # X in the lab frame -Yl = np.ascontiguousarray(I3[:, 1].reshape(3, 1)) # Y in the lab frame -Zl = np.ascontiguousarray(I3[:, 2].reshape(3, 1)) # Z in the lab frame - -# reference stretch -vInv_ref = np.array([[1.0, 1.0, 1.0, 0.0, 0.0, 0.0]], order='C').T - -# reference beam direction and eta=0 ref in LAB FRAME for standard geometry -bVec_ref = -Zl -eta_ref = Xl - -# ###################################################################### -# Funtions - - -def anglesToGVec(angs, bHat_l=bVec_ref, eHat_l=eta_ref, chi=0.0, rMat_c=I3): - """ - from 'eta' frame out to lab (with handy kwargs to go to crystal or sample) - - * setting omega to zero in ang imput and omitting rMat_c leaves - in the lab frame in accordance with beam frame specs. - """ - angs = np.ascontiguousarray(np.atleast_2d(angs)) - bHat_l = np.ascontiguousarray(bHat_l.flatten()) - eHat_l = np.ascontiguousarray(eHat_l.flatten()) - rMat_c = np.ascontiguousarray(rMat_c) - chi = float(chi) - return _transforms_CAPI.anglesToGVec(angs, bHat_l, eHat_l, chi, rMat_c) - - -def anglesToDVec(angs, bHat_l=bVec_ref, eHat_l=eta_ref, chi=0.0, rMat_c=I3): - """ - from 'eta' frame out to lab (with handy kwargs to go to crystal or sample) - - * setting omega to zero in ang imput and omitting rMat_c leaves - in the lab frame in accordance with beam frame specs. - """ - angs = np.ascontiguousarray(np.atleast_2d(angs)) - bHat_l = np.ascontiguousarray(bHat_l.flatten()) - eHat_l = np.ascontiguousarray(eHat_l.flatten()) - rMat_c = np.ascontiguousarray(rMat_c) - chi = float(chi) - return _transforms_CAPI.anglesToDVec(angs, bHat_l, eHat_l, chi, rMat_c) - - -def makeGVector(hkl, bMat): - """ - take a CRYSTAL RELATIVE B matrix onto a list of hkls to output unit - reciprocal lattice vectors (a.k.a. lattice plane normals) - - Required Arguments: - hkls -- (3, n) ndarray of n hstacked reciprocal lattice vector component - triplets - bMat -- (3, 3) ndarray representing the matirix taking reciprocal lattice - vectors to the crystal reference frame - - Output: - gVecs -- (3, n) ndarray of n unit reciprocal lattice vectors - (a.k.a. lattice plane normals) - - To Do: - * might benefit from some assert statements to catch improperly shaped - input. - """ - assert hkl.shape[0] == 3, 'hkl input must be (3, n)' - return unitRowVector(np.dot(bMat, hkl)) - - -def gvecToDetectorXY( - gVec_c, rMat_d, rMat_s, rMat_c, tVec_d, tVec_s, tVec_c, beamVec=bVec_ref -): - """ - Takes a list of unit reciprocal lattice vectors in crystal frame to the - specified detector-relative frame, subject to the conditions: - - 1) the reciprocal lattice vector must be able to satisfy a bragg condition - 2) the associated diffracted beam must intersect the detector plane - - Required Arguments: - gVec_c -- (n, 3) ndarray of n reciprocal lattice vector components - in the CRYSTAL FRAME - rMat_d -- (3, 3) ndarray, the COB taking DETECTOR FRAME components - to LAB FRAME - rMat_s -- (3, 3) ndarray, the COB taking SAMPLE FRAME components - to LAB FRAME - rMat_c -- (3, 3) ndarray, the COB taking CRYSTAL FRAME components - to SAMPLE FRAME - tVec_d -- (3, 1) ndarray, the translation vector connecting - LAB to DETECTOR - tVec_s -- (3, 1) ndarray, the translation vector connecting - LAB to SAMPLE - tVec_c -- (3, 1) ndarray, the translation vector connecting - SAMPLE to CRYSTAL - - Outputs: - (m, 2) ndarray containing the intersections of m <= n diffracted beams - associated with gVecs - """ - rMat_d = np.ascontiguousarray(rMat_d) - rMat_s = np.ascontiguousarray(rMat_s) - rMat_c = np.ascontiguousarray(rMat_c) - gVec_c = np.ascontiguousarray(np.atleast_2d(gVec_c)) - tVec_d = np.ascontiguousarray(tVec_d.flatten()) - tVec_s = np.ascontiguousarray(tVec_s.flatten()) - tVec_c = np.ascontiguousarray(tVec_c.flatten()) - beamVec = np.ascontiguousarray(beamVec.flatten()) - return _transforms_CAPI.gvecToDetectorXY( - gVec_c, rMat_d, rMat_s, rMat_c, tVec_d, tVec_s, tVec_c, beamVec - ) - - -def gvecToDetectorXYArray( - gVec_c, rMat_d, rMat_s, rMat_c, tVec_d, tVec_s, tVec_c, beamVec=bVec_ref -): - """ - Takes a list of unit reciprocal lattice vectors in crystal frame to the - specified detector-relative frame, subject to the conditions: - - 1) the reciprocal lattice vector must be able to satisfy a bragg condition - 2) the associated diffracted beam must intersect the detector plane - - Required Arguments: - gVec_c -- (n, 3) ndarray of n reciprocal lattice vector components - in the CRYSTAL FRAME - rMat_d -- (3, 3) ndarray, the COB taking DETECTOR FRAME components - to LAB FRAME - rMat_s -- (n, 3, 3) ndarray of n COB taking SAMPLE FRAME components - to LAB FRAME - rMat_c -- (3, 3) ndarray, the COB taking CRYSTAL FRAME components - to SAMPLE FRAME - tVec_d -- (3, 1) ndarray, the translation vector connecting - LAB to DETECTOR in LAB - tVec_s -- (3, 1) ndarray, the translation vector connecting - LAB to SAMPLE in LAB - tVec_c -- (3, 1) ndarray, the translation vector connecting - SAMPLE to CRYSTAL in SAMPLE - - Outputs: - (m, 2) ndarray containing the intersections of m <= n diffracted beams - associated with gVecs - """ - gVec_c = np.ascontiguousarray(gVec_c) - rMat_d = np.ascontiguousarray(rMat_d) - rMat_s = np.ascontiguousarray(rMat_s) - rMat_c = np.ascontiguousarray(rMat_c) - tVec_d = np.ascontiguousarray(tVec_d.flatten()) - tVec_s = np.ascontiguousarray(tVec_s.flatten()) - tVec_c = np.ascontiguousarray(tVec_c.flatten()) - beamVec = np.ascontiguousarray(beamVec.flatten()) - return _transforms_CAPI.gvecToDetectorXYArray( - gVec_c, rMat_d, rMat_s, rMat_c, tVec_d, tVec_s, tVec_c, beamVec - ) - - -def detectorXYToGvec( - xy_det, - rMat_d, - rMat_s, - tVec_d, - tVec_s, - tVec_c, - beamVec=bVec_ref, - etaVec=eta_ref, -): - """ - Takes a list cartesian (x, y) pairs in the detector coordinates and - calculates the associated reciprocal lattice (G) vectors and - (bragg angle, azimuth) pairs with respect to the specified beam and - azimth (eta) reference directions - - Required Arguments: - xy_det -- (n, 2) ndarray or list-like input of n detector (x, y) points - rMat_d -- (3, 3) ndarray, the COB taking DETECTOR FRAME components - to LAB FRAME - rMat_s -- (3, 3) ndarray, the COB taking SAMPLE FRAME components - to LAB FRAME - tVec_d -- (3, 1) ndarray, the translation vector connecting - LAB to DETECTOR in LAB - tVec_s -- (3, 1) ndarray, the translation vector connecting - LAB to SAMPLE in LAB - tVec_c -- (3, 1) ndarray, the translation vector connecting - SAMPLE to CRYSTAL in SAMPLE - - Optional Keyword Arguments: - beamVec -- (3, 1) mdarray containing the incident beam direction - components in the LAB FRAME - etaVec -- (3, 1) mdarray containing the reference azimuth direction - components in the LAB FRAME - - Outputs: - (n, 2) ndarray containing the (tTh, eta) pairs associated with each (x, y) - (n, 3) ndarray containing the associated G vector directions - in the LAB FRAME associated with gVecs - """ - xy_det = np.ascontiguousarray(np.atleast_2d(xy_det)) - rMat_d = np.ascontiguousarray(rMat_d) - rMat_s = np.ascontiguousarray(rMat_s) - tVec_d = np.ascontiguousarray(tVec_d.flatten()) - tVec_s = np.ascontiguousarray(tVec_s.flatten()) - tVec_c = np.ascontiguousarray(tVec_c.flatten()) - beamVec = np.ascontiguousarray(beamVec.flatten()) - etaVec = np.ascontiguousarray(etaVec.flatten()) - return _transforms_CAPI.detectorXYToGvec( - xy_det, rMat_d, rMat_s, tVec_d, tVec_s, tVec_c, beamVec, etaVec - ) - - -def detectorXYToGvecArray( - xy_det, - rMat_d, - rMat_s, - tVec_d, - tVec_s, - tVec_c, - beamVec=bVec_ref, - etaVec=eta_ref, -): - """ - Takes a list cartesian (x, y) pairs in the detector coordinates and - calculates the associated reciprocal lattice (G) vectors and - (bragg angle, azimuth) pairs with respect to the specified beam and azimth - (eta) reference directions - - Required Arguments: - xy_det -- (n, 2) ndarray or list-like input of n detector (x, y) points - rMat_d -- (3, 3) ndarray, the COB taking DETECTOR FRAME components - to LAB FRAME - rMat_s -- (n, 3, 3) ndarray, the COB taking SAMPLE FRAME components - to LAB FRAME - tVec_d -- (3, 1) ndarray, the translation vector connecting - LAB to DETECTOR in LAB - tVec_s -- (3, 1) ndarray, the translation vector connecting - LAB to SAMPLE in LAB - tVec_c -- (3, 1) ndarray, the translation vector connecting - SAMPLE to CRYSTAL in SAMPLE - - Optional Keyword Arguments: - beamVec -- (3, 1) mdarray containing the incident beam direction - components in the LAB FRAME - etaVec -- (3, 1) mdarray containing the reference azimuth direction - components in the LAB FRAME - - Outputs: - (n, 2) ndarray containing the (tTh, eta) pairs associated with each (x, y) - (n, 3) ndarray containing the associated G vector directions - in the LAB FRAME associated with gVecs - """ - xy_det = np.ascontiguousarray(np.atleast_2d(xy_det)) - rMat_d = np.ascontiguousarray(rMat_d) - rMat_s = np.ascontiguousarray(rMat_s) - tVec_d = np.ascontiguousarray(tVec_d.flatten()) - tVec_s = np.ascontiguousarray(tVec_s.flatten()) - tVec_c = np.ascontiguousarray(tVec_c.flatten()) - beamVec = np.ascontiguousarray(beamVec.flatten()) - etaVec = np.ascontiguousarray(etaVec.flatten()) - return _transforms_CAPI.detectorXYToGvec( - xy_det, rMat_d, rMat_s, tVec_d, tVec_s, tVec_c, beamVec, etaVec - ) - - -def oscillAnglesOfHKLs( - hkls, - chi, - rMat_c, - bMat, - wavelength, - vInv=None, - beamVec=bVec_ref, - etaVec=eta_ref, -): - """ - Takes a list of unit reciprocal lattice vectors in crystal frame to the - specified detector-relative frame, subject to the conditions: - - 1) the reciprocal lattice vector must be able to satisfy a bragg condition - 2) the associated diffracted beam must intersect the detector plane - - Required Arguments: - hkls -- (n, 3) ndarray of n reciprocal lattice vectors - in the CRYSTAL FRAME - chi -- float representing the inclination angle of the - oscillation axis (std coords) - rMat_c -- (3, 3) ndarray, the COB taking CRYSTAL FRAME components - to SAMPLE FRAME - bMat -- (3, 3) ndarray, the COB taking RECIPROCAL LATTICE components - to CRYSTAL FRAME - wavelength -- float representing the x-ray wavelength in Angstroms - - Optional Keyword Arguments: - beamVec -- (3, 1) mdarray containing the incident beam direction - components in the LAB FRAME - etaVec -- (3, 1) mdarray containing the reference azimuth direction - components in the LAB FRAME - - Outputs: - ome0 -- (n, 3) ndarray containing the feasible (tTh, eta, ome) triplets - for each input hkl (first solution) - ome1 -- (n, 3) ndarray containing the feasible (tTh, eta, ome) triplets - for each input hkl (second solution) - - Notes: - ------------------------------------------------------------------------ - The reciprocal lattice vector, G, will satisfy the the Bragg condition - when: - - b.T * G / ||G|| = -sin(theta) - - where b is the incident beam direction (k_i) and theta is the Bragg - angle consistent with G and the specified wavelength. The components of - G in the lab frame in this case are obtained using the crystal - orientation, Rc, and the single-parameter oscillation matrix, Rs(ome): - - Rs(ome) * Rc * G / ||G|| - - The equation above can be rearranged to yield an expression of the form: - - a*sin(ome) + b*cos(ome) = c - - which is solved using the relation: - - a*sin(x) + b*cos(x) = sqrt(a**2 + b**2) * sin(x + alpha) - - --> sin(x + alpha) = c / sqrt(a**2 + b**2) - - where: - - alpha = atan2(b, a) - - The solutions are: - - / - | arcsin(c / sqrt(a**2 + b**2)) - alpha - x = < - | pi - arcsin(c / sqrt(a**2 + b**2)) - alpha - \ - - There is a double root in the case the reflection is tangent to the - Debye-Scherrer cone (c**2 = a**2 + b**2), and no solution if the - Laue condition cannot be satisfied (filled with NaNs in the results - array here) - """ - hkls = np.array(hkls, dtype=float, order='C') - if vInv is None: - vInv = np.ascontiguousarray(vInv_ref.flatten()) - else: - vInv = np.ascontiguousarray(vInv.flatten()) - beamVec = np.ascontiguousarray(beamVec.flatten()) - etaVec = np.ascontiguousarray(etaVec.flatten()) - bMat = np.ascontiguousarray(bMat) - return _transforms_CAPI.oscillAnglesOfHKLs( - hkls, chi, rMat_c, bMat, wavelength, vInv, beamVec, etaVec - ) - - -""" -####################################################################### -###### Utility Functions ###### -####################################################################### - -""" - - -def arccosSafe(temp): - """ - Protect against numbers slightly larger than 1 due to round-off - """ - - # Oh, the tricks we must play to make this overloaded and robust... - if type(temp) is list: - temp = np.asarray(temp) - elif type(temp) is np.ndarray: - if len(temp.shape) == 0: - temp = temp.reshape(1) - - if (temp > 1.00001).any() or (temp < -1.00001).any(): - raise RuntimeError(f"Failed to take arccos of {temp}") - - gte1 = temp >= 1.0 - lte1 = temp <= -1.0 - - temp[gte1] = 1 - temp[lte1] = -1 - - return np.arccos(temp) - - -def angularDifference(angList0, angList1, units=angularUnits): - """ - Do the proper (acute) angular difference in the context of a branch cut. - - *) Default angular range is [-pi, pi] - """ - period = periodDict[units] - # take difference as arrays - diffAngles = np.atleast_1d(angList0) - np.atleast_1d(angList1) - - return abs(np.remainder(diffAngles + 0.5 * period, period) - 0.5 * period) - - -def unitRowVector(vecIn): - vecIn = np.ascontiguousarray(vecIn) - if vecIn.ndim == 1: - return _transforms_CAPI.unitRowVector(vecIn) - elif vecIn.ndim == 2: - return _transforms_CAPI.unitRowVectors(vecIn) - else: - assert vecIn.ndim in [ - 1, - 2, - ], "arg shape must be 1-d or 2-d, yours is %d-d" % (vecIn.ndim) - - -def makeDetectorRotMat(tiltAngles): - """ - Form the (3, 3) tilt rotations from the tilt angle list: - - tiltAngles = [gamma_Xl, gamma_Yl, gamma_Zl] in radians - """ - arg = np.ascontiguousarray(np.r_[tiltAngles].flatten()) - return _transforms_CAPI.makeDetectorRotMat(arg) - - -def makeOscillRotMat(oscillAngles): - """ - oscillAngles = [chi, ome] - """ - arg = np.ascontiguousarray(np.r_[oscillAngles].flatten()) - return _transforms_CAPI.makeOscillRotMat(arg) - - -def makeOscillRotMatArray(chi, omeArray): - """ - Applies makeOscillAngles multiple times, for one - chi value and an array of omega values. - """ - arg = np.ascontiguousarray(omeArray) - return _transforms_CAPI.makeOscillRotMatArray(chi, arg) - - -def makeRotMatOfExpMap(expMap): - """ - make a rotation matrix from an exponential map - """ - arg = np.ascontiguousarray(expMap.flatten()) - return _transforms_CAPI.makeRotMatOfExpMap(arg) - - -def makeRotMatOfQuat(quats): - """ - make rotation matrix from vstacked unit quaternions - - """ - arg = np.ascontiguousarray(quats) - return _transforms_CAPI.makeRotMatOfQuat(arg) - - -def makeBinaryRotMat(axis): - arg = np.ascontiguousarray(axis.flatten()) - return _transforms_CAPI.makeBinaryRotMat(arg) - - -def makeEtaFrameRotMat(bHat_l, eHat_l): - arg1 = np.ascontiguousarray(bHat_l.flatten()) - arg2 = np.ascontiguousarray(eHat_l.flatten()) - return _transforms_CAPI.makeEtaFrameRotMat(arg1, arg2) - - -def validateAngleRanges(angList, angMin, angMax, ccw=True): - # FIXME: broken - angList = np.asarray(angList, dtype=np.double, order="C") - angMin = np.asarray(angMin, dtype=np.double, order="C") - angMax = np.asarray(angMax, dtype=np.double, order="C") - return _transforms_CAPI.validateAngleRanges(angList, angMin, angMax, ccw) - - -def rotate_vecs_about_axis(angle, axis, vecs): - return _transforms_CAPI.rotate_vecs_about_axis(angle, axis, vecs) - - -def quat_distance(q1, q2, qsym): - """ - qsym coming from hexrd.crystallogray.PlaneData.getQSym() is C-contiguous - """ - q1 = np.ascontiguousarray(q1.flatten()) - q2 = np.ascontiguousarray(q2.flatten()) - return _transforms_CAPI.quat_distance(q1, q2, qsym) - - -def homochoricOfQuat(quats): - """ - Compute homochoric parameters of unit quaternions - - quats is (4, n) - """ - q = np.ascontiguousarray(quats.T) - return _transforms_CAPI.homochoricOfQuat(q) - - -# def rotateVecsAboutAxis(angle, axis, vecs): -# return _transforms_CAPI.rotateVecsAboutAxis(angle, axis, vecs) diff --git a/hexrd/core/transforms/stdbool.h b/hexrd/core/transforms/stdbool.h deleted file mode 100644 index 3def25f93..000000000 --- a/hexrd/core/transforms/stdbool.h +++ /dev/null @@ -1,6 +0,0 @@ -#pragma once - -#define false 0 -#define true 1 - -#define bool int diff --git a/hexrd/core/transforms/transforms_CAPI.c b/hexrd/core/transforms/transforms_CAPI.c deleted file mode 100644 index 68beb05b1..000000000 --- a/hexrd/core/transforms/transforms_CAPI.c +++ /dev/null @@ -1,1221 +0,0 @@ -/* - * gcc -c -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/include/python2.7 -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/include/numpy transforms_CAPI.c - * - * gcc -bundle -flat_namespace -undefined suppress -o _transforms_CAPI.so transforms_CAPI.o - */ - -#include "transforms_CAPI.h" -#include "transforms_CFUNC.h" - -static PyMethodDef _transform_methods[] = { - {"anglesToGVec",anglesToGVec,METH_VARARGS,"take angle tuples to G-vectors"}, - {"anglesToDVec",anglesToDVec,METH_VARARGS,"take angle tuples to unit diffraction vectors"}, - {"makeGVector",makeGVector,METH_VARARGS,"Make G-vectors from hkls and B-matrix"}, - {"gvecToDetectorXY",gvecToDetectorXY,METH_VARARGS,""}, - {"gvecToDetectorXYArray",gvecToDetectorXYArray,METH_VARARGS,""}, - {"detectorXYToGvec",detectorXYToGvec,METH_VARARGS,"take cartesian coordinates to G-vectors"}, - {"detectorXYToGvecArray",detectorXYToGvecArray,METH_VARARGS,"take cartesian coordinates to G-vectors"}, - {"oscillAnglesOfHKLs",oscillAnglesOfHKLs,METH_VARARGS,"solve angle specs for G-vectors"}, - {"arccosSafe",arccosSafe,METH_VARARGS,""}, - {"angularDifference",angularDifference,METH_VARARGS,"difference for cyclical angles"}, - {"mapAngle",mapAngle,METH_VARARGS,"map cyclical angles to specified range"}, - {"columnNorm",columnNorm,METH_VARARGS,""}, - {"rowNorm",rowNorm,METH_VARARGS,""}, - {"unitRowVector",unitRowVector,METH_VARARGS,"Normalize a single row vector"}, - {"unitRowVectors",unitRowVectors,METH_VARARGS,"Normalize a collection of row vectors"}, - {"makeDetectorRotMat",makeDetectorRotMat,METH_VARARGS,""}, - {"makeOscillRotMat",makeOscillRotMat,METH_VARARGS,""}, - {"makeOscillRotMatArray",makeOscillRotMatArray,METH_VARARGS,""}, - {"makeRotMatOfExpMap",makeRotMatOfExpMap,METH_VARARGS,""}, - {"makeRotMatOfQuat",makeRotMatOfQuat,METH_VARARGS,""}, - {"makeBinaryRotMat",makeBinaryRotMat,METH_VARARGS,""}, - {"makeEtaFrameRotMat",makeEtaFrameRotMat,METH_VARARGS,"Make eta basis COB matrix"}, - {"validateAngleRanges",validateAngleRanges,METH_VARARGS,""}, - {"rotate_vecs_about_axis",rotate_vecs_about_axis,METH_VARARGS,"Rotate vectors about an axis"}, - {"quat_distance",quat_distance,METH_VARARGS,"Compute distance between two unit quaternions"}, - {"homochoricOfQuat",homochoricOfQuat,METH_VARARGS,"Compute homochoric parameterization of list of unit quaternions"}, - {NULL,NULL} -}; - -static struct PyModuleDef transforms_capi_pymodule = -{ - PyModuleDef_HEAD_INIT, - "_transforms_CAPI", /* name of module */ - "", /* module documentation, may be NULL */ - -1, /* size of per-interpreter state of the module, or -1 if the module keeps state in global variables. */ - _transform_methods -}; - -PyMODINIT_FUNC PyInit__transforms_CAPI(void) -{ - import_array(); - return PyModule_Create(&transforms_capi_pymodule); -} - -/******************************************************************************/ -/* Funtions */ - -static PyObject * anglesToGVec(PyObject * self, PyObject * args) -{ - PyArrayObject *angs, *bHat_l, *eHat_l, *rMat_c; - PyArrayObject *gVec_c; - double chi; - npy_intp nvecs, rdims[2]; - - int nangs, nbhat, nehat, nrmat; - int da1, db1, de1, dr1, dr2; - - double *angs_ptr, *bHat_l_ptr, *eHat_l_ptr, *rMat_c_ptr; - double *gVec_c_ptr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OOOdO", - &angs, - &bHat_l, &eHat_l, - &chi, &rMat_c)) return(NULL); - if ( angs == NULL ) return(NULL); - - /* Verify shape of input arrays */ - nangs = PyArray_NDIM(angs); - nbhat = PyArray_NDIM(bHat_l); - nehat = PyArray_NDIM(eHat_l); - nrmat = PyArray_NDIM(rMat_c); - - assert( nangs==2 && nbhat==1 && nehat==1 && nrmat==2 ); - - /* Verify dimensions of input arrays */ - nvecs = PyArray_DIMS(angs)[0]; //rows - da1 = PyArray_DIMS(angs)[1]; //cols - - db1 = PyArray_DIMS(bHat_l)[0]; - de1 = PyArray_DIMS(eHat_l)[0]; - dr1 = PyArray_DIMS(rMat_c)[0]; - dr2 = PyArray_DIMS(rMat_c)[1]; - - assert( da1 == 3 ); - assert( db1 == 3 && de1 == 3); - assert( dr1 == 3 && dr2 == 3); - - /* Allocate C-style array for return data */ - rdims[0] = nvecs; rdims[1] = 3; - gVec_c = (PyArrayObject*)PyArray_EMPTY(2,rdims,NPY_DOUBLE,0); - - /* Grab pointers to the various data arrays */ - angs_ptr = (double*)PyArray_DATA(angs); - bHat_l_ptr = (double*)PyArray_DATA(bHat_l); - eHat_l_ptr = (double*)PyArray_DATA(eHat_l); - rMat_c_ptr = (double*)PyArray_DATA(rMat_c); - gVec_c_ptr = (double*)PyArray_DATA(gVec_c); - - /* Call the actual function */ - anglesToGvec_cfunc(nvecs, angs_ptr, - bHat_l_ptr, eHat_l_ptr, - chi, rMat_c_ptr, - gVec_c_ptr); - - /* Build and return the nested data structure */ - return((PyObject*)gVec_c); -} - -static PyObject * anglesToDVec(PyObject * self, PyObject * args) -{ - PyArrayObject *angs, *bHat_l, *eHat_l, *rMat_c; - PyArrayObject *dVec_c; - double chi; - npy_intp nvecs, rdims[2]; - - int nangs, nbhat, nehat, nrmat; - int da1, db1, de1, dr1, dr2; - - double *angs_ptr, *bHat_l_ptr, *eHat_l_ptr, *rMat_c_ptr; - double *dVec_c_ptr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OOOdO", - &angs, - &bHat_l, &eHat_l, - &chi, &rMat_c)) return(NULL); - if ( angs == NULL ) return(NULL); - - /* Verify shape of input arrays */ - nangs = PyArray_NDIM(angs); - nbhat = PyArray_NDIM(bHat_l); - nehat = PyArray_NDIM(eHat_l); - nrmat = PyArray_NDIM(rMat_c); - - assert( nangs==2 && nbhat==1 && nehat==1 && nrmat==2 ); - - /* Verify dimensions of input arrays */ - nvecs = PyArray_DIMS(angs)[0]; //rows - da1 = PyArray_DIMS(angs)[1]; //cols - - db1 = PyArray_DIMS(bHat_l)[0]; - de1 = PyArray_DIMS(eHat_l)[0]; - dr1 = PyArray_DIMS(rMat_c)[0]; - dr2 = PyArray_DIMS(rMat_c)[1]; - - assert( da1 == 3 ); - assert( db1 == 3 && de1 == 3); - assert( dr1 == 3 && dr2 == 3); - - /* Allocate C-style array for return data */ - rdims[0] = nvecs; rdims[1] = 3; - dVec_c = (PyArrayObject*)PyArray_EMPTY(2,rdims,NPY_DOUBLE,0); - - /* Grab pointers to the various data arrays */ - angs_ptr = (double*)PyArray_DATA(angs); - bHat_l_ptr = (double*)PyArray_DATA(bHat_l); - eHat_l_ptr = (double*)PyArray_DATA(eHat_l); - rMat_c_ptr = (double*)PyArray_DATA(rMat_c); - dVec_c_ptr = (double*)PyArray_DATA(dVec_c); - - /* Call the actual function */ - anglesToDvec_cfunc(nvecs, angs_ptr, - bHat_l_ptr, eHat_l_ptr, - chi, rMat_c_ptr, - dVec_c_ptr); - - /* Build and return the nested data structure */ - return((PyObject*)dVec_c); -} - -static PyObject * makeGVector(PyObject * self, PyObject * args) -{ - return(NULL); -} - -/* - Takes a list of unit reciprocal lattice vectors in crystal frame to the - specified detector-relative frame, subject to the conditions: - - 1) the reciprocal lattice vector must be able to satisfy a bragg condition - 2) the associated diffracted beam must intersect the detector plane - - Required Arguments: - gVec_c -- (n, 3) ndarray of n reciprocal lattice vectors in the CRYSTAL FRAME - rMat_d -- (3, 3) ndarray, the COB taking DETECTOR FRAME components to LAB FRAME - rMat_s -- (3, 3) ndarray, the COB taking SAMPLE FRAME components to LAB FRAME - rMat_c -- (3, 3) ndarray, the COB taking CRYSTAL FRAME components to SAMPLE FRAME - tVec_d -- (3, 1) ndarray, the translation vector connecting LAB to DETECTOR - tVec_s -- (3, 1) ndarray, the translation vector connecting LAB to SAMPLE - tVec_c -- (3, 1) ndarray, the translation vector connecting SAMPLE to CRYSTAL - - Outputs: - (m, 2) ndarray containing the intersections of m <= n diffracted beams - associated with gVecs -*/ -static PyObject * gvecToDetectorXY(PyObject * self, PyObject * args) -{ - PyArrayObject *gVec_c, - *rMat_d, *rMat_s, *rMat_c, - *tVec_d, *tVec_s, *tVec_c, - *beamVec; - PyArrayObject *result; - - int dgc, drd, drs, drc, dtd, dts, dtc, dbv; - npy_intp npts, dims[2]; - - double *gVec_c_Ptr, - *rMat_d_Ptr, *rMat_s_Ptr, *rMat_c_Ptr, - *tVec_d_Ptr, *tVec_s_Ptr, *tVec_c_Ptr, - *beamVec_Ptr; - double *result_Ptr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OOOOOOOO", - &gVec_c, - &rMat_d, &rMat_s, &rMat_c, - &tVec_d, &tVec_s, &tVec_c, - &beamVec)) return(NULL); - if ( gVec_c == NULL || - rMat_d == NULL || rMat_s == NULL || rMat_c == NULL || - tVec_d == NULL || tVec_s == NULL || tVec_c == NULL || - beamVec == NULL ) return(NULL); - - /* Verify shape of input arrays */ - dgc = PyArray_NDIM(gVec_c); - drd = PyArray_NDIM(rMat_d); - drs = PyArray_NDIM(rMat_s); - drc = PyArray_NDIM(rMat_c); - dtd = PyArray_NDIM(tVec_d); - dts = PyArray_NDIM(tVec_s); - dtc = PyArray_NDIM(tVec_c); - dbv = PyArray_NDIM(beamVec); - assert( dgc == 2 ); - assert( drd == 2 && drs == 2 && drc == 2 ); - assert( dtd == 1 && dts == 1 && dtc == 1 ); - assert( dbv == 1 ); - - /* Verify dimensions of input arrays */ - npts = PyArray_DIMS(gVec_c)[0]; - - assert( PyArray_DIMS(gVec_c)[1] == 3 ); - assert( PyArray_DIMS(rMat_d)[0] == 3 && PyArray_DIMS(rMat_d)[1] == 3 ); - assert( PyArray_DIMS(rMat_s)[0] == 3 && PyArray_DIMS(rMat_s)[1] == 3 ); - assert( PyArray_DIMS(rMat_c)[0] == 3 && PyArray_DIMS(rMat_c)[1] == 3 ); - assert( PyArray_DIMS(tVec_d)[0] == 3 ); - assert( PyArray_DIMS(tVec_s)[0] == 3 ); - assert( PyArray_DIMS(tVec_c)[0] == 3 ); - assert( PyArray_DIMS(beamVec)[0] == 3 ); - - /* Allocate C-style array for return data */ - dims[0] = npts; dims[1] = 2; - result = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - - /* Grab data pointers into various arrays */ - gVec_c_Ptr = (double*)PyArray_DATA(gVec_c); - - rMat_d_Ptr = (double*)PyArray_DATA(rMat_d); - rMat_s_Ptr = (double*)PyArray_DATA(rMat_s); - rMat_c_Ptr = (double*)PyArray_DATA(rMat_c); - - tVec_d_Ptr = (double*)PyArray_DATA(tVec_d); - tVec_s_Ptr = (double*)PyArray_DATA(tVec_s); - tVec_c_Ptr = (double*)PyArray_DATA(tVec_c); - - beamVec_Ptr = (double*)PyArray_DATA(beamVec); - - result_Ptr = (double*)PyArray_DATA(result); - - /* Call the computational routine */ - gvecToDetectorXY_cfunc(npts, gVec_c_Ptr, - rMat_d_Ptr, rMat_s_Ptr, rMat_c_Ptr, - tVec_d_Ptr, tVec_s_Ptr, tVec_c_Ptr, - beamVec_Ptr, - result_Ptr); - - /* Build and return the nested data structure */ - return((PyObject*)result); -} - -/* - Takes a list of unit reciprocal lattice vectors in crystal frame to the - specified detector-relative frame, subject to the conditions: - - 1) the reciprocal lattice vector must be able to satisfy a bragg condition - 2) the associated diffracted beam must intersect the detector plane - - Required Arguments: - gVec_c -- (n, 3) ndarray of n reciprocal lattice vectors in the CRYSTAL FRAME - rMat_d -- (3, 3) ndarray, the COB taking DETECTOR FRAME components to LAB FRAME - rMat_s -- (n, 3, 3) ndarray, the COB taking SAMPLE FRAME components to LAB FRAME - rMat_c -- (3, 3) ndarray, the COB taking CRYSTAL FRAME components to SAMPLE FRAME - tVec_d -- (3, 1) ndarray, the translation vector connecting LAB to DETECTOR - tVec_s -- (3, 1) ndarray, the translation vector connecting LAB to SAMPLE - tVec_c -- (3, 1) ndarray, the translation vector connecting SAMPLE to CRYSTAL - - Outputs: - (m, 2) ndarray containing the intersections of m <= n diffracted beams - associated with gVecs -*/ -static PyObject * gvecToDetectorXYArray(PyObject * self, PyObject * args) -{ - PyArrayObject *gVec_c, - *rMat_d, *rMat_s, *rMat_c, - *tVec_d, *tVec_s, *tVec_c, - *beamVec; - PyArrayObject *result; - - int dgc, drd, drs, drc, dtd, dts, dtc, dbv; - npy_intp npts, dims[2]; - - double *gVec_c_Ptr, - *rMat_d_Ptr, *rMat_s_Ptr, *rMat_c_Ptr, - *tVec_d_Ptr, *tVec_s_Ptr, *tVec_c_Ptr, - *beamVec_Ptr; - double *result_Ptr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OOOOOOOO", - &gVec_c, - &rMat_d, &rMat_s, &rMat_c, - &tVec_d, &tVec_s, &tVec_c, - &beamVec)) return(NULL); - if ( gVec_c == NULL || - rMat_d == NULL || rMat_s == NULL || rMat_c == NULL || - tVec_d == NULL || tVec_s == NULL || tVec_c == NULL || - beamVec == NULL ) return(NULL); - - /* Verify shape of input arrays */ - dgc = PyArray_NDIM(gVec_c); - drd = PyArray_NDIM(rMat_d); - drs = PyArray_NDIM(rMat_s); - drc = PyArray_NDIM(rMat_c); - dtd = PyArray_NDIM(tVec_d); - dts = PyArray_NDIM(tVec_s); - dtc = PyArray_NDIM(tVec_c); - dbv = PyArray_NDIM(beamVec); - assert( dgc == 2 ); - assert( drd == 2 && drs == 3 && drc == 2 ); - assert( dtd == 1 && dts == 1 && dtc == 1 ); - assert( dbv == 1 ); - - /* Verify dimensions of input arrays */ - npts = PyArray_DIMS(gVec_c)[0]; - - if (npts != PyArray_DIM(rMat_s, 0)) { - PyErr_Format(PyExc_ValueError, "gVec_c and rMat_s length mismatch %d vs %d", - (int)PyArray_DIM(gVec_c, 0), (int)PyArray_DIM(rMat_s, 0)); - return NULL; - } - assert( PyArray_DIMS(gVec_c)[1] == 3 ); - assert( PyArray_DIMS(rMat_d)[0] == 3 && PyArray_DIMS(rMat_d)[1] == 3 ); - assert( PyArray_DIMS(rMat_s)[1] == 3 && PyArray_DIMS(rMat_s)[2] == 3 ); - assert( PyArray_DIMS(rMat_c)[0] == 3 && PyArray_DIMS(rMat_c)[1] == 3 ); - assert( PyArray_DIMS(tVec_d)[0] == 3 ); - assert( PyArray_DIMS(tVec_s)[0] == 3 ); - assert( PyArray_DIMS(tVec_c)[0] == 3 ); - assert( PyArray_DIMS(beamVec)[0] == 3 ); - - /* Allocate C-style array for return data */ - dims[0] = npts; dims[1] = 2; - result = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - - /* Grab data pointers into various arrays */ - gVec_c_Ptr = (double*)PyArray_DATA(gVec_c); - - rMat_d_Ptr = (double*)PyArray_DATA(rMat_d); - rMat_s_Ptr = (double*)PyArray_DATA(rMat_s); - rMat_c_Ptr = (double*)PyArray_DATA(rMat_c); - - tVec_d_Ptr = (double*)PyArray_DATA(tVec_d); - tVec_s_Ptr = (double*)PyArray_DATA(tVec_s); - tVec_c_Ptr = (double*)PyArray_DATA(tVec_c); - - beamVec_Ptr = (double*)PyArray_DATA(beamVec); - - result_Ptr = (double*)PyArray_DATA(result); - - /* Call the computational routine */ - gvecToDetectorXYArray_cfunc(npts, gVec_c_Ptr, - rMat_d_Ptr, rMat_s_Ptr, rMat_c_Ptr, - tVec_d_Ptr, tVec_s_Ptr, tVec_c_Ptr, - beamVec_Ptr, - result_Ptr); - - /* Build and return the nested data structure */ - return((PyObject*)result); -} - -/* - Takes a list cartesian (x, y) pairs in the detector coordinates and calculates - the associated reciprocal lattice (G) vectors and (bragg angle, azimuth) pairs - with respect to the specified beam and azimth (eta) reference directions - - Required Arguments: - xy_det -- (n, 2) ndarray or list-like input of n detector (x, y) points - rMat_d -- (3, 3) ndarray, the COB taking DETECTOR FRAME components to LAB FRAME - rMat_s -- (3, 3) ndarray, the COB taking SAMPLE FRAME components to LAB FRAME - tVec_d -- (3, 1) ndarray, the translation vector connecting LAB to DETECTOR - tVec_s -- (3, 1) ndarray, the translation vector connecting LAB to SAMPLE - tVec_c -- (3, 1) ndarray, the translation vector connecting SAMPLE to CRYSTAL - - Optional Keyword Arguments: - beamVec -- (1, 3) mdarray containing the incident beam direction components in the LAB FRAME - etaVec -- (1, 3) mdarray containing the reference azimuth direction components in the LAB FRAME - - Outputs: - (n, 2) ndarray containing the (tTh, eta) pairs associated with each (x, y) - (n, 3) ndarray containing the associated G vector directions in the LAB FRAME - associated with gVecs -*/ -static PyObject * detectorXYToGvec(PyObject * self, PyObject * args) -{ - PyArrayObject *xy_det, *rMat_d, *rMat_s, - *tVec_d, *tVec_s, *tVec_c, - *beamVec, *etaVec; - PyArrayObject *tTh, *eta, *gVec_l; - PyObject *inner_tuple, *outer_tuple; - - int dxy, drd, drs, dtd, dts, dtc, dbv, dev; - npy_intp npts, dims[2]; - - double *xy_Ptr, *rMat_d_Ptr, *rMat_s_Ptr, - *tVec_d_Ptr, *tVec_s_Ptr, *tVec_c_Ptr, - *beamVec_Ptr, *etaVec_Ptr; - double *tTh_Ptr, *eta_Ptr, *gVec_l_Ptr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OOOOOOOO", - &xy_det, - &rMat_d, &rMat_s, - &tVec_d, &tVec_s, &tVec_c, - &beamVec, &etaVec)) return(NULL); - if ( xy_det == NULL || rMat_d == NULL || rMat_s == NULL || - tVec_d == NULL || tVec_s == NULL || tVec_c == NULL || - beamVec == NULL || etaVec == NULL ) return(NULL); - - /* Verify shape of input arrays */ - dxy = PyArray_NDIM(xy_det); - drd = PyArray_NDIM(rMat_d); - drs = PyArray_NDIM(rMat_s); - dtd = PyArray_NDIM(tVec_d); - dts = PyArray_NDIM(tVec_s); - dtc = PyArray_NDIM(tVec_c); - dbv = PyArray_NDIM(beamVec); - dev = PyArray_NDIM(etaVec); - assert( dxy == 2 && drd == 2 && drs == 2 && - dtd == 1 && dts == 1 && dtc == 1 && - dbv == 1 && dev == 1); - - /* Verify dimensions of input arrays */ - npts = PyArray_DIMS(xy_det)[0]; - - assert( PyArray_DIMS(xy_det)[1] == 2 ); - assert( PyArray_DIMS(rMat_d)[0] == 3 && PyArray_DIMS(rMat_d)[1] == 3 ); - assert( PyArray_DIMS(rMat_s)[0] == 3 && PyArray_DIMS(rMat_s)[1] == 3 ); - assert( PyArray_DIMS(tVec_d)[0] == 3 ); - assert( PyArray_DIMS(tVec_s)[0] == 3 ); - assert( PyArray_DIMS(tVec_c)[0] == 3 ); - assert( PyArray_DIMS(beamVec)[0] == 3 ); - assert( PyArray_DIMS(etaVec)[0] == 3 ); - - /* Allocate arrays for return values */ - dims[0] = npts; dims[1] = 3; - gVec_l = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - - tTh = (PyArrayObject*)PyArray_EMPTY(1,&npts,NPY_DOUBLE,0); - eta = (PyArrayObject*)PyArray_EMPTY(1,&npts,NPY_DOUBLE,0); - - /* Grab data pointers into various arrays */ - xy_Ptr = (double*)PyArray_DATA(xy_det); - gVec_l_Ptr = (double*)PyArray_DATA(gVec_l); - - tTh_Ptr = (double*)PyArray_DATA(tTh); - eta_Ptr = (double*)PyArray_DATA(eta); - - rMat_d_Ptr = (double*)PyArray_DATA(rMat_d); - rMat_s_Ptr = (double*)PyArray_DATA(rMat_s); - - tVec_d_Ptr = (double*)PyArray_DATA(tVec_d); - tVec_s_Ptr = (double*)PyArray_DATA(tVec_s); - tVec_c_Ptr = (double*)PyArray_DATA(tVec_c); - - beamVec_Ptr = (double*)PyArray_DATA(beamVec); - etaVec_Ptr = (double*)PyArray_DATA(etaVec); - - /* Call the computational routine */ - detectorXYToGvec_cfunc(npts, xy_Ptr, - rMat_d_Ptr, rMat_s_Ptr, - tVec_d_Ptr, tVec_s_Ptr, tVec_c_Ptr, - beamVec_Ptr, etaVec_Ptr, - tTh_Ptr, eta_Ptr, gVec_l_Ptr); - - /* Build and return the nested data structure */ - /* Note that Py_BuildValue with 'O' increases reference count */ - inner_tuple = Py_BuildValue("OO",tTh,eta); - outer_tuple = Py_BuildValue("OO", inner_tuple, gVec_l); - Py_DECREF(inner_tuple); - Py_DECREF(tTh); - Py_DECREF(eta); - Py_DECREF(gVec_l); - return outer_tuple; -} - -/* - Takes a list cartesian (x, y) pairs in the detector coordinates and calculates - the associated reciprocal lattice (G) vectors and (bragg angle, azimuth) pairs - with respect to the specified beam and azimth (eta) reference directions - - Required Arguments: - xy_det -- (n, 2) ndarray or list-like input of n detector (x, y) points - rMat_d -- (3, 3) ndarray, the COB taking DETECTOR FRAME components to LAB FRAME - rMat_s -- (n, 3, 3) ndarray, the COB taking SAMPLE FRAME components to LAB FRAME - tVec_d -- (3, 1) ndarray, the translation vector connecting LAB to DETECTOR - tVec_s -- (3, 1) ndarray, the translation vector connecting LAB to SAMPLE - tVec_c -- (3, 1) ndarray, the translation vector connecting SAMPLE to CRYSTAL - - Optional Keyword Arguments: - beamVec -- (1, 3) mdarray containing the incident beam direction components in the LAB FRAME - etaVec -- (1, 3) mdarray containing the reference azimuth direction components in the LAB FRAME - - Outputs: - (n, 2) ndarray containing the (tTh, eta) pairs associated with each (x, y) - (n, 3) ndarray containing the associated G vector directions in the LAB FRAME - associated with gVecs -*/ -static PyObject * detectorXYToGvecArray(PyObject * self, PyObject * args) -{ - PyArrayObject *xy_det, *rMat_d, *rMat_s, - *tVec_d, *tVec_s, *tVec_c, - *beamVec, *etaVec; - PyArrayObject *tTh, *eta, *gVec_l; - PyObject *inner_tuple, *outer_tuple; - - int dxy, drd, drs, dtd, dts, dtc, dbv, dev; - npy_intp npts, dims[2]; - - double *xy_Ptr, *rMat_d_Ptr, *rMat_s_Ptr, - *tVec_d_Ptr, *tVec_s_Ptr, *tVec_c_Ptr, - *beamVec_Ptr, *etaVec_Ptr; - double *tTh_Ptr, *eta_Ptr, *gVec_l_Ptr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OOOOOOOO", - &xy_det, - &rMat_d, &rMat_s, - &tVec_d, &tVec_s, &tVec_c, - &beamVec, &etaVec)) return(NULL); - if ( xy_det == NULL || rMat_d == NULL || rMat_s == NULL || - tVec_d == NULL || tVec_s == NULL || tVec_c == NULL || - beamVec == NULL || etaVec == NULL ) return(NULL); - - /* Verify shape of input arrays */ - dxy = PyArray_NDIM(xy_det); - drd = PyArray_NDIM(rMat_d); - drs = PyArray_NDIM(rMat_s); - dtd = PyArray_NDIM(tVec_d); - dts = PyArray_NDIM(tVec_s); - dtc = PyArray_NDIM(tVec_c); - dbv = PyArray_NDIM(beamVec); - dev = PyArray_NDIM(etaVec); - assert( dxy == 2 && drd == 2 && drs == 2 && - dtd == 1 && dts == 1 && dtc == 1 && - dbv == 1 && dev == 1); - - /* Verify dimensions of input arrays */ - npts = PyArray_DIMS(xy_det)[0]; - if (npts != PyArray_DIM(rMat_s, 0)) { - PyErr_Format(PyExc_ValueError, "xy_det and rMat_s length mismatch %d vs %d", - (int)PyArray_DIM(xy_det, 0), (int)PyArray_DIM(rMat_s, 0)); - return NULL; - } - - assert( PyArray_DIMS(xy_det)[1] == 2 ); - assert( PyArray_DIMS(rMat_d)[0] == 3 && PyArray_DIMS(rMat_d)[1] == 3 ); - assert( PyArray_DIMS(rMat_s)[0] == 3 && PyArray_DIMS(rMat_s)[1] == 3 ); - assert( PyArray_DIMS(tVec_d)[0] == 3 ); - assert( PyArray_DIMS(tVec_s)[0] == 3 ); - assert( PyArray_DIMS(tVec_c)[0] == 3 ); - assert( PyArray_DIMS(beamVec)[0] == 3 ); - assert( PyArray_DIMS(etaVec)[0] == 3 ); - - /* Allocate arrays for return values */ - dims[0] = npts; dims[1] = 3; - gVec_l = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - - tTh = (PyArrayObject*)PyArray_EMPTY(1,&npts,NPY_DOUBLE,0); - eta = (PyArrayObject*)PyArray_EMPTY(1,&npts,NPY_DOUBLE,0); - - /* Grab data pointers into various arrays */ - xy_Ptr = (double*)PyArray_DATA(xy_det); - gVec_l_Ptr = (double*)PyArray_DATA(gVec_l); - - tTh_Ptr = (double*)PyArray_DATA(tTh); - eta_Ptr = (double*)PyArray_DATA(eta); - - rMat_d_Ptr = (double*)PyArray_DATA(rMat_d); - rMat_s_Ptr = (double*)PyArray_DATA(rMat_s); - - tVec_d_Ptr = (double*)PyArray_DATA(tVec_d); - tVec_s_Ptr = (double*)PyArray_DATA(tVec_s); - tVec_c_Ptr = (double*)PyArray_DATA(tVec_c); - - beamVec_Ptr = (double*)PyArray_DATA(beamVec); - etaVec_Ptr = (double*)PyArray_DATA(etaVec); - - /* Call the computational routine */ - detectorXYToGvecArray_cfunc(npts, xy_Ptr, - rMat_d_Ptr, rMat_s_Ptr, - tVec_d_Ptr, tVec_s_Ptr, tVec_c_Ptr, - beamVec_Ptr, etaVec_Ptr, - tTh_Ptr, eta_Ptr, gVec_l_Ptr); - - /* Build and return the nested data structure */ - /* Note that Py_BuildValue with 'O' increases reference count */ - inner_tuple = Py_BuildValue("OO",tTh,eta); - outer_tuple = Py_BuildValue("OO", inner_tuple, gVec_l); - Py_DECREF(inner_tuple); - Py_DECREF(tTh); - Py_DECREF(eta); - Py_DECREF(gVec_l); - return outer_tuple; -} - -static PyObject * oscillAnglesOfHKLs(PyObject * self, PyObject * args) -{ - PyArrayObject *hkls, *rMat_c, *bMat, - *vInv_s, *beamVec, *etaVec; - PyFloatObject *chi, *wavelength; - PyArrayObject *oangs0, *oangs1; - PyObject *return_tuple; - - int dhkls, drc, dbm, dvi, dbv, dev; - npy_intp npts, dims[2]; - - double *hkls_Ptr, chi_d, - *rMat_c_Ptr, *bMat_Ptr, wavelen_d, - *vInv_s_Ptr, *beamVec_Ptr, *etaVec_Ptr; - double *oangs0_Ptr, *oangs1_Ptr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OOOOOOOO", - &hkls, &chi, - &rMat_c, &bMat, &wavelength, - &vInv_s, &beamVec, &etaVec)) return(NULL); - if ( hkls == NULL || chi == NULL || - rMat_c == NULL || bMat == NULL || wavelength == NULL || - vInv_s == NULL || beamVec == NULL || etaVec == NULL ) return(NULL); - - /* Verify shape of input arrays */ - dhkls = PyArray_NDIM(hkls); - drc = PyArray_NDIM(rMat_c); - dbm = PyArray_NDIM(bMat); - dvi = PyArray_NDIM(vInv_s); - dbv = PyArray_NDIM(beamVec); - dev = PyArray_NDIM(etaVec); - assert( dhkls == 2 && drc == 2 && dbm == 2 && - dvi == 1 && dbv == 1 && dev == 1); - - /* Verify dimensions of input arrays */ - npts = PyArray_DIMS(hkls)[0]; - - assert( PyArray_DIMS(hkls)[1] == 3 ); - assert( PyArray_DIMS(rMat_c)[0] == 3 && PyArray_DIMS(rMat_c)[1] == 3 ); - assert( PyArray_DIMS(bMat)[0] == 3 && PyArray_DIMS(bMat)[1] == 3 ); - assert( PyArray_DIMS(vInv_s)[0] == 6 ); - assert( PyArray_DIMS(beamVec)[0] == 3 ); - assert( PyArray_DIMS(etaVec)[0] == 3 ); - - /* Allocate arrays for return values */ - dims[0] = npts; dims[1] = 3; - oangs0 = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - oangs1 = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - - /* Grab data pointers into various arrays */ - hkls_Ptr = (double*)PyArray_DATA(hkls); - - chi_d = PyFloat_AsDouble((PyObject*)chi); - wavelen_d = PyFloat_AsDouble((PyObject*)wavelength); - - rMat_c_Ptr = (double*)PyArray_DATA(rMat_c); - bMat_Ptr = (double*)PyArray_DATA(bMat); - - vInv_s_Ptr = (double*)PyArray_DATA(vInv_s); - - beamVec_Ptr = (double*)PyArray_DATA(beamVec); - etaVec_Ptr = (double*)PyArray_DATA(etaVec); - - oangs0_Ptr = (double*)PyArray_DATA(oangs0); - oangs1_Ptr = (double*)PyArray_DATA(oangs1); - - /* Call the computational routine */ - oscillAnglesOfHKLs_cfunc(npts, hkls_Ptr, chi_d, - rMat_c_Ptr, bMat_Ptr, wavelen_d, - vInv_s_Ptr, beamVec_Ptr, etaVec_Ptr, - oangs0_Ptr, oangs1_Ptr); - - /* Build and return the list data structure */ - return_tuple = Py_BuildValue("OO",oangs0,oangs1); - Py_DECREF(oangs1); - Py_DECREF(oangs0); - - return return_tuple; -} - -/******************************************************************************/ -/* Utility Funtions */ - -static PyObject * arccosSafe(PyObject * self, PyObject * args) -{ - return(NULL); -} - -static PyObject * angularDifference(PyObject * self, PyObject * args) -{ - return(NULL); -} - -static PyObject * mapAngle(PyObject * self, PyObject * args) -{ - return(NULL); -} - -static PyObject * columnNorm(PyObject * self, PyObject * args) -{ - return(NULL); -} - -static PyObject * rowNorm(PyObject * self, PyObject * args) -{ - return(NULL); -} - -static PyObject * unitRowVector(PyObject * self, PyObject * args) -{ - PyArrayObject *vecIn, *vecOut; - double *cIn, *cOut; - int d; - npy_intp n; - - if ( !PyArg_ParseTuple(args,"O", &vecIn)) return(NULL); - if ( vecIn == NULL ) return(NULL); - - assert( PyArray_ISCONTIGUOUS(vecIn) ); - assert( PyArray_ISALIGNED(vecIn) ); - - d = PyArray_NDIM(vecIn); - - assert(d == 1); - - n = PyArray_DIMS(vecIn)[0]; - - vecOut = (PyArrayObject*)PyArray_EMPTY(d,PyArray_DIMS(vecIn),NPY_DOUBLE,0); - - cIn = (double*)PyArray_DATA(vecIn); - cOut = (double*)PyArray_DATA(vecOut); - - unitRowVector_cfunc(n,cIn,cOut); - - return((PyObject*)vecOut); -} - -static PyObject * unitRowVectors(PyObject *self, PyObject *args) -{ - PyArrayObject *vecIn, *vecOut; - double *cIn, *cOut; - int d; - npy_intp m,n; - - if ( !PyArg_ParseTuple(args,"O", &vecIn)) return(NULL); - if ( vecIn == NULL ) return(NULL); - - assert( PyArray_ISCONTIGUOUS(vecIn) ); - assert( PyArray_ISALIGNED(vecIn) ); - - d = PyArray_NDIM(vecIn); - - assert(d == 2); - - m = PyArray_DIMS(vecIn)[0]; - n = PyArray_DIMS(vecIn)[1]; - - vecOut = (PyArrayObject*)PyArray_EMPTY(d,PyArray_DIMS(vecIn),NPY_DOUBLE,0); - - cIn = (double*)PyArray_DATA(vecIn); - cOut = (double*)PyArray_DATA(vecOut); - - unitRowVectors_cfunc(m,n,cIn,cOut); - - return((PyObject*)vecOut); -} - -static PyObject * makeDetectorRotMat(PyObject * self, PyObject * args) -{ - PyArrayObject *tiltAngles, *rMat; - int dt; - npy_intp nt, dims[2]; - double *tPtr, *rPtr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"O", &tiltAngles)) return(NULL); - if ( tiltAngles == NULL ) return(NULL); - - /* Verify shape of input arrays */ - dt = PyArray_NDIM(tiltAngles); - assert( dt == 1 ); - - /* Verify dimensions of input arrays */ - nt = PyArray_DIMS(tiltAngles)[0]; - assert( nt == 3 ); - - /* Allocate the result matrix with appropriate dimensions and type */ - dims[0] = 3; dims[1] = 3; - rMat = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - - /* Grab pointers to the various data arrays */ - tPtr = (double*)PyArray_DATA(tiltAngles); - rPtr = (double*)PyArray_DATA(rMat); - - /* Call the actual function */ - makeDetectorRotMat_cfunc(tPtr,rPtr); - - return((PyObject*)rMat); -} - -static PyObject * makeOscillRotMat(PyObject * self, PyObject * args) -{ - PyArrayObject *oscillAngles, *rMat; - int doa; - npy_intp no, dims[2]; - double *oPtr, *rPtr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"O", &oscillAngles)) return(NULL); - if ( oscillAngles == NULL ) return(NULL); - - /* Verify shape of input arrays */ - doa = PyArray_NDIM(oscillAngles); - assert( doa == 1 ); - - /* Verify dimensions of input arrays */ - no = PyArray_DIMS(oscillAngles)[0]; - assert( no == 2 ); - - /* Allocate the result matrix with appropriate dimensions and type */ - dims[0] = 3; dims[1] = 3; - rMat = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - - /* Grab pointers to the various data arrays */ - oPtr = (double*)PyArray_DATA(oscillAngles); - rPtr = (double*)PyArray_DATA(rMat); - - /* Call the actual function */ - makeOscillRotMat_cfunc(oPtr[0], oPtr[1], rPtr); - - return((PyObject*)rMat); -} - -static PyObject * makeOscillRotMatArray(PyObject * self, PyObject * args) -{ - PyObject *chiObj; - double chi; - PyArrayObject *omeArray, *rMat; - int doa; - npy_intp i, no, dims[3]; - double *oPtr, *rPtr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OO", &chiObj, &omeArray)) return(NULL); - if ( chiObj == NULL ) return(NULL); - if ( omeArray == NULL ) return(NULL); - - /* Get chi */ - chi = PyFloat_AsDouble(chiObj); - if (chi == -1 && PyErr_Occurred()) return(NULL); - - /* Verify shape of input arrays */ - doa = PyArray_NDIM(omeArray); - assert( doa == 1 ); - - /* Verify dimensions of input arrays */ - no = PyArray_DIMS(omeArray)[0]; - - /* Allocate the result matrix with appropriate dimensions and type */ - dims[0] = no; dims[1] = 3; dims[2] = 3; - rMat = (PyArrayObject*)PyArray_EMPTY(3,dims,NPY_DOUBLE,0); - - /* Grab pointers to the various data arrays */ - oPtr = (double*)PyArray_DATA(omeArray); - rPtr = (double*)PyArray_DATA(rMat); - - /* Call the actual function repeatedly */ - for (i = 0; i < no; ++i) { - makeOscillRotMat_cfunc(chi, oPtr[i], rPtr + i*9); - } - - return((PyObject*)rMat); -} - -static PyObject * makeRotMatOfExpMap(PyObject * self, PyObject * args) -{ - PyArrayObject *expMap, *rMat; - int de; - npy_intp ne, dims[2]; - double *ePtr, *rPtr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"O", &expMap)) return(NULL); - if ( expMap == NULL ) return(NULL); - - /* Verify shape of input arrays */ - de = PyArray_NDIM(expMap); - assert( de == 1 ); - - /* Verify dimensions of input arrays */ - ne = PyArray_DIMS(expMap)[0]; - assert( ne == 3 ); - - /* Allocate the result matrix with appropriate dimensions and type */ - dims[0] = 3; dims[1] = 3; - rMat = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - - /* Grab pointers to the various data arrays */ - ePtr = (double*)PyArray_DATA(expMap); - rPtr = (double*)PyArray_DATA(rMat); - - /* Call the actual function */ - makeRotMatOfExpMap_cfunc(ePtr,rPtr); - - return((PyObject*)rMat); -} - -static PyObject * makeRotMatOfQuat(PyObject * self, PyObject * args) -{ - PyArrayObject *quat, *rMat; - int nq, ne, de; - npy_intp dims2[2]={3, 3}, dims3[3]; - double *qPtr, *rPtr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"O", &quat)) return(NULL); - if ( quat == NULL ) return(NULL); - - /* Verify shape of input arrays */ - de = PyArray_NDIM(quat); - if (de == 1) { - ne = PyArray_DIMS(quat)[0]; - assert( ne == 4 ); - nq = 1; - /* Allocate the result matrix with appropriate dimensions and type */ - rMat = (PyArrayObject*)PyArray_EMPTY(2,dims2,NPY_DOUBLE,0); - } else { - assert( de == 2 ); - nq = PyArray_DIMS(quat)[0]; - ne = PyArray_DIMS(quat)[1]; - assert( ne == 4 ); - dims3[0] = nq; dims3[1] = 3; dims3[2] = 3; - /* Allocate the result matrix with appropriate dimensions and type */ - rMat = (PyArrayObject*)PyArray_EMPTY(3,dims3,NPY_DOUBLE,0); - } - - /* Grab pointers to the various data arrays */ - qPtr = (double*)PyArray_DATA(quat); - rPtr = (double*)PyArray_DATA(rMat); - - /* Call the actual function */ - makeRotMatOfQuat_cfunc(nq, qPtr, rPtr); - - return((PyObject*)rMat); -} - -static PyObject * makeBinaryRotMat(PyObject * self, PyObject * args) -{ - PyArrayObject *axis, *rMat; - int da; - npy_intp na, dims[2]; - double *aPtr, *rPtr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"O", &axis)) return(NULL); - if ( axis == NULL ) return(NULL); - - /* Verify shape of input arrays */ - da = PyArray_NDIM(axis); - assert( da == 1 ); - - /* Verify dimensions of input arrays */ - na = PyArray_DIMS(axis)[0]; - assert( na == 3 ); - - /* Allocate the result matrix with appropriate dimensions and type */ - dims[0] = 3; dims[1] = 3; - rMat = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - - /* Grab pointers to the various data arrays */ - aPtr = (double*)PyArray_DATA(axis); - rPtr = (double*)PyArray_DATA(rMat); - - /* Call the actual function */ - makeBinaryRotMat_cfunc(aPtr,rPtr); - - return((PyObject*)rMat); -} - -static PyObject * makeEtaFrameRotMat(PyObject * self, PyObject * args) -{ - PyArrayObject *bHat, *eHat, *rMat; - int db, de; - npy_intp nb, ne, dims[2]; - double *bPtr, *ePtr, *rPtr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OO", &bHat,&eHat)) return(NULL); - if ( bHat == NULL || eHat == NULL ) return(NULL); - - /* Verify shape of input arrays */ - db = PyArray_NDIM(bHat); - de = PyArray_NDIM(eHat); - assert( db == 1 && de == 1); - - /* Verify dimensions of input arrays */ - nb = PyArray_DIMS(bHat)[0]; - ne = PyArray_DIMS(eHat)[0]; - assert( nb == 3 && ne == 3 ); - - /* Allocate the result matrix with appropriate dimensions and type */ - dims[0] = 3; dims[1] = 3; - rMat = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - - /* Grab pointers to the various data arrays */ - bPtr = (double*)PyArray_DATA(bHat); - ePtr = (double*)PyArray_DATA(eHat); - rPtr = (double*)PyArray_DATA(rMat); - - /* Call the actual function */ - makeEtaFrameRotMat_cfunc(bPtr,ePtr,rPtr); - - return((PyObject*)rMat); -} - -static PyObject * validateAngleRanges(PyObject * self, PyObject * args) -{ - PyArrayObject *angList, *angMin, *angMax, *reflInRange; - PyObject *ccw; - int ccwVal = 1; /* ccwVal set to True by default */ - int da, dmin, dmax; - npy_intp na, nmin, nmax; - double *aPtr, *minPtr, *maxPtr; - bool *rPtr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OOOO", &angList,&angMin,&angMax,&ccw)) return(NULL); - if ( angList == NULL || angMin == NULL || angMax == NULL ) return(NULL); - - /* Verify shape of input arrays */ - da = PyArray_NDIM(angList); - dmin = PyArray_NDIM(angMin); - dmax = PyArray_NDIM(angMax); - assert( da == 1 && dmin == 1 && dmax ==1 ); - - /* Verify dimensions of input arrays */ - na = PyArray_DIMS(angList)[0]; - nmin = PyArray_DIMS(angMin)[0]; - nmax = PyArray_DIMS(angMax)[0]; - assert( nmin == nmax ); - - /* Check the value of ccw */ - if ( ccw == Py_True ) - ccwVal = 1; - else - ccwVal = 0; - - /* Allocate the result matrix with appropriate dimensions and type */ - reflInRange = (PyArrayObject*)PyArray_EMPTY(1,PyArray_DIMS(angList),NPY_BOOL,false); - assert( reflInRange != NULL ); - - /* Grab pointers to the various data arrays */ - aPtr = (double*)PyArray_DATA(angList); - minPtr = (double*)PyArray_DATA(angMin); - maxPtr = (double*)PyArray_DATA(angMax); - rPtr = (bool*)PyArray_DATA(reflInRange); - - /* Call the actual function */ - validateAngleRanges_cfunc(na,aPtr,nmin,minPtr,maxPtr,rPtr,ccwVal); - - return((PyObject*)reflInRange); -} - -static PyObject * rotate_vecs_about_axis(PyObject * self, PyObject * args) -{ - PyArrayObject *angles, *axes, *vecs; - PyArrayObject *rVecs; - int da, dax, dv; - npy_intp na, nax0, nax1, nv0, nv1; - double *aPtr, *axesPtr, *vecsPtr; - double *rPtr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OOO", &angles,&axes,&vecs)) return(NULL); - if ( angles == NULL || axes == NULL || vecs == NULL ) return(NULL); - - /* Verify shape of input arrays */ - da = PyArray_NDIM(angles); - dax = PyArray_NDIM(axes); - dv = PyArray_NDIM(vecs); - assert( da == 1 && dax == 2 && dv == 2 ); - - /* Verify dimensions of input arrays */ - na = PyArray_DIMS(angles)[0]; - nax0 = PyArray_DIMS(axes)[0]; - nax1 = PyArray_DIMS(axes)[1]; - nv0 = PyArray_DIMS(vecs)[0]; - nv1 = PyArray_DIMS(vecs)[1]; - assert( na == 1 || na == nv0 ); - assert( nax0 == 1 || nax0 == nv0 ); - assert( nax1 == 3 && nv1 == 3 ); - - /* Allocate the result vectors with appropriate dimensions and type */ - rVecs = (PyArrayObject*)PyArray_EMPTY(2,PyArray_DIMS(vecs),NPY_DOUBLE,0.0); - assert( rVecs != NULL ); - - /* Grab pointers to the various data arrays */ - aPtr = (double*)PyArray_DATA(angles); - axesPtr = (double*)PyArray_DATA(axes); - vecsPtr = (double*)PyArray_DATA(vecs); - rPtr = (double*)PyArray_DATA(rVecs); - - /* Call the actual function */ - rotate_vecs_about_axis_cfunc(na,aPtr,nax0,axesPtr,nv0,vecsPtr,rPtr); - - return((PyObject*)rVecs); -} - -static PyObject * quat_distance(PyObject * self, PyObject * args) -{ - PyArrayObject *q1, *q2, *qsym; - double *q1Ptr, *q2Ptr, *qsymPtr; - int dq1, dq2, dqsym; - int nq1, nq2, nqsym, nsym; - double dist = 0.0; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"OOO", &q1,&q2,&qsym)) return(NULL); - if ( q1 == NULL || q2 == NULL || qsym == NULL ) return(NULL); - - /* Verify shape of input arrays */ - dq1 = PyArray_NDIM(q1); - dq2 = PyArray_NDIM(q2); - dqsym = PyArray_NDIM(qsym); - assert( dq1 == 1 && dq2 == 1 && dqsym == 2 ); - - /* Verify dimensions of input arrays */ - nq1 = PyArray_DIMS(q1)[0]; - nq2 = PyArray_DIMS(q2)[0]; - nqsym = PyArray_DIMS(qsym)[0]; - nsym = PyArray_DIMS(qsym)[1]; - assert( nq1 == 4 && nq2 == 4 && nqsym == 4 ); - - /* Grab pointers to the various data arrays */ - q1Ptr = (double*)PyArray_DATA(q1); - q2Ptr = (double*)PyArray_DATA(q2); - qsymPtr = (double*)PyArray_DATA(qsym); - - /* Call the actual function */ - dist = quat_distance_cfunc(nsym,q1Ptr,q2Ptr,qsymPtr); - if (dist < 0) { - PyErr_SetString(PyExc_RuntimeError, "Could not allocate memory"); - return NULL; - } - return(PyFloat_FromDouble(dist)); -} - -static PyObject * homochoricOfQuat(PyObject * self, PyObject * args) -{ - PyArrayObject *quat, *hVec; - int nq, dq, ne; - npy_intp dims[2]; - double *qPtr, *hPtr; - - /* Parse arguments */ - if ( !PyArg_ParseTuple(args,"O", &quat)) return(NULL); - if ( quat == NULL ) return(NULL); - - /* Verify shape of input arrays */ - dq = PyArray_NDIM(quat); - if (dq == 1) { - ne = PyArray_DIMS(quat)[0]; - assert( ne == 4 ); - nq = 1; - dims[0] = nq; dims[1] = 3; - /* Allocate the result matrix with appropriate dimensions and type */ - hVec = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - } else { - assert( dq == 2 ); - nq = PyArray_DIMS(quat)[0]; - ne = PyArray_DIMS(quat)[1]; - assert( ne == 4 ); - dims[0] = nq; dims[1] = 3; - /* Allocate the result matrix with appropriate dimensions and type */ - hVec = (PyArrayObject*)PyArray_EMPTY(2,dims,NPY_DOUBLE,0); - } - - /* Grab pointers to the various data arrays */ - qPtr = (double*)PyArray_DATA(quat); - hPtr = (double*)PyArray_DATA(hVec); - - /* Call the actual function */ - homochoricOfQuat_cfunc(nq, qPtr, hPtr); - - return((PyObject*)hVec); -} diff --git a/hexrd/core/transforms/transforms_CAPI.h b/hexrd/core/transforms/transforms_CAPI.h deleted file mode 100644 index ae0aacf75..000000000 --- a/hexrd/core/transforms/transforms_CAPI.h +++ /dev/null @@ -1,79 +0,0 @@ -#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION - -#include -#include - -/******************************************************************************/ -/** The functions declared in this header make use of Python's C API to **/ -/** parse input arguments and verify their types (and sizes when **/ -/** appropriate). They also allocte new Python objects for use as return **/ -/** values. **/ -/** **/ -/** In short, these functions perform all of the necessary C API calls. **/ -/** **/ -/** The actual computations are handled by functions declared in the file **/ -/** "transforms_CFUNC.h". **/ -/** **/ -/** This separation of Python C API calls and C implementations allows the C **/ -/** functions to call eachother without the unnecessary overhead of passing **/ -/** arguments and return values as Python objects. **/ -/** **/ -/******************************************************************************/ - -/******************************************************************************/ -/* Funtions */ - -static PyObject * anglesToGVec(PyObject * self, PyObject * args); - -static PyObject * anglesToDVec(PyObject * self, PyObject * args); - -static PyObject * makeGVector(PyObject * self, PyObject * args); - -static PyObject * gvecToDetectorXY(PyObject * self, PyObject * args); - -static PyObject * gvecToDetectorXYArray(PyObject * self, PyObject * args); - -static PyObject * detectorXYToGvec(PyObject * self, PyObject * args); - -static PyObject * detectorXYToGvecArray(PyObject * self, PyObject * args); - -static PyObject * oscillAnglesOfHKLs(PyObject * self, PyObject * args); - -/******************************************************************************/ -/* Utility Funtions */ - -static PyObject * arccosSafe(PyObject * self, PyObject * args); - -static PyObject * angularDifference(PyObject * self, PyObject * args); - -static PyObject * mapAngle(PyObject * self, PyObject * args); - -static PyObject * columnNorm(PyObject * self, PyObject * args); - -static PyObject * rowNorm(PyObject * self, PyObject * args); - -static PyObject * unitRowVector(PyObject * self, PyObject * args); - -static PyObject * unitRowVectors(PyObject * self, PyObject * args); - -static PyObject * makeDetectorRotMat(PyObject * self, PyObject * args); - -static PyObject * makeOscillRotMat(PyObject * self, PyObject * args); - -static PyObject * makeOscillRotMatArray(PyObject * self, PyObject * args); - -static PyObject * makeRotMatOfExpMap(PyObject * self, PyObject * args); - -static PyObject * makeRotMatOfQuat(PyObject * self, PyObject * args); - -static PyObject * makeBinaryRotMat(PyObject * self, PyObject * args); - -static PyObject * makeEtaFrameRotMat(PyObject * self, PyObject * args); - -static PyObject * validateAngleRanges(PyObject * self, PyObject * args); - -static PyObject * rotate_vecs_about_axis(PyObject * self, PyObject * args); - -static PyObject * quat_distance(PyObject * self, PyObject * args); - -static PyObject * homochoricOfQuat(PyObject * self, PyObject * args); diff --git a/hexrd/core/transforms/transforms_CFUNC.c b/hexrd/core/transforms/transforms_CFUNC.c deleted file mode 100644 index f9b233c33..000000000 --- a/hexrd/core/transforms/transforms_CFUNC.c +++ /dev/null @@ -1,1521 +0,0 @@ -#define _USE_MATH_DEFINES -#include -#include -#include -#include - -#include "transforms_CFUNC.h" - -/* - * Microsoft's C compiler, when running in C mode, does not support the inline - * keyword. However it does support an __inline one. - * - * So if compiling with MSC, just use __inline as inline - */ -#if defined(_MSC_VER) -# define inline __inline -#endif - -/* - * For now, disable C99 codepaths - */ -#define USE_C99_CODE 0 -#if ! defined(USE_C99_CODE) -# if defined(__STDC__) -# if (__STD_VERSION__ >= 199901L) -# define USE_C99_CODE 1 -# else -# define USE_C99_CODE 0 -# endif -# endif -#endif - -#if ! USE_C99_CODE -/* - * Just remove any "restrict" keyword that may be present. - */ -#define restrict -#endif - -static double epsf = 2.2e-16; -static double sqrt_epsf = 1.5e-8; -static double Zl[3] = {0.0,0.0,1.0}; - - -/******************************************************************************/ -/* Functions */ -#if USE_C99_CODE -static inline -double * -v3_v3s_inplace_add(double *dst_src1, - const double *src2, int stride) -{ - dst_src1[0] += src2[0]; - dst_src1[1] += src2[1*stride]; - dst_src1[2] += src2[2*stride]; - return dst_src1; -} - -static inline -double * -v3_v3s_add(const double *src1, - const double *src2, int stride, - double * restrict dst) -{ - dst[0] = src1[0] + src2[0]; - dst[1] = src1[1] + src2[1*stride]; - dst[2] = src1[2] + src2[2*stride]; - - return dst; -} - -static inline -double * -v3_v3s_inplace_sub(double *dst_src1, - const double *src2, int stride) -{ - dst_src1[0] -= src2[0]; - dst_src1[1] -= src2[1*stride]; - dst_src1[2] -= src2[2*stride]; - return dst_src1; -} - -static inline -double * -v3_v3s_sub(const double *src1, - const double *src2, int stride, - double * restrict dst) -{ - dst[0] = src1[0] - src2[0]; - dst[1] = src1[1] - src2[1*stride]; - dst[2] = src1[2] - src2[2*stride]; - - return dst; -} - -static inline -double * -v3_inplace_normalize(double * restrict v) -{ - double sqr_norm = v[0]*v[0] + v[1]*v[1] + v[2]*v[2]; - - if (sqr_norm > epsf) { - double normalize_factor = 1./sqrt(sqr_norm); - v[0] *= normalize_factor; - v[1] *= normalize_factor; - v[2] *= normalize_factor; - } - - return v; -} - -static inline -double * -v3_normalize(const double *in, - double * restrict out) -{ - double in0 = in[0], in1 = in[1], in2 = in[2]; - double sqr_norm = in0*in0 + in1*in1 + in2*in2; - - if (sqr_norm > epsf) { - double normalize_factor = 1./sqrt(sqr_norm); - out[0] = in0 * normalize_factor; - out[1] = in1 * normalize_factor; - out[2] = in2 * normalize_factor; - } else { - out[0] = in0; - out[1] = in1; - out[2] = in2; - } - - return out; -} - -static inline -double * -m33_inplace_transpose(double * restrict m) -{ - double e1 = m[1]; - double e2 = m[2]; - double e5 = m[5]; - m[1] = m[3]; - m[2] = m[6]; - m[5] = m[7]; - m[3] = e1; - m[6] = e2; - m[7] = e5; - - return m; -} - -static inline -double * -m33_transpose(const double *m, - double * restrict dst) -{ - dst[0] = m[0]; dst[1] = m[3]; dst[2] = m[6]; - dst[3] = m[1]; dst[4] = m[4]; dst[5] = m[7]; - dst[7] = m[2]; dst[8] = m[5]; dst[9] = m[9]; - - return dst; -} - -static inline -double -v3_v3s_dot(const double *v1, - const double *v2, int stride) -{ - return v1[0]*v2[0] + v1[1]*v2[stride] + v1[2]*v2[2*stride]; -} - - -/* 3x3 matrix by strided 3 vector product ------------------------------------- - hopefully a constant stride will be optimized - */ -static inline -double * -m33_v3s_multiply(const double *m, - const double *v, int stride, - double * restrict dst) -{ - dst[0] = m[0]*v[0] + m[1]*v[stride] + m[2]*v[2*stride]; - dst[1] = m[3]*v[0] + m[4]*v[stride] + m[5]*v[2*stride]; - dst[2] = m[6]*v[0] + m[7]*v[stride] + m[8]*v[2*stride]; - - return dst; -} - -/* transposed 3x3 matrix by strided 3 vector product -------------------------- - */ -static inline -double * -v3s_m33t_multiply(const double *v, int stride, - const double *m, - double * restrict dst) -{ - double v0 = v[0]; double v1 = v[stride]; double v2 = v[2*stride]; - dst[0] = v0*m[0] + v1*m[1] + v2*m[2]; - dst[1] = v0*m[3] + v1*m[4] + v2*m[5]; - dst[2] = v0*m[6] + v1*m[7] + v2*m[8]; - - return dst; -} - -static inline -double * -v3s_m33_multiply(const double *v, int stride, - const double *m, - double * restrict dst) -{ - double v0 = v[0]; double v1 = v[stride]; double v2 = v[2*stride]; - dst[0] = v0*m[0] + v1*m[3] + v2*m[6]; - dst[1] = v0*m[1] + v1*m[4] + v2*m[7]; - dst[2] = v0*m[2] + v1*m[5] + v2*m[8]; - - return dst; -} - -static inline -double * -m33t_v3s_multiply(const double *m, - const double *v, int stride, - double * restrict dst) -{ - dst[0] = m[0]*v[0] + m[3]*v[stride] + m[6]*v[2*stride]; - dst[1] = m[1]*v[0] + m[4]*v[stride] + m[7]*v[2*stride]; - dst[2] = m[2]*v[0] + m[5]*v[stride] + m[8]*v[2*stride]; - - return dst; -} - -static inline -double * -m33_m33_multiply(const double *src1, - const double *src2, - double * restrict dst) -{ - v3s_m33_multiply(src1 + 0, 1, src2, dst+0); - v3s_m33_multiply(src1 + 3, 1, src2, dst+3); - v3s_m33_multiply(src1 + 6, 1, src2, dst+6); - - return dst; -} - -static inline -double * -m33t_m33_multiply(const double *src1, - const double *src2, - double * restrict dst) -{ - v3s_m33_multiply(src1 + 0, 3, src2, dst+0); - v3s_m33_multiply(src1 + 1, 3, src2, dst+3); - v3s_m33_multiply(src1 + 2, 3, src2, dst+6); - - return dst; -} - -static inline -double * -m33_m33t_multiply(const double *src1, - const double *src2, - double * restrict dst) -{ - return m33_inplace_transpose(m33t_m33_multiply(src2, src1, dst)); -} - -static inline -double * -m33t_m33t_multiply(const double *src1, - const double *src2, - double * restrict dst) -{ - return m33_inplace_transpose(m33_m33_multiply(src2, src1, dst)); -} - -#endif - -#if USE_C99_CODE -static inline -void anglesToGvec_single(double *v3_ang, double *m33_e, - double chi, double *m33_c, - double * restrict v3_c) -{ - double v3_g[3], v3_tmp1[3], v3_tmp2[3], m33_s[9], m33_ctst[9]; - - /* build g */ - double cx = cos(0.5*v3_ang[0]); - double sx = sin(0.5*v3_ang[0]); - double cy = cos(v3_ang[1]); - double sy = sin(v3_ang[1]); - v3_g[0] = cx*cy; - v3_g[1] = cx*sy; - v3_g[2] = sx; - - /* build S */ - makeOscillRotMat_cfunc(chi, v3_ang[2], m33_s); - - /* beam frame to lab frame */ - /* eval the chain: - C.T _dot_ S.T _dot_ E _dot_ g - */ - m33_v3s_multiply (m33_e, v3_g, 1, v3_tmp1); /* E _dot_ g */ - m33t_v3s_multiply(m33_s, v3_tmp1, 1, v3_tmp2); /* S.T _dot_ E _dot_ g */ - m33t_v3s_multiply(m33_c, v3_tmp2, 1, v3_c); /* the whole lot */ -} - -void anglesToGvec_cfunc(long int nvecs, double * angs, - double * bHat_l, double * eHat_l, - double chi, double * rMat_c, - double * gVec_c) -{ - double m33_e[9]; - - makeEtaFrameRotMat_cfunc(bHat_l, eHat_l, m33_e); - - for (int i = 0; i= ztol && bDot <= 1.0-ztol ) { - /* - * If we are here diffraction is possible so increment the number of - * admissable vectors - */ - double brMat[9]; - makeBinaryRotMat_cfunc(gVec_l, brMat); - - double dVec_l[3]; - m33_v3s_multiply(brMat, bHat_l, 1, dVec_l); - double denom = v3_v3s_dot(nVec_l, dVec_l, 1); - - if (denom > ztol) { - double u = num/denom; - double v3_tmp[3]; - - /* v3_tmp = P0_l + u*dVec_l - tVec_d */ - for (int j=0; j<3; j++) - v3_tmp[j] = P0_l[j] + u*dVec_l[j] - tVec_d[j]; - - result[0] = v3_v3s_dot(v3_tmp, rMat_d + 0, 3); - result[1] = v3_v3s_dot(v3_tmp, rMat_d + 1, 3); - - /* result when computation can be finished */ - return; - } - } - - /* default result when computation can't be finished */ - result[0] = NAN; - result[1] = NAN; -} - -/* - * The only difference between this and the non-Array version - * is that rMat_s is an array of matrices of length npts instead - * of a single matrix. - */ -void gvecToDetectorXYArray_cfunc(long int npts, double * gVec_c_array, - double * rMat_d, double * rMat_s_array, double * rMat_c, - double * tVec_d, double * tVec_s, double * tVec_c, - double * beamVec, double * result_array) -{ - /* Normalize the beam vector */ - double bHat_l[3]; - v3_normalize(beamVec, bHat_l); - double nVec_l[3]; - m33_v3s_multiply(rMat_d, Zl, 1, nVec_l); - - for (size_t i = 0; i < npts; i++) { - double *rMat_s = rMat_s_array + 9*i; - double *gVec_c = gVec_c_array + 3*i; - double * restrict result = result_array + 2*i; - /* Initialize the detector normal and frame origins */ - - double P0_l[3]; - m33_v3s_multiply(rMat_s, tVec_c, 1, P0_l); - v3_v3s_inplace_add(P0_l, tVec_s, 1); - - double P3_l_minus_P0_l[3]; - v3_v3s_sub(tVec_d, P0_l, 1, P3_l_minus_P0_l); - double num = v3_v3s_dot(nVec_l, P3_l_minus_P0_l, 1); - - double gHat_c[3]; - v3_normalize(gVec_c, gHat_c); - /* - double rMat_sc[9]; - m33_m33_multiply(rMat_s, rMat_c, rMat_sc); - double gVec_l[3]; - m33_v3s_multiply(rMat_sc, gHat_c, 1, gVec_l); - */ - double tmp_vec[3], gVec_l[3]; - m33_v3s_multiply(rMat_c, gHat_c, 1, tmp_vec); - m33_v3s_multiply(rMat_s, tmp_vec, 1, gVec_l); - - double bDot = -v3_v3s_dot(bHat_l, gVec_l, 1); - double ztol = epsf; - - if (bDot < ztol || bDot > 1.0-ztol) { - result[0] = NAN; result[1] = NAN; - continue; - } - - double brMat[9]; - makeBinaryRotMat_cfunc(gVec_l, brMat); - - double dVec_l[3]; - m33_v3s_multiply(brMat, bHat_l, 1, dVec_l); - double denom = v3_v3s_dot(nVec_l, dVec_l, 1); - if (denom < ztol) { - result[0] = NAN; result[1] = NAN; - continue; - } - - double u = num/denom; - double v3_tmp[3]; - for (int j=0; j < 3; j++) - v3_tmp[j] = u*dVec_l[j] - P3_l_minus_P0_l[j]; - - result[0] = v3_v3s_dot(v3_tmp, rMat_d + 0, 3); - result[1] = v3_v3s_dot(v3_tmp, rMat_d + 1, 3); - } -} - -#else -void gvecToDetectorXYOne_cfunc(double * gVec_c, double * rMat_d, - double * rMat_sc, double * tVec_d, - double * bHat_l, - double * nVec_l, double num, double * P0_l, - double * result) -{ - int j, k; - double bDot, ztol, denom, u; - double gHat_c[3], gVec_l[3], dVec_l[3], P2_l[3], P2_d[3]; - double brMat[9]; - - ztol = epsf; - - /* Compute unit reciprocal lattice vector in crystal frame w/o - translation */ - unitRowVector_cfunc(3, gVec_c, gHat_c); - - /* Compute unit reciprocal lattice vector in lab frame and dot with beam - vector */ - bDot = 0.0; - for (j=0; j<3; j++) { - gVec_l[j] = 0.0; - for (k=0; k<3; k++) - gVec_l[j] += rMat_sc[3*j+k]*gHat_c[k]; - - bDot -= bHat_l[j]*gVec_l[j]; - } - - if ( bDot >= ztol && bDot <= 1.0-ztol ) { - /* If we are here diffraction is possible so increment the number of - admissable vectors */ - makeBinaryRotMat_cfunc(gVec_l, brMat); - - denom = 0.0; - for (j=0; j<3; j++) { - dVec_l[j] = 0.0; - for (k=0; k<3; k++) - dVec_l[j] -= brMat[3*j+k]*bHat_l[k]; - - denom += nVec_l[j]*dVec_l[j]; - } - - if ( denom < -ztol ) { - - u = num/denom; - - for (j=0; j<3; j++) - P2_l[j] = P0_l[j]+u*dVec_l[j]; - - for (j=0; j<2; j++) { - P2_d[j] = 0.0; - for (k=0; k<3; k++) - P2_d[j] += rMat_d[3*k+j]*(P2_l[k]-tVec_d[k]); - result[j] = P2_d[j]; - } - /* result when computation can be finished */ - return; - } - } - /* default result when computation can't be finished */ - result[0] = NAN; - result[1] = NAN; -} - -/* - * The only difference between this and the non-Array version - * is that rMat_s is an array of matrices of length npts instead - * of a single matrix. - */ -void gvecToDetectorXYArray_cfunc(long int npts, double * gVec_c, - double * rMat_d, double * rMat_s, - double * rMat_c, double * tVec_d, - double * tVec_s, double * tVec_c, - double * beamVec, double * result) -{ - long int i; - int j, k, l; - double num; - double nVec_l[3], bHat_l[3], P0_l[3], P3_l[3]; - double rMat_sc[9]; - - /* Normalize the beam vector */ - unitRowVector_cfunc(3,beamVec,bHat_l); - - for (i=0L; i < npts; i++) { - /* Initialize the detector normal and frame origins */ - num = 0.0; - for (j=0; j<3; j++) { - nVec_l[j] = 0.0; - P0_l[j] = tVec_s[j]; - - for (k=0; k<3; k++) { - nVec_l[j] += rMat_d[3*j+k]*Zl[k]; - P0_l[j] += rMat_s[9*i + 3*j+k]*tVec_c[k]; - } - - P3_l[j] = tVec_d[j]; - - num += nVec_l[j]*(P3_l[j]-P0_l[j]); - } - - /* Compute the matrix product of rMat_s and rMat_c */ - for (j=0; j<3; j++) { - for (k=0; k<3; k++) { - rMat_sc[3*j+k] = 0.0; - for (l=0; l<3; l++) { - rMat_sc[3*j+k] += rMat_s[9*i + 3*j+l]*rMat_c[3*l+k]; - } - } - } - - gvecToDetectorXYOne_cfunc(gVec_c + 3*i, rMat_d, rMat_sc, - tVec_d, bHat_l, nVec_l, num, - P0_l, result + 2*i); - } -} - -#endif -void gvecToDetectorXY_cfunc(long int npts, double * gVec_c, - double * rMat_d, double * rMat_s, double * rMat_c, - double * tVec_d, double * tVec_s, double * tVec_c, - double * beamVec, double * result) -{ - long int i; - int j, k, l; - - double num; - double nVec_l[3], bHat_l[3], P0_l[3], P3_l[3]; - double rMat_sc[9]; - - /* Normalize the beam vector */ - unitRowVector_cfunc(3,beamVec,bHat_l); - - /* Initialize the detector normal and frame origins */ - num = 0.0; - for (j=0; j<3; j++) { - nVec_l[j] = 0.0; - P0_l[j] = tVec_s[j]; - - for (k=0; k<3; k++) { - nVec_l[j] += rMat_d[3*j+k]*Zl[k]; - P0_l[j] += rMat_s[3*j+k]*tVec_c[k]; - } - - P3_l[j] = tVec_d[j]; - - num += nVec_l[j]*(P3_l[j]-P0_l[j]); - } - - /* Compute the matrix product of rMat_s and rMat_c */ - for (j=0; j<3; j++) { - for (k=0; k<3; k++) { - rMat_sc[3*j+k] = 0.0; - for (l=0; l<3; l++) { - rMat_sc[3*j+k] += rMat_s[3*j+l]*rMat_c[3*l+k]; - } - } - } - - for (i=0L; i epsf ) { - double nrm_factor = 1.0/sqrt(nrm); - for (j=0; j<3; j++) { - dHat_l[j] *= nrm_factor; - } - } - - /* Compute tTh */ - b_dot_dHat_l = 0.0; - for (j=0; j<3; j++) { - b_dot_dHat_l += bVec[j]*dHat_l[j]; - } - tTh = acos(b_dot_dHat_l); - - /* Compute eta */ - for (j=0; j<2; j++) { - tVec2[j] = 0.0; - for (k=0; k<3; k++) { - tVec2[j] += rMat_e[3*k+j]*dHat_l[k]; - } - } - eta = atan2(tVec2[1], tVec2[0]); - - /* Compute n_g vector */ - nrm = 0.0; - for (j=0; j<3; j++) { - double val; - int j1 = j < 2 ? j+1 : 0; - int j2 = j > 0 ? j-1 : 2; - val = bVec[j1] * dHat_l[j2] - bVec[j2] * dHat_l[j1]; - nrm += val*val; - n_g[j] = val; - } - if ( nrm > epsf ) { - double nrm_factor = 1.0/sqrt(nrm); - for (j=0; j<3; j++) { - n_g[j] *= nrm_factor; - } - } - - /* Rotate dHat_l vector */ - phi = 0.5*(M_PI-tTh); - *tTh_out = tTh; - *eta_out = eta; - rotate_vecs_about_axis_cfunc(1, &phi, 1, n_g, 1, dHat_l, gVec_l_out); -} - -void detectorXYToGvec_cfunc(long int npts, double * xy, - double * rMat_d, double * rMat_s, - double * tVec_d, double * tVec_s, double * tVec_c, - double * beamVec, double * etaVec, - double * tTh, double * eta, double * gVec_l) -{ - long int i; - int j, k; - double nrm, bVec[3], tVec1[3]; - double rMat_e[9]; - - /* Fill rMat_e */ - makeEtaFrameRotMat_cfunc(beamVec,etaVec,rMat_e); - - /* Normalize the beam vector */ - nrm = 0.0; - for (j=0; j<3; j++) { - nrm += beamVec[j]*beamVec[j]; - } - - if ( nrm > epsf ) { - double nrm_factor = 1.0/sqrt(nrm); - for (j=0; j<3; j++) - bVec[j] = beamVec[j]*nrm_factor; - } else { - for (j=0; j<3; j++) - bVec[j] = beamVec[j]; - } - - /* Compute shift vector */ - for (j=0; j<3; j++) { - tVec1[j] = tVec_d[j]-tVec_s[j]; - for (k=0; k<3; k++) { - tVec1[j] -= rMat_s[3*j+k]*tVec_c[k]; - } - } - - for (i=0; i epsf ) { - double nrm_factor = 1.0/sqrt(nrm); - for (j=0; j<3; j++) - bVec[j] = beamVec[j]*nrm_factor; - } else { - for (j=0; j<3; j++) - bVec[j] = beamVec[j]; - } - - for (j=0; j<3; j++) { - tVec1[j] = tVec_d[j]-tVec_s[j]; - for (k=0; k<3; k++) { - tVec1[j] -= rMat_s[3*j+k]*tVec_c[k]; - } - } - - for (i=0; i epsf ) { - for (j=0; j<3; j++) - bHat_l[j] = beamVec[j]/nrm0; - } else { - for (j=0; j<3; j++) - bHat_l[j] = beamVec[j]; - } - - /* Normalize the eta vector */ - nrm0 = 0.0; - for (j=0; j<3; j++) { - nrm0 += etaVec[j]*etaVec[j]; - } - nrm0 = sqrt(nrm0); - if ( nrm0 > epsf ) { - for (j=0; j<3; j++) - eHat_l[j] = etaVec[j]/nrm0; - } else { - for (j=0; j<3; j++) - eHat_l[j] = etaVec[j]; - } - - /* Check for consistent reference coordiantes */ - nrm0 = 0.0; - for (j=0; j<3; j++) { - nrm0 += bHat_l[j]*eHat_l[j]; - } - if ( fabs(nrm0) < 1.0-sqrt_epsf ) crc = true; - - /* Compute the sine and cosine of the oscillation axis tilt */ - cchi = cos(chi); - schi = sin(chi); - - for (i=0; i epsf ) { - for (j=0; j<3; j++) { - gHat_c[j] /= nrm0; - gHat_s[j] = tmpVec[j]/nrm0; - } - } - - /* Compute the sine of the Bragg angle */ - sintht = 0.5*wavelength*nrm0; - - /* Compute the coefficients of the harmonic equation */ - a = gHat_s[2]*bHat_l[0] + schi*gHat_s[0]*bHat_l[1] - cchi*gHat_s[0]*bHat_l[2]; - b = gHat_s[0]*bHat_l[0] - schi*gHat_s[2]*bHat_l[1] + cchi*gHat_s[2]*bHat_l[2]; - c = - sintht - cchi*gHat_s[1]*bHat_l[1] - schi*gHat_s[1]*bHat_l[2]; - - /* Form solution */ - abMag = sqrt(a*a + b*b); assert( abMag > 0.0 ); - phaseAng = atan2(b,a); - rhs = c/abMag; - - if ( fabs(rhs) > 1.0 ) { - for (j=0; j<3; j++) - oangs0[3L*i+j] = NAN; - for (j=0; j<3; j++) - oangs1[3L*i+j] = NAN; - continue; - } - - rhsAng = asin(rhs); - - /* Write ome angles */ - oangs0[3L*i+2] = rhsAng - phaseAng; - oangs1[3L*i+2] = M_PI - rhsAng - phaseAng; - - if ( crc ) { - makeEtaFrameRotMat_cfunc(bHat_l,eHat_l,rMat_e); - - oVec[0] = chi; - - oVec[1] = oangs0[3L*i+2]; - makeOscillRotMat_cfunc(oVec[0], oVec[1], rMat_s); - - for (j=0; j<3; j++) { - tVec0[j] = 0.0; - for (k=0; k<3; k++) { - tVec0[j] += rMat_s[3*j+k]*gHat_s[k]; - } - } - for (j=0; j<2; j++) { - gVec_e[j] = 0.0; - for (k=0; k<3; k++) { - gVec_e[j] += rMat_e[3*k+j]*tVec0[k]; - } - } - oangs0[3L*i+1] = atan2(gVec_e[1],gVec_e[0]); - - oVec[1] = oangs1[3L*i+2]; - makeOscillRotMat_cfunc(oVec[0], oVec[1], rMat_s); - - for (j=0; j<3; j++) { - tVec0[j] = 0.0; - for (k=0; k<3; k++) { - tVec0[j] += rMat_s[3*j+k]*gHat_s[k]; - } - } - for (j=0; j<2; j++) { - gVec_e[j] = 0.0; - for (k=0; k<3; k++) { - gVec_e[j] += rMat_e[3*k+j]*tVec0[k]; - } - } - oangs1[3L*i+1] = atan2(gVec_e[1],gVec_e[0]); - - oangs0[3L*i+0] = 2.0*asin(sintht); - oangs1[3L*i+0] = oangs0[3L*i+0]; - } - } -} - -/******************************************************************************/ -/* Utility Funtions */ - -void unitRowVector_cfunc(int n, double * cIn, double * cOut) -{ - int j; - double nrm; - - nrm = 0.0; - for (j=0; j epsf ) { - for (j=0; j epsf ) { - for (j=0; j epsf ) { - s = sin(phi)/phi; - c = (1.0-cos(phi))/(phi*phi); - - rPtr[1] -= s*ePtr[2]; - rPtr[2] += s*ePtr[1]; - rPtr[3] += s*ePtr[2]; - rPtr[5] -= s*ePtr[0]; - rPtr[6] -= s*ePtr[1]; - rPtr[7] += s*ePtr[0]; - - rPtr[1] += c*ePtr[0]*ePtr[1]; - rPtr[2] += c*ePtr[0]*ePtr[2]; - rPtr[3] += c*ePtr[1]*ePtr[0]; - rPtr[5] += c*ePtr[1]*ePtr[2]; - rPtr[6] += c*ePtr[2]*ePtr[0]; - rPtr[7] += c*ePtr[2]*ePtr[1]; - - rPtr[0] -= c*(ePtr[1]*ePtr[1]+ePtr[2]*ePtr[2]); - rPtr[4] -= c*(ePtr[2]*ePtr[2]+ePtr[0]*ePtr[0]); - rPtr[8] -= c*(ePtr[0]*ePtr[0]+ePtr[1]*ePtr[1]); - } -} - -void makeRotMatOfQuat_cfunc(int nq, double * qPtr, double * rPtr) -{ - int i, j; - double c, s, phi, n[3]={0.0,0.0,0.0}; - - for (i=0; i epsf) { - n[0] = (1. / sin(0.5*phi)) * qPtr[4*i+1]; - n[1] = (1. / sin(0.5*phi)) * qPtr[4*i+2]; - n[2] = (1. / sin(0.5*phi)) * qPtr[4*i+3]; - - s = sin(phi); - c = cos(phi); - - rPtr[9*i+0] = c + n[0]*n[0]*(1. - c); - rPtr[9*i+1] = n[0]*n[1]*(1. - c) - n[2]*s; - rPtr[9*i+2] = n[0]*n[2]*(1. - c) + n[1]*s; - rPtr[9*i+3] = n[1]*n[0]*(1. - c) + n[2]*s; - rPtr[9*i+4] = c + n[1]*n[1]*(1. - c); - rPtr[9*i+5] = n[1]*n[2]*(1. - c) - n[0]*s; - rPtr[9*i+6] = n[2]*n[0]*(1. - c) - n[1]*s; - rPtr[9*i+7] = n[2]*n[1]*(1. - c) + n[0]*s; - rPtr[9*i+8] = c + n[2]*n[2]*(1. - c); - } - else { - for (j=0; j<9; j++) { - if ( j%4 == 0 ) - rPtr[9*i+j] = 1.0; - else - rPtr[9*i+j] = 0.0; - } - } - } -} - -void makeBinaryRotMat_cfunc(double * aPtr, double * rPtr) -{ - int i, j; - - for (i=0; i<3; i++) { - for (j=0; j<3; j++) { - rPtr[3*i+j] = 2.0*aPtr[i]*aPtr[j]; - } - rPtr[3*i+i] -= 1.0; - } -} - -void makeEtaFrameRotMat_cfunc(double * bPtr, double * ePtr, double * rPtr) -{ - /* - * This function generates a COB matrix that takes components in BEAM frame to - * LAB frame - * - * NOTE: the beam and eta vectors MUST NOT BE COLINEAR!!!! - */ - int i; - - double yPtr[3], bHat[3], yHat[3], xHat[3]; - - /* find Y as e ^ b */ - yPtr[0] = ePtr[1]*bPtr[2] - bPtr[1]*ePtr[2]; - yPtr[1] = ePtr[2]*bPtr[0] - bPtr[2]*ePtr[0]; - yPtr[2] = ePtr[0]*bPtr[1] - bPtr[0]*ePtr[1]; - - /* Normalize beam (-Z) and Y vectors */ - unitRowVector_cfunc(3, bPtr, bHat); - unitRowVector_cfunc(3, yPtr, yHat); - - /* Find X as b ^ Y */ - xHat[0] = bHat[1]*yHat[2] - yHat[1]*bHat[2]; - xHat[1] = bHat[2]*yHat[0] - yHat[2]*bHat[0]; - xHat[2] = bHat[0]*yHat[1] - yHat[0]*bHat[1]; - - /* Assign columns */ - /* Assign Y column */ - for (i=0; i<3; i++) - { - rPtr[3*i+0] = xHat[i]; - rPtr[3*i+1] = yHat[i]; - rPtr[3*i+2] = -bHat[i]; - } -} - -void validateAngleRanges_old_cfunc(int na, double * aPtr, int nr, double * minPtr, double * maxPtr, bool * rPtr) -{ - int i, j; - double thetaMax, theta; - - /* Each angle should only be examined once. Any more is a waste of time. */ - for (i=0; i 2.0*M_PI ) - thetaMax -= 2.0*M_PI; - - while ( theta < 0.0 ) - theta += 2.0*M_PI; - while ( theta > 2.0*M_PI ) - theta -= 2.0*M_PI; - - if ( theta > -sqrt_epsf && theta < thetaMax + sqrt_epsf ) { - rPtr[i] = true; - - /* No need to check other ranges */ - break; - } - } - } -} - -void validateAngleRanges_cfunc(int na, double * aPtr, int nr, - double * minPtr, double * maxPtr, - bool * rPtr, int ccw) -{ - int i, j; - double thetaMax, theta; - double *startPtr, *stopPtr; - - if ( ccw ) { - startPtr = minPtr; - stopPtr = maxPtr; - } else { - startPtr = maxPtr; - stopPtr = minPtr; - } - - /* Each angle should only be examined once. Any more is a waste of time. */ - for (i=0; i 2.0*M_PI ) - thetaMax -= 2.0*M_PI; - - /* Check for an empty range */ - if ( fabs(thetaMax) < sqrt_epsf ) { - rPtr[i] = true; - - /* No need to check other ranges */ - break; - } - - /* Check for a range which spans a full circle */ - if ( fabs(thetaMax-2.0*M_PI) < sqrt_epsf ) { - - /* Double check the initial range */ - if ( (ccw && maxPtr[j] > minPtr[j]) || ((!ccw) && maxPtr[j] < minPtr[j]) ) { - rPtr[i] = true; - - /* No need to check other ranges */ - break; - } - } - - while ( theta < 0.0 ) - theta += 2.0*M_PI; - while ( theta > 2.0*M_PI ) - theta -= 2.0*M_PI; - - if ( theta >= -sqrt_epsf && theta <= thetaMax+sqrt_epsf ) { - rPtr[i] = true; - - /* No need to check other ranges */ - break; - } - } - } -} - - -void rotate_vecs_about_axis_cfunc(long int na, double * angles, - long int nax, double * axes, - long int nv, double * vecs, - double * rVecs) -{ - int i, j, sa, sax; - double c, s, nrm, proj, aCrossV[3]; - - if ( na == 1 ) sa = 0; - else sa = 1; - if ( nax == 1 ) sax = 0; - else sax = 3; - - for (i=0; i 1 || i == 0 ) { - nrm = 0.0; - for (j=0; j<3; j++) - nrm += axes[sax*i+j]*axes[sax*i+j]; - nrm = sqrt(nrm); - } - - /* Compute projection of vec along axis */ - proj = 0.0; - for (j=0; j<3; j++) - proj += axes[sax*i+j]*vecs[3*i+j]; - - /* Compute the cross product of the axis with vec */ - for (j=0; j<3; j++) - aCrossV[j] = axes[sax*i+(j+1)%3]*vecs[3*i+(j+2)%3]-axes[sax*i+(j+2)%3]*vecs[3*i+(j+1)%3]; - - /* Combine the three terms to compute the rotated vector */ - for (j=0; j<3; j++) { - rVecs[3*i+j] = c*vecs[3*i+j]+(s/nrm)*aCrossV[j]+(1.0-c)*proj*axes[sax*i+j]/(nrm*nrm); - } - } -} - -double quat_distance_cfunc(int nsym, double * q1, double * q2, double * qsym) -{ - int i; - double q0, q0_max = 0.0, dist = 0.0; - double *q2s; - - if ( NULL == (q2s = (double *)malloc(4*nsym*sizeof(double))) ) { - printf("malloc failed\n"); - return(-1); - } - - /* For each symmetry in qsym compute its inner product with q2 */ - for (i=0; i q0_max ) { - q0_max = fabs(q0); - } - } - - if ( q0_max <= 1.0 ) - dist = 2.0*acos(q0_max); - else if ( q0_max - 1. < 1e-12 ) - /* in case of quats loaded from single precision file */ - dist = 0.; - else - dist = NAN; - - free(q2s); - - return(dist); -} - -void homochoricOfQuat_cfunc(int nq, double * qPtr, double * hPtr) -{ - int i; - double arg, f, s, phi; - - for (i=0; i epsf) { - arg = 0.75*(phi - sin(phi)); - if (arg < 0.) { - f = -pow(-arg, 1./3.); - } else { - f = pow(arg, 1./3.); - } - s = 1. / sin(0.5*phi); - - hPtr[3*i+0] = f * s * qPtr[4*i+1]; - hPtr[3*i+1] = f * s * qPtr[4*i+2]; - hPtr[3*i+2] = f * s * qPtr[4*i+3]; - } - else { - hPtr[3*i+0] = 0.; - hPtr[3*i+1] = 0.; - hPtr[3*i+2] = 0.; - } - } -} diff --git a/hexrd/core/transforms/transforms_CFUNC.h b/hexrd/core/transforms/transforms_CFUNC.h deleted file mode 100644 index 628883f7d..000000000 --- a/hexrd/core/transforms/transforms_CFUNC.h +++ /dev/null @@ -1,95 +0,0 @@ -/******************************************************************************/ -/** The functions declared in this header use standard C to implement the **/ -/** various computations needed by the functions declared in the file **/ -/** "transforms_CAPI.h". **/ -/** **/ -/** This separation of Python C API calls and C implementations allows the C **/ -/** functions to call eachother without the unnecessary overhead of passing **/ -/** arguments and return values as Python objects. **/ -/** **/ -/******************************************************************************/ - - -#ifdef _WIN32 - #define false 0 - #define true 1 - #define bool int - #ifndef NAN - static const unsigned long __nan[2] = {0xffffffff, 0x7fffffff}; - #define NAN (*(const float *) __nan) - #endif -#else - #include -#endif - -/******************************************************************************/ -/* Funtions */ - -void anglesToGvec_cfunc(long int nvecs, double * angs, - double * bHat_l, double * eHat_l, - double chi, double * rMat_c, - double * gVec_c); - -void anglesToDvec_cfunc(long int nvecs, double * angs, - double * bHat_l, double * eHat_l, - double chi, double * rMat_c, - double * gVec_c); - -void gvecToDetectorXY_cfunc(long int npts, double * gVec_c_Ptr, - double * rMat_d_Ptr, double * rMat_s_Ptr, double * rMat_c_Ptr, - double * tVec_d_Ptr, double * tVec_s_Ptr, double * tVec_c_Ptr, - double * beamVec_Ptr, - double * result_Ptr); - -void gvecToDetectorXYArray_cfunc(long int npts, double * gVec_c_Ptr, - double * rMat_d_Ptr, double * rMat_s_Ptr, double * rMat_c_Ptr, - double * tVec_d_Ptr, double * tVec_s_Ptr, double * tVec_c_Ptr, - double * beamVec_Ptr, - double * result_Ptr); - -void detectorXYToGvec_cfunc(long int npts, double * xy, - double * rMat_d, double * rMat_s, - double * tVec_d, double * tVec_s, double * tVec_c, - double * beamVec, double * etaVec, - double * tTh, double * eta, double * gVec_l); - -void detectorXYToGvecArray_cfunc(long int npts, double * xy, - double * rMat_d, double * rMat_s, - double * tVec_d, double * tVec_s, double * tVec_c, - double * beamVec, double * etaVec, - double * tTh, double * eta, double * gVec_l); - -void oscillAnglesOfHKLs_cfunc(long int npts, double * hkls, double chi, - double * rMat_c, double * bMat, double wavelength, - double * vInv_s, double * beamVec, double * etaVec, - double * oangs0, double * oangs1); - -/******************************************************************************/ -/* Utility Funtions */ - -void unitRowVector_cfunc(int n, double * cIn, double * cOut); - -void unitRowVectors_cfunc(int m, int n, double * cIn, double * cOut); - -void makeDetectorRotMat_cfunc(double * tPtr, double * rPtr); - -void makeOscillRotMat_cfunc(double chi, double ome, double * rPtr); - -void makeRotMatOfExpMap_cfunc(double * ePtr, double * rPtr); - -void makeRotMatOfQuat_cfunc(int nq, double * qPtr, double * rPtr); - -void makeBinaryRotMat_cfunc(double * aPtr, double * rPtr); - -void makeEtaFrameRotMat_cfunc(double * bPtr, double * ePtr, double * rPtr); - -void validateAngleRanges_cfunc(int na, double * aPtr, int nr, double * minPtr, double * maxPtr, bool * rPtr, int ccw); - -void rotate_vecs_about_axis_cfunc(long int na, double * angles, - long int nax, double * axes, - long int nv, double * vecs, - double * rVecs); - -double quat_distance_cfunc(int nsym, double * q1, double * q2, double * qsym); - -void homochoricOfQuat_cfunc(int nq, double * qPtr, double * hPtr); diff --git a/hexrd/core/transforms/xf.py b/hexrd/core/transforms/xf.py deleted file mode 100644 index 10ae71fc4..000000000 --- a/hexrd/core/transforms/xf.py +++ /dev/null @@ -1,1191 +0,0 @@ -#! /usr/bin/env python -# ============================================================ -# Copyright (c) 2012, Lawrence Livermore National Security, LLC. -# Produced at the Lawrence Livermore National Laboratory. -# Written by Joel Bernier and others. -# LLNL-CODE-529294. -# All rights reserved. -# -# This file is part of HEXRD. For details on dowloading the source, -# see the file COPYING. -# -# Please also see the file LICENSE. -# -# This program is free software; you can redistribute it and/or modify it under -# the terms of the GNU Lesser General Public License (as published by the Free -# Software Foundation) version 2.1 dated February 1999. -# -# This program is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the IMPLIED WARRANTY OF MERCHANTABILITY -# or FITNESS FOR A PARTICULAR PURPOSE. See the terms and conditions of the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU Lesser General Public -# License along with this program (see file LICENSE); if not, write to -# the Free Software Foundation, Inc., 59 Temple Place, Suite 330, -# Boston, MA 02111-1307 USA or visit . -# ============================================================ - -import logging - -import numpy as np -import numba - -# np.seterr(invalid='ignore') - -import scipy.sparse as sparse - -from hexrd.core import matrixutil as mutil - - -# Added to not break people importing these methods -from hexrd.core.rotations import ( - mapAngle, - quatProductMatrix as quat_product_matrix, - arccosSafe, - angularDifference, -) -from hexrd.core.matrixutil import columnNorm, rowNorm - -logger = logging.getLogger(__name__) - -# ============================================================================= -# Module Data -# ============================================================================= - -epsf = np.finfo(float).eps # ~2.2e-16 -ten_epsf = 10 * epsf # ~2.2e-15 -sqrt_epsf = np.sqrt(epsf) # ~1.5e-8 - -periodDict = {'degrees': 360.0, 'radians': 2 * np.pi} -angularUnits = 'radians' # module-level angle units -d2r = np.pi / 180.0 - -# basis vectors -I3 = np.eye(3) # (3, 3) identity -Xl = np.array([[1.0, 0.0, 0.0]], order='C').T # X in the lab frame -Yl = np.array([[0.0, 1.0, 0.0]], order='C').T # Y in the lab frame -Zl = np.array([[0.0, 0.0, 1.0]], order='C').T # Z in the lab frame - -zeroVec = np.zeros(3, order='C') - -# reference beam direction and eta=0 ref in LAB FRAME for standard geometry -bVec_ref = -Zl -eta_ref = Xl - -# reference stretch -vInv_ref = np.array([[1.0, 1.0, 1.0, 0.0, 0.0, 0.0]], order='C').T - - -# ============================================================================= -# Functions -# ============================================================================= - - -def makeGVector(hkl, bMat): - """ - take a CRYSTAL RELATIVE B matrix onto a list of hkls to output unit - reciprocal lattice vectors (a.k.a. lattice plane normals) - - Required Arguments: - hkls -- (3, n) ndarray of n hstacked reciprocal lattice vector component - triplets - bMat -- (3, 3) ndarray representing the matirix taking reciprocal lattice - vectors to the crystal reference frame - - Output: - gVecs -- (3, n) ndarray of n unit reciprocal lattice vectors - (a.k.a. lattice plane normals) - - To Do: - * might benefit from some assert statements to catch improperly shaped - input. - """ - assert hkl.shape[0] == 3, 'hkl input must be (3, n)' - return unitVector(np.dot(bMat, hkl)) - - -@numba.njit(nogil=True, cache=True) -def _anglesToGVecHelper(angs, out): - # gVec_e = np.vstack([[np.cos(0.5*angs[:, 0]) * np.cos(angs[:, 1])], - # [np.cos(0.5*angs[:, 0]) * np.sin(angs[:, 1])], - # [np.sin(0.5*angs[:, 0])]]) - n = angs.shape[0] - for i in range(n): - ca0 = np.cos(0.5 * angs[i, 0]) - sa0 = np.sin(0.5 * angs[i, 0]) - ca1 = np.cos(angs[i, 1]) - sa1 = np.sin(angs[i, 1]) - out[i, 0] = ca0 * ca1 - out[i, 1] = ca0 * sa1 - out[i, 2] = sa0 - - -def anglesToGVec(angs, bHat_l, eHat_l, rMat_s=I3, rMat_c=I3): - """ - from 'eta' frame out to lab - (with handy kwargs to go to crystal or sample) - """ - rMat_e = makeEtaFrameRotMat(bHat_l, eHat_l) - gVec_e = np.empty((angs.shape[0], 3)) - _anglesToGVecHelper(angs, gVec_e) - mat = np.dot(rMat_c.T, np.dot(rMat_s.T, rMat_e)) - return np.dot(mat, gVec_e.T) - - -def gvecToDetectorXY( - gVec_c, rMat_d, rMat_s, rMat_c, tVec_d, tVec_s, tVec_c, beamVec=bVec_ref -): - """ - Takes a list of unit reciprocal lattice vectors in crystal frame to the - specified detector-relative frame. - - The following conditions must be satisfied: - - 1) the reciprocal lattice vector must be able to satisfy a bragg condition - 2) the associated diffracted beam must intersect the detector plane - - Parameters - ---------- - gVec_c : numpy.ndarray - (3, n) array of n reciprocal lattice vectors in the CRYSTAL FRAME. - rMat_d : numpy.ndarray - (3, 3) array, the COB taking DETECTOR FRAME components to LAB FRAME. - rMat_s : numpy.ndarray - (3, 3) array, the COB taking SAMPLE FRAME components to LAB FRAME. - rMat_c : numpy.ndarray - (3, 3) array, the COB taking CRYSTAL FRAME components to SAMPLE FRAME. - tVec_d : numpy.ndarray - (3, 1) array, the translation vector connecting LAB to DETECTOR. - tVec_s : numpy.ndarray - (3, 1) array, the translation vector connecting LAB to SAMPLE. - tVec_c : numpy.ndarray - (3, 1) array, the translation vector connecting SAMPLE to CRYSTAL. - beamVec : numpy.ndarray, optional - (3, 1) array, the vector pointing to the X-ray source in the LAB FRAME. - The default is bVec_ref. - - Returns - ------- - numpy.ndarray - (3, m)darray containing the intersections of m <= n diffracted beams - associated with gVecs. - - """ - ztol = epsf - - nVec_l = np.dot(rMat_d, Zl) # detector plane normal - bHat_l = unitVector(beamVec.reshape(3, 1)) # make sure beam vector is unit - P0_l = tVec_s + np.dot(rMat_s, tVec_c) # origin of CRYSTAL FRAME - P3_l = tVec_d # origin of DETECTOR FRAME - - # form unit reciprocal lattice vectors in lab frame (w/o translation) - gVec_l = np.dot(rMat_s, np.dot(rMat_c, unitVector(gVec_c))) - - # dot with beam vector (upstream direction) - bDot = np.dot(-bHat_l.T, gVec_l).squeeze() - - # see who can diffract; initialize output array with NaNs - canDiffract = np.atleast_1d(np.logical_and(bDot >= ztol, bDot <= 1.0 - ztol)) - npts = sum(canDiffract) - retval = np.nan * np.ones_like(gVec_l) - if np.any(canDiffract): - # subset of admissable reciprocal lattice vectors - adm_gVec_l = gVec_l[:, canDiffract].reshape(3, npts) - - # initialize diffracted beam vector array - dVec_l = np.empty((3, npts)) - for ipt in range(npts): - dVec_l[:, ipt] = np.dot( - makeBinaryRotMat(adm_gVec_l[:, ipt]), -bHat_l - ).squeeze() - - # ############################################################### - # displacement vector calculation - - # first check for non-instersections - denom = np.dot(nVec_l.T, dVec_l).flatten() - dzero = abs(denom) < ztol - denom[dzero] = 1.0 # mitigate divide-by-zero - cantIntersect = denom > 0.0 # index to dVec_l that can't hit det - - # displacement scaling (along dVec_l) - u = np.dot(nVec_l.T, P3_l - P0_l).flatten() / denom - - # filter out non-intersections, fill with NaNs - u[np.logical_or(dzero, cantIntersect)] = np.nan - - # diffracted beam points IN DETECTOR FRAME - P2_l = P0_l + np.tile(u, (3, 1)) * dVec_l - P2_d = np.dot(rMat_d.T, P2_l - tVec_d) - - # put feasible transformed gVecs into return array - retval[:, canDiffract] = P2_d - return retval[:2, :].T - - -def detectorXYToGvec( - xy_det, - rMat_d, - rMat_s, - tVec_d, - tVec_s, - tVec_c, - distortion=None, - beamVec=bVec_ref, - etaVec=eta_ref, - output_ref=False, -): - """ - Takes a list cartesian (x, y) pairs in the detector coordinates and - calculates the associated reciprocal lattice (G) vectors and - (bragg angle, azimuth) pairs with respect to the specified beam and - azimth (eta) reference directions. - - Parameters - ---------- - xy_det : numpy.ndarray - DESCRIPTION. - rMat_d : numpy.ndarray - (3, 3) array, the COB taking DETECTOR FRAME components to LAB FRAME. - rMat_s : numpy.ndarray - (3, 3) array, the COB taking SAMPLE FRAME components to LAB FRAME. - tVec_d : numpy.ndarray - (3, 1) array, the translation vector connecting LAB to DETECTOR. - tVec_s : numpy.ndarray - (3, 1) array, the translation vector connecting LAB to SAMPLE. - tVec_c : numpy.ndarray - (3, 1) array, the translation vector connecting SAMPLE to CRYSTAL. - distortion : class, optional - the distortion operator, if applicable. The default is None. - beamVec : numpy.ndarray, optional - (3, 1) array, the vector pointing to the X-ray source in the LAB FRAME. - The default is bVec_ref. - etaVec : numpy.ndarray, optional - (3, 1) array, the vector defining eta=0 in the LAB FRAME. - The default is eta_ref. - output_ref : bool, optional - if true, will output calculated angles in the LAB FRAME. - The default is False. - - Returns - ------- - tth_eta_ref : tuple, optional - if output_ref is True, (n, 2) ndarray containing the (tTh, eta) pairs - in the LAB REFERENCE FRAME associated with each (x, y) - tth_eta : tuple - (n, 2) ndarray containing the (tTh, eta) pairs associated with - each (x, y). - gVec_l : numpy.ndarray - (3, n) ndarray containing the associated G vector directions in the - LAB FRAME associated with gVecs - """ - - npts = len(xy_det) # number of input (x, y) pairs - - # force unit vectors for beam vector and ref eta - bHat_l = unitVector(beamVec.reshape(3, 1)) - eHat_l = unitVector(etaVec.reshape(3, 1)) - - if distortion is not None: - xy_det = distortion.apply(xy_det) - - # form in-plane vectors for detector points list in DETECTOR FRAME - P2_d = np.hstack([np.atleast_2d(xy_det), np.zeros((npts, 1))]).T - - # in LAB FRAME - P2_l = np.dot(rMat_d, P2_d) + tVec_d - P0_l = tVec_s + np.dot(rMat_s, tVec_c) # origin of CRYSTAL FRAME - - # diffraction unit vector components in LAB FRAME - dHat_l = unitVector(P2_l - P0_l) - - # ############################################################### - # generate output - - # DEBUGGING - assert ( - abs(np.dot(bHat_l.T, eHat_l)) < 1.0 - sqrt_epsf - ), "eta ref and beam cannot be parallel!" - - rMat_e = makeEtaFrameRotMat(bHat_l, eHat_l) - dHat_e = np.dot(rMat_e.T, dHat_l) - - tTh = np.arccos(np.dot(bHat_l.T, dHat_l)).flatten() - eta = np.arctan2(dHat_e[1, :], dHat_e[0, :]).flatten() - - # angles for reference frame - dHat_ref_l = unitVector(P2_l) - dHat_ref_e = np.dot(rMat_e.T, dHat_ref_l) - tTh_ref = np.arccos(np.dot(bHat_l.T, unitVector(P2_l))).flatten() - eta_ref = np.arctan2(dHat_ref_e[1, :], dHat_ref_e[0, :]).flatten() - - # get G-vectors by rotating d by 90-theta about b x d - # (numpy 'cross' works on row vectors) - n_g = unitVector(np.cross(bHat_l.T, dHat_l.T).T) - - gVec_l = rotate_vecs_about_axis(0.5 * (np.pi - tTh), n_g, dHat_l) - - if output_ref: - return (tTh_ref, eta_ref), (tTh, eta), gVec_l - return (tTh, eta), gVec_l - - -def oscillAnglesOfHKLs( - hkls, - chi, - rMat_c, - bMat, - wavelength, - vInv=vInv_ref, - beamVec=bVec_ref, - etaVec=eta_ref, -): - """ - - - Parameters - ---------- - hkls : numpy.ndarray - (3, n_) array of reciprocal lattice vector components. - chi : scalar - the inclination angle of the SAMPLE FRAME about the LAB X. - rMat_c : numpy.ndarray - (3, 3) array, the COB taking CRYSTAL FRAME components to SAMPLE FRAME. - bMat : numpy.ndarray - (3, 3) array, the COB matrix taking RECIPROCAL FRAME components to the - CRYSTAL FRAME. - wavelength : scalar - the X-ray wavelength in Ångstroms. - vInv : numpy.ndarray, optional - (6, 1) array representing the components of the inverse stretch tensor - in the SAMPLE FRAME. The default is vInv_ref. - beamVec : numpy.ndarray, optional - (3, 1) array, the vector pointing to the X-ray source in the LAB FRAME. - The default is bVec_ref. - etaVec : numpy.ndarray, optional - (3, 1) array, the vector defining eta=0 in the LAB FRAME. - The default is eta_ref. - - Returns - ------- - retval : tuple - Two (3, n) numpy.ndarrays representing the pair of feasible - (tTh, eta, ome) solutions for the n input hkls. - - Notes - ----- - The reciprocal lattice vector, G, will satisfy the the Bragg condition - when: - - b.T * G / ||G|| = -sin(theta) - - where b is the incident beam direction (k_i) and theta is the Bragg - angle consistent with G and the specified wavelength. The components of - G in the lab frame in this case are obtained using the crystal - orientation, Rc, and the single-parameter oscillation matrix, Rs(ome): - - Rs(ome) * Rc * G / ||G|| - - The equation above can be rearranged to yeild an expression of the form: - - a*sin(ome) + b*cos(ome) = c - - which is solved using the relation: - - a*sin(x) + b*cos(x) = sqrt(a**2 + b**2) * sin(x + alpha) - - --> sin(x + alpha) = c / sqrt(a**2 + b**2) - - where: - - alpha = arctan2(b, a) - - The solutions are: - - / - | arcsin(c / sqrt(a**2 + b**2)) - alpha - x = < - | pi - arcsin(c / sqrt(a**2 + b**2)) - alpha - \\ - - There is a double root in the case the reflection is tangent to the - Debye-Scherrer cone (c**2 = a**2 + b**2), and no solution if the - Laue condition cannot be satisfied (filled with NaNs in the results - array here) - """ - # reciprocal lattice vectors in CRYSTAL frame - gVec_c = np.dot(bMat, hkls) - # stretch tensor in SAMPLE frame - vMat_s = mutil.vecMVToSymm(vInv) - # reciprocal lattice vectors in SAMPLE frame - gVec_s = np.dot(vMat_s, np.dot(rMat_c, gVec_c)) - # unit reciprocal lattice vectors in SAMPLE frame - gHat_s = unitVector(gVec_s) - # unit reciprocal lattice vectors in CRYSTAL frame - gHat_c = np.dot(rMat_c.T, gHat_s) - # make sure beam direction is a unit vector - bHat_l = unitVector(beamVec.reshape(3, 1)) - # make sure eta = 0 direction is a unit vector - eHat_l = unitVector(etaVec.reshape(3, 1)) - - # sin of the Bragg angle assoc. with wavelength - sintht = 0.5 * wavelength * columnNorm(gVec_s) - - # sin and cos of the oscillation axis tilt - cchi = np.cos(chi) - schi = np.sin(chi) - - # coefficients for harmonic equation - a = ( - gHat_s[2, :] * bHat_l[0] - + schi * gHat_s[0, :] * bHat_l[1] - - cchi * gHat_s[0, :] * bHat_l[2] - ) - b = ( - gHat_s[0, :] * bHat_l[0] - - schi * gHat_s[2, :] * bHat_l[1] - + cchi * gHat_s[2, :] * bHat_l[2] - ) - c = -sintht - cchi * gHat_s[1, :] * bHat_l[1] - schi * gHat_s[1, :] * bHat_l[2] - - # should all be 1-d: a = a.flatten(); b = b.flatten(); c = c.flatten() - - # form solution - abMag = np.sqrt(a * a + b * b) - assert np.all(abMag > 0), "Beam vector specification is infealible!" - phaseAng = np.arctan2(b, a) - rhs = c / abMag - rhs[abs(rhs) > 1.0] = np.nan - rhsAng = np.arcsin(rhs) # will give NaN for abs(rhs) > 1. + 0.5*epsf - - # write ome angle output arrays (NaNs persist here) - ome0 = rhsAng - phaseAng - ome1 = np.pi - rhsAng - phaseAng - - goodOnes_s = -np.isnan(ome0) - - # DEBUGGING - assert np.all( - np.isnan(ome0) == np.isnan(ome1) - ), "infeasible hkls do not match for ome0, ome1!" - - # do etas -- ONLY COMPUTE IN CASE CONSISTENT REFERENCE COORDINATES - if abs(np.dot(bHat_l.T, eHat_l)) < 1.0 - sqrt_epsf and np.any(goodOnes_s): - eta0 = np.nan * np.ones_like(ome0) - eta1 = np.nan * np.ones_like(ome1) - - # make eta basis COB with beam antiparallel with Z - rMat_e = makeEtaFrameRotMat(bHat_l, eHat_l) - - goodOnes = np.tile(goodOnes_s, (1, 2)).flatten() - - numGood_s = sum(goodOnes_s) - numGood = 2 * numGood_s - tmp_eta = np.empty(numGood) - tmp_gvec = np.tile(gHat_c, (1, 2))[:, goodOnes] - allome = np.hstack([ome0, ome1]) - - for i in range(numGood): - rMat_s = makeOscillRotMat([chi, allome[goodOnes][i]]) - gVec_e = np.dot( - rMat_e.T, - np.dot(rMat_s, np.dot(rMat_c, tmp_gvec[:, i].reshape(3, 1))), - ) - tmp_eta[i] = np.arctan2(gVec_e[1], gVec_e[0]) - eta0[goodOnes_s] = tmp_eta[:numGood_s] - eta1[goodOnes_s] = tmp_eta[numGood_s:] - - # make assoc tTh array - tTh = 2.0 * np.arcsin(sintht).flatten() - tTh0 = tTh - tTh0[-goodOnes_s] = np.nan - retval = ( - np.vstack([tTh0.flatten(), eta0.flatten(), ome0.flatten()]), - np.vstack([tTh0.flatten(), eta1.flatten(), ome1.flatten()]), - ) - else: - retval = (ome0.flatten(), ome1.flatten()) - return retval - - -def polarRebin( - thisFrame, - npdiv=2, - mmPerPixel=(0.2, 0.2), - convertToTTh=False, - rMat_d=I3, - tVec_d=np.r_[0.0, 0.0, -1000.0], - beamVec=bVec_ref, - etaVec=eta_ref, - rhoRange=np.r_[20, 200], - numRho=1000, - etaRange=(d2r * np.r_[-5, 355]), - numEta=36, - verbose=True, - log=None, -): - """ - Performs polar rebinning of an input image. - - Parameters - ---------- - thisFrame : array_like - An (n, m) array representing an image. - npdiv : int, optional - The number of pixel subdivision. The default is 2. - mmPerPixel : array_like, optional - A 2-element sequence with the pixel pitch in mm for the row and column - dimensions, respectively. The default is (0.2, 0.2). - convertToTTh : bool, optional - If true, output abcissa is in angular (two-theta) coordinates. - The default is False. - rMat_d : numpy.ndarray, optional - (3, 3) array, the COB taking DETECTOR FRAME components to LAB FRAME. - The default is I3. - tVec_d : numpy.ndarray, optional - (3, 1) array, the translation vector connecting LAB to DETECTOR. - The default is np.r_[0., 0., -1000.]. - beamVec : TYPE, optional - DESCRIPTION. The default is bVec_ref. - etaVec : TYPE, optional - DESCRIPTION. The default is eta_ref. - rhoRange : array_like, optional - the min/max limits of the rebinning region in pixels. - The default is np.r_[20, 200]. - numRho : int, optional - the number of bins to use in the radial dimension. The default is 1000. - etaRange : array_like, optional - A 2-element sequence describing the min/max azimuthal angles in - radians. The default is (-0.08726646259971647, 6.19591884457987). - numEta : int, optional - The number of azimuthal bins. The default is 36. - verbose : bool, optional - DESCRIPTION. The default is True. - log : TYPE, optional - DESCRIPTION. The default is None. - - Raises - ------ - RuntimeError - DESCRIPTION. - - Returns - ------- - polImg : TYPE - DESCRIPTION. - - """ - """ - Caking algorithm - - INPUTS - - thisFrame - npdiv=2, pixel subdivision (n x n) to determine bin membership - rhoRange=[100, 1000] - radial range in pixels - numRho=1200 - number of radial bins - etaRange=np.pi*np.r_[-5, 355]/180. -- range of eta - numEta=36 - number of eta subdivisions - ROI=None - region of interest (four vector) - corrected=False - uses 2-theta instead of rho - verbose=True, - - """ - - startEta = etaRange[0] - stopEta = etaRange[1] - - startRho = rhoRange[0] - stopRho = rhoRange[1] - - subPixArea = 1 / float(npdiv) ** 2 # areal rescaling for subpixel intensities - - # MASTER COORDINATES - # - in pixel indices, UPPER LEFT PIXEL is [0, 0] --> (row, col) - # - in fractional pixels, UPPER LEFT CORNER is [-0.5, -0.5] --> (row, col) - # - in cartesian frame, the LOWER LEFT CORNER is [0, 0] --> (col, row) - x = thisFrame[0, :, :].flatten() - y = thisFrame[1, :, :].flatten() - roiData = thisFrame[2, :, :].flatten() - - # need rhos (or tThs) and etas) - if convertToTTh: - dAngs = detectorXYToGvec( - np.vstack([x, y]).T, - rMat_d, - I3, - tVec_d, - zeroVec, - zeroVec, - beamVec=beamVec, - etaVec=etaVec, - ) - rho = dAngs[0][0] # this is tTh now - eta = dAngs[0][1] - else: - # in here, we are vanilla cartesian - rho = np.sqrt(x * x + y * y) - eta = np.arctan2(y, x) - eta = mapAngle(eta, [startEta, 2 * np.pi + startEta], units='radians') - - # MAKE POLAR BIN CENTER ARRAY - deltaEta = (stopEta - startEta) / float(numEta) - deltaRho = (stopRho - startRho) / float(numRho) - - rowEta = startEta + deltaEta * (np.arange(numEta) + 0.5) - colRho = startRho + deltaRho * (np.arange(numRho) + 0.5) - - # initialize output dictionary - polImg = {} - polImg['radius'] = colRho - polImg['azimuth'] = rowEta - polImg['intensity'] = np.zeros((numEta, numRho)) - polImg['deltaRho'] = deltaRho - - if verbose: - msg = "INFO: Masking pixels\n" - if log: - log.write(msg) - else: - logger.info(msg) - - rhoI = startRho - 10 * deltaRho - rhoF = stopRho + 10 * deltaRho - inAnnulus = np.where((rho >= rhoI) & (rho <= rhoF))[0] - for i in range(numEta): - if verbose: - msg = f"INFO: Processing sector {i + 1} of {numEta}\n" - if log: - log.write(msg) - else: - logger.info(msg) - - # import pdb;pdb.set_trace() - etaI1 = rowEta[i] - 10.5 * deltaEta - etaF1 = rowEta[i] + 10.5 * deltaEta - - tmpEta = eta[inAnnulus] - inSector = np.where((tmpEta >= etaI1) & (tmpEta <= etaF1))[0] - - nptsIn = len(inSector) - - tmpX = x[inAnnulus[inSector]] - tmpY = y[inAnnulus[inSector]] - tmpI = roiData[inAnnulus[inSector]] - - # subdivide pixels - # - note that these are in fractional pixel coordinates (centered) - # - must convert to working units (see 'self.pixelPitchUnits') - subCrds = (np.arange(npdiv) + 0.5) / npdiv - - intX, intY = np.meshgrid(subCrds, subCrds) - - intX = np.tile(intX.flatten(), (nptsIn, 1)).T.flatten() - intY = np.tile(intY.flatten(), (nptsIn, 1)).T.flatten() - - # expand coords using pixel subdivision - tmpX = np.tile(tmpX, (npdiv**2, 1)).flatten() + (intX - 0.5) * mmPerPixel[0] - tmpY = np.tile(tmpY, (npdiv**2, 1)).flatten() + (intY - 0.5) * mmPerPixel[1] - tmpI = np.tile(tmpI, (npdiv**2, 1)).flatten() / subPixArea - - if convertToTTh: - dAngs = detectorXYToGvec( - np.vstack([tmpX, tmpY]).T, - rMat_d, - I3, - tVec_d, - zeroVec, - zeroVec, - beamVec=beamVec, - etaVec=etaVec, - ) - tmpRho = dAngs[0][0] # this is tTh now - tmpEta = dAngs[0][1] - else: - tmpRho = np.sqrt(tmpX * tmpX + tmpY * tmpY) - tmpEta = np.arctan2(tmpY, tmpX) - tmpEta = mapAngle(tmpEta, [startEta, 2 * np.pi + startEta], units='radians') - - etaI2 = rowEta[i] - 0.5 * deltaEta - etaF2 = rowEta[i] + 0.5 * deltaEta - - inSector2 = ((tmpRho >= startRho) & (tmpRho <= stopRho)) & ( - (tmpEta >= etaI2) & (tmpEta <= etaF2) - ) - - tmpRho = tmpRho[inSector2] - tmpI = tmpI[inSector2] - - binId = np.floor((tmpRho - startRho) / deltaRho) - nSubpixelsIn = len(binId) - - if nSubpixelsIn > 0: - tmpI = sparse.csc_matrix( - (tmpI, (binId, np.arange(nSubpixelsIn))), - shape=(numRho, nSubpixelsIn), - ) - binId = sparse.csc_matrix( - (np.ones(nSubpixelsIn), (binId, np.arange(nSubpixelsIn))), - shape=(numRho, nSubpixelsIn), - ) - - # Normalized contribution to the ith sector's radial bins - binIdSum = np.asarray(binId.sum(1)).flatten() - whereNZ = np.asarray(np.not_equal(polImg['intensity'][i, :], binIdSum)) - polImg['intensity'][i, whereNZ] = ( - np.asarray(tmpI.sum(1))[whereNZ].flatten() / binIdSum[whereNZ] - ) - return polImg - - -# ============================================================================= -# Utility Functions -# ============================================================================= - - -def _z_project(x, y): - return np.cos(x) * np.sin(y) - np.sin(x) * np.cos(y) - - -def reg_grid_indices(edges, points_1d): - """ - get indices in a 1-d regular grid. - - edges are just that: - - point: x (2.5) - | - edges: |1 |2 |3 |4 |5 - ------------------------- - indices: | 0 | 1 | 2 | 3 | - ------------------------- - - above the deltas are + and the index for the point is 1 - - point: x (2.5) - | - edges: |5 |4 |3 |2 |1 - ------------------------- - indices: | 0 | 1 | 2 | 3 | - ------------------------- - - here the deltas are - and the index for the point is 2 - - * can handle grids with +/- deltas - * be careful when using with a cyclical angular array! edges and points - must be mapped to the same branch cut, and - abs(edges[0] - edges[-1]) = 2*pi - """ - ztol = 1e-12 - - assert len(edges) >= 2, "must have at least 2 edges" - - points_1d = np.r_[points_1d].flatten() - delta = float(edges[1] - edges[0]) - - on_last_rhs = points_1d >= edges[-1] - ztol - points_1d[on_last_rhs] = points_1d[on_last_rhs] - ztol - - if delta > 0: - on_last_rhs = points_1d >= edges[-1] - ztol - points_1d[on_last_rhs] = points_1d[on_last_rhs] - ztol - idx = np.floor((points_1d - edges[0]) / delta) - elif delta < 0: - on_last_rhs = points_1d <= edges[-1] + ztol - points_1d[on_last_rhs] = points_1d[on_last_rhs] + ztol - idx = np.ceil((points_1d - edges[0]) / delta) - 1 - else: - raise RuntimeError("edges array gives delta of 0") - return np.array(idx, dtype=int) - - -@numba.njit(nogil=True, cache=True) -def _unitVectorSingle(a, b): - n = a.shape[0] - nrm = 0.0 - for i in range(n): - nrm += a[i] * a[i] - nrm = np.sqrt(nrm) - # prevent divide by zero - if nrm > epsf: - for i in range(n): - b[i] = a[i] / nrm - else: - for i in range(n): - b[i] = a[i] - - -@numba.njit(nogil=True, cache=True) -def _unitVectorMulti(a, b): - n = a.shape[0] - m = a.shape[1] - for j in range(m): - nrm = 0.0 - for i in range(n): - nrm += a[i, j] * a[i, j] - nrm = np.sqrt(nrm) - # prevent divide by zero - if nrm > epsf: - for i in range(n): - b[i, j] = a[i, j] / nrm - else: - for i in range(n): - b[i, j] = a[i, j] - - -def unitVector(a): - """ - normalize array of column vectors (hstacked, axis = 0) - """ - result = np.empty_like(a) - if a.ndim == 1: - _unitVectorSingle(a, result) - elif a.ndim == 2: - _unitVectorMulti(a, result) - else: - raise ValueError( - "incorrect arg shape; must be 1-d or 2-d, " + "yours is %d-d" % (a.ndim) - ) - return result - - -def makeDetectorRotMat(tiltAngles): - """ - Form the (3, 3) tilt rotations from the tilt angle list: - - tiltAngles = [gamma_Xl, gamma_Yl, gamma_Zl] in radians - """ - # form rMat_d from parameter list - cos_gX = np.cos(tiltAngles[0]) - sin_gX = np.sin(tiltAngles[0]) - - cos_gY = np.cos(tiltAngles[1]) - sin_gY = np.sin(tiltAngles[1]) - - cos_gZ = np.cos(tiltAngles[2]) - sin_gZ = np.sin(tiltAngles[2]) - - rotXl = np.array([[1.0, 0.0, 0.0], [0.0, cos_gX, -sin_gX], [0.0, sin_gX, cos_gX]]) - rotYl = np.array([[cos_gY, 0.0, sin_gY], [0.0, 1.0, 0.0], [-sin_gY, 0.0, cos_gY]]) - rotZl = np.array([[cos_gZ, -sin_gZ, 0.0], [sin_gZ, cos_gZ, 0.0], [0.0, 0.0, 1.0]]) - return np.dot(rotZl, np.dot(rotYl, rotXl)) - - -def makeOscillRotMat(oscillAngles): - """ - oscillAngles = [chi, ome] - """ - cchi = np.cos(oscillAngles[0]) - schi = np.sin(oscillAngles[0]) - - come = np.cos(oscillAngles[1]) - some = np.sin(oscillAngles[1]) - - rchi = np.array([[1.0, 0.0, 0.0], [0.0, cchi, -schi], [0.0, schi, cchi]]) - - rome = np.array([[come, 0.0, some], [0.0, 1.0, 0.0], [-some, 0.0, come]]) - - return np.dot(rchi, rome) - - -def makeRotMatOfExpMap(expMap): - """ - Calculate the rotation matrix of an exponential map. - - Parameters - ---------- - expMap : array_like - A 3-element sequence representing the exponential map n*phi. - - Returns - ------- - rMat : np.ndarray - The (3, 3) array representing the rotation matrix of the input - exponential map parameters. - """ - expMap = np.asarray(expMap).flatten() - phi = np.norm(expMap) - if phi > epsf: - wMat = np.array( - [ - [0.0, -expMap[2], expMap[1]], - [expMap[2], 0.0, -expMap[0]], - [-expMap[1], expMap[0], 0.0], - ] - ) - - rMat = ( - I3 - + (np.sin(phi) / phi) * wMat - + ((1.0 - np.cos(phi)) / (phi * phi)) * np.dot(wMat, wMat) - ) - else: - rMat = I3 - - return rMat - - -def makeBinaryRotMat(axis): - """ """ - n = np.asarray(axis).flatten() - assert len(n) == 3, 'Axis input does not have 3 components' - return 2 * np.dot(n.reshape(3, 1), n.reshape(1, 3)) - I3 - - -@numba.njit(nogil=True, cache=True) -def _makeEtaFrameRotMat(bHat_l, eHat_l, out): - # bHat_l and eHat_l CANNOT have 0 magnitude! - # must catch this case as well as colinear bHat_l/eHat_l elsewhere... - bHat_mag = np.sqrt(bHat_l[0] ** 2 + bHat_l[1] ** 2 + bHat_l[2] ** 2) - - # assign Ze as -bHat_l - for i in range(3): - out[i, 2] = -bHat_l[i] / bHat_mag - - # find Ye as Ze ^ eHat_l - Ye0 = out[1, 2] * eHat_l[2] - eHat_l[1] * out[2, 2] - Ye1 = out[2, 2] * eHat_l[0] - eHat_l[2] * out[0, 2] - Ye2 = out[0, 2] * eHat_l[1] - eHat_l[0] * out[1, 2] - - Ye_mag = np.sqrt(Ye0**2 + Ye1**2 + Ye2**2) - - out[0, 1] = Ye0 / Ye_mag - out[1, 1] = Ye1 / Ye_mag - out[2, 1] = Ye2 / Ye_mag - - # find Xe as Ye ^ Ze - out[0, 0] = out[1, 1] * out[2, 2] - out[1, 2] * out[2, 1] - out[1, 0] = out[2, 1] * out[0, 2] - out[2, 2] * out[0, 1] - out[2, 0] = out[0, 1] * out[1, 2] - out[0, 2] * out[1, 1] - - -def makeEtaFrameRotMat(bHat_l, eHat_l): - """ - make eta basis COB matrix with beam antiparallel with Z - - takes components from ETA frame to LAB - - **NO EXCEPTION HANDLING FOR COLINEAR ARGS IN NUMBA VERSION! - - ...put checks for non-zero magnitudes and non-colinearity in wrapper? - """ - result = np.empty((3, 3)) - _makeEtaFrameRotMat(bHat_l.reshape(3), eHat_l.reshape(3), result) - return result - - -def angles_in_range(angles, starts, stops, degrees=True): - """Determine whether angles lie in or out of specified ranges - - *angles* - a list/array of angles - *starts* - a list of range starts - *stops* - a list of range stops - - OPTIONAL ARGS: - *degrees* - [True] angles & ranges in degrees (or radians) - """ - TAU = 360.0 if degrees else 2 * np.pi - nw = len(starts) - na = len(angles) - in_range = np.zeros((na), dtype=bool) - for i in range(nw): - amin = starts[i] - amax = stops[i] - for j in range(na): - a = angles[j] - acheck = amin + np.mod(a - amin, TAU) - if acheck <= amax: - in_range[j] = True - - return in_range - - -def validateAngleRanges(angList, startAngs, stopAngs, ccw=True): - """ - A better way to go. find out if an angle is in the range - CCW or CW from start to stop - - There is, of course, an ambigutiy if the start and stop angle are - the same; we treat them as implying 2*pi having been mapped - """ - # Prefer ravel over flatten because flatten never skips the copy - angList = np.asarray(angList).ravel() # needs to have len - startAngs = np.asarray(startAngs).ravel() # needs to have len - stopAngs = np.asarray(stopAngs).ravel() # needs to have len - - n_ranges = len(startAngs) - assert len(stopAngs) == n_ranges, "length of min and max angular limits must match!" - - # to avoid warnings in >=, <= later down, mark nans; - # need these to trick output to False in the case of nan input - nan_mask = np.isnan(angList) - - reflInRange = np.zeros(angList.shape, dtype=bool) - - # bin length for chunking - binLen = np.pi / 2.0 - - # in plane vectors defining wedges - x0 = np.vstack([np.cos(startAngs), np.sin(startAngs)]) - x1 = np.vstack([np.cos(stopAngs), np.sin(stopAngs)]) - - # dot products - dp = np.sum(x0 * x1, axis=0) - if np.any(dp >= 1.0 - sqrt_epsf) and n_ranges > 1: - # ambiguous case - raise RuntimeError("At least one of your ranges is alread 360 degrees!") - elif dp[0] >= 1.0 - sqrt_epsf and n_ranges == 1: - # trivial case! - reflInRange = np.ones(angList.shape, dtype=bool) - reflInRange[nan_mask] = False - else: - # solve for arc lengths - # ...note: no zeros should have made it here - a = x0[0, :] * x1[1, :] - x0[1, :] * x1[0, :] - b = x0[0, :] * x1[0, :] + x0[1, :] * x1[1, :] - phi = np.arctan2(b, a) - - arclen = 0.5 * np.pi - phi # these are clockwise - cw_phis = arclen < 0 - arclen[cw_phis] = 2 * np.pi + arclen[cw_phis] # all positive (CW) now - if not ccw: - arclen = 2 * np.pi - arclen - - if sum(arclen) > 2 * np.pi: - raise RuntimeWarning("Specified angle ranges sum to > 360 degrees") - - # check that there are no more thandp = np.zeros(n_ranges) - for i in range(n_ranges): - # number or subranges using 'binLen' - numSubranges = int(np.ceil(arclen[i] / binLen)) - - # check remaider - binrem = np.remainder(arclen[i], binLen) - if binrem == 0: - finalBinLen = binLen - else: - finalBinLen = binrem - - # if clockwise, negate bin length - if not ccw: - binLen = -binLen - finalBinLen = -finalBinLen - - # Create sub ranges on the fly to avoid ambiguity in dot product - # for wedges >= 180 degrees - subRanges = np.array( - [startAngs[i] + binLen * j for j in range(numSubranges)] - + [startAngs[i] + binLen * (numSubranges - 1) + finalBinLen] - ) - - for k in range(numSubranges): - zStart = _z_project(angList, subRanges[k]) - zStop = _z_project(angList, subRanges[k + 1]) - if ccw: - zStart[nan_mask] = 999.0 - zStop[nan_mask] = -999.0 - reflInRange = reflInRange | np.logical_and(zStart <= 0, zStop >= 0) - else: - zStart[nan_mask] = -999.0 - zStop[nan_mask] = 999.0 - reflInRange = reflInRange | np.logical_and(zStart >= 0, zStop <= 0) - return reflInRange - - -def rotate_vecs_about_axis(angle, axis, vecs): - """ - Rotate vectors about an axis - - INPUTS - *angle* - array of angles (len == 1 or n) - *axis* - array of unit vectors (shape == (3, 1) or (3, n)) - *vec* - array of vectors to be rotated (shape = (3, n)) - - Quaternion formula: - if we split v into parallel and perpedicular components w.r.t. the - axis of quaternion q, - - v = a + n - - then the action of rotating the vector dot(R(q), v) becomes - - v_rot = (q0**2 - |q|**2)(a + n) + 2*dot(q, a)*q + 2*q0*cross(q, n) - - """ - angle = np.atleast_1d(angle) - - # quaternion components - q0 = np.cos(0.5 * angle) - q1 = np.sin(0.5 * angle) - qv = np.tile(q1, (3, 1)) * axis - - # component perpendicular to axes (inherits shape of vecs) - vp0 = ( - vecs[0, :] - - axis[0, :] * axis[0, :] * vecs[0, :] - - axis[0, :] * axis[1, :] * vecs[1, :] - - axis[0, :] * axis[2, :] * vecs[2, :] - ) - vp1 = ( - vecs[1, :] - - axis[1, :] * axis[1, :] * vecs[1, :] - - axis[1, :] * axis[0, :] * vecs[0, :] - - axis[1, :] * axis[2, :] * vecs[2, :] - ) - vp2 = ( - vecs[2, :] - - axis[2, :] * axis[2, :] * vecs[2, :] - - axis[2, :] * axis[0, :] * vecs[0, :] - - axis[2, :] * axis[1, :] * vecs[1, :] - ) - - # dot product with components along; cross product with components normal - qdota = ( - axis[0, :] * vecs[0, :] + axis[1, :] * vecs[1, :] + axis[2, :] * vecs[2, :] - ) * (axis[0, :] * qv[0, :] + axis[1, :] * qv[1, :] + axis[2, :] * qv[2, :]) - qcrossn = np.vstack( - [ - qv[1, :] * vp2 - qv[2, :] * vp1, - qv[2, :] * vp0 - qv[0, :] * vp2, - qv[0, :] * vp1 - qv[1, :] * vp0, - ] - ) - - # quaternion formula - v_rot = ( - np.tile(q0 * q0 - q1 * q1, (3, 1)) * vecs - + 2.0 * np.tile(qdota, (3, 1)) * qv - + 2.0 * np.tile(q0, (3, 1)) * qcrossn - ) - return v_rot - - -def quat_distance(q1, q2, qsym): - """ """ - # qsym from PlaneData objects are (4, nsym) - # convert symmetries to (4, 4) qprod matrices - nsym = qsym.shape[1] - rsym = np.zeros((nsym, 4, 4)) - for i in range(nsym): - rsym[i, :, :] = quat_product_matrix(qsym[:, i], mult='right') - - # inverse of q1 in matrix form - q1i = quat_product_matrix( - np.r_[1, -1, -1, -1] * np.atleast_1d(q1).flatten(), mult='right' - ) - - # Do R * Gc, store as vstacked equivalent quaternions (nsym, 4) - q2s = np.dot(rsym, q2) - - # Calculate the class of misorientations for full symmetrically equivalent - # q1 and q2: (4, ) * (4, nsym) - eqv_mis = np.dot(q1i, q2s.T) - - # find the largest scalar component - q0_max = np.argmax(abs(eqv_mis[0, :])) - - # compute the distance - qmin = eqv_mis[:, q0_max] - - return 2 * arccosSafe(qmin[0] * np.sign(qmin[0])) diff --git a/hexrd/core/transforms/xfcapi.py b/hexrd/core/transforms/xfcapi.py index 1e96661b4..fb7e1cc7f 100644 --- a/hexrd/core/transforms/xfcapi.py +++ b/hexrd/core/transforms/xfcapi.py @@ -1,43 +1,691 @@ -# We will replace these functions with the new versions as we -# add and test them. -# NOTE: we are only importing what is currently being used in hexrd -# and hexrdgui. This is so that we can see clearly what is in use. -from .old_xfcapi import ( - anglesToDVec, - anglesToGVec, - detectorXYToGvec, - gvecToDetectorXY, - gvecToDetectorXYArray, - oscillAnglesOfHKLs, - angularDifference, - makeDetectorRotMat, - makeEtaFrameRotMat, - makeOscillRotMat, - makeOscillRotMatArray, - makeRotMatOfExpMap, - makeRotMatOfQuat, - mapAngle, - rowNorm, - unitRowVector, - bVec_ref, - eta_ref, - Xl, - Yl, -) - - -from .new_capi.xf_new_capi import ( - angles_to_dvec, - angles_to_gvec, - gvec_to_xy, - make_beam_rmat, - make_detector_rmat, - make_rmat_of_expmap, - make_sample_rmat, - oscill_angles_of_hkls, - quat_distance, - rotate_vecs_about_axis, - unit_vector, - validate_angle_ranges, - xy_to_gvec, -) +# -*- coding: utf-8 -*- +""" +Transforms module implementation using a support C extension module. + +Currently, this implementation contains code for the following functions: + + - angles_to_gvec + - angles_to_dvec + - gvec_to_xy + - xy_to_gvec (partial) + + - unit_vector + - make_rmat_of_exp_map + - make_binary_rmat + - make_beam_rmat + - validate_angle_ranges + - rotate_vecs_about_axis + - quat_distance + +There are also some functions that maybe would be needed in the transforms +module: + - makeGVector + - makeRotMatOfQuat + - homochoricOfQuat +""" +from typing import Optional, Tuple, Sequence +import numpy as np +from numpy.typing import NDArray + +from hexrd.core.extensions import transforms_c_api +from hexrd.core.extensions import transforms as cpp_transforms +from hexrd.core.distortion.distortionabc import DistortionABC +from hexrd.core import constants as cnst + + +def angles_to_gvec( + angs: NDArray[np.float64], + beam_vec: NDArray[np.float64] = cnst.beam_vec, + eta_vec: NDArray[np.float64] = cnst.eta_vec, + chi: float = 0.0, + rmat_c: NDArray[np.float64] = cnst.identity_3x3, +) -> NDArray[np.float64]: + """ + + Takes triplets of angles in the beam frame (2*theta, eta[, omega]) + to components of unit G-vectors in the LAB frame. If the omega + values are not trivial (i.e. angs[:, 2] = 0.), then the components + are in the SAMPLE frame. If the crystal rmat is specified and + is not the identity, then the components are in the CRYSTAL frame. + + G vectors here referes to the reciprocal lattice vectors. + + Parameters + ---------- + angs : ndarray + The euler angles of diffraction. The last dimension must be 2 or 3. In + (2*theta, eta[, omega]) format. + beam_vec : ndarray, optional + Unit vector pointing towards the X-ray source in the lab frame. + Defaults to [0,0,-1] + eta_vec : ndarray, optional + Vector defining eta=0 in the lab frame. Defaults to [1,0,0] + chi : float, optional + The inclination angle of the sample frame about the lab frame X. + rmat_c : ndarray, optional + The change of basis matrix from the reciprocal frame to the crystal + frame. Defaults to the identity. + + Returns + ------- + ndarray + (3,n) array of unit reciprocal lattice vectors, frame depends on the + input parameters + """ + + orig_ndim = angs.ndim + + # if only a pair is provided... convert to a triplet with omegas == 0 + # so that behavior is preserved. + if angs.shape[-1] == 2: + angs = np.hstack((angs, np.zeros(angs.shape[:-1] + (1,)))) + + angs = np.ascontiguousarray(np.atleast_2d(angs)) + beam_vec = np.ascontiguousarray(beam_vec.flatten()) + eta_vec = np.ascontiguousarray(eta_vec.flatten()) + rmat_c = np.ascontiguousarray(rmat_c) + + result = cpp_transforms.anglesToGVec(angs, beam_vec, eta_vec, chi, rmat_c) + + return result[0] if orig_ndim == 1 else result + + +def angles_to_dvec( + angs: NDArray[np.float64], + beam_vec: NDArray[np.float64] = cnst.beam_vec, + eta_vec: NDArray[np.float64] = cnst.eta_vec, + chi: float = 0.0, + rmat_c: NDArray[np.float64] = cnst.identity_3x3, +) -> NDArray[np.float64]: + """Calculate diffraction vectors from beam frame angles. + + Takes triplets of angles in the beam frame (2*theta, eta[, omega]) + to components of unit diffraction vectors in the LAB frame. If the omega + values are not trivial (i.e. angs[:, 2] = 0.), then the components + are in the SAMPLE frame. If the crystal rmat is specified and + is not the identity, then the components are in the CRYSTAL frame. + + + Parameters + ---------- + angs : ndarray + The euler angles of diffraction. The last dimension must be 2 or 3. In + (2*theta, eta[, omega]) format. + beam_vec : ndarray, optional + Unit vector pointing towards the X-ray source in the lab frame. + Defaults to [0,0,-1] + eta_vec : ndarray, optional + Vector defining eta=0 in the lab frame. Defaults to [1,0,0] + chi : float, optional + The inclination angle of the sample frame about the lab frame X. + rmat_c : ndarray, optional + The change of basis matrix from the reciprocal frame to the crystal + frame. Defaults to the identity. + + Returns + ------- + ndarray + (3,n) array of unit diffraction vectors, frame depends on the input + parameters + """ + # TODO: Improve capi to avoid multiplications when rmat_c is None + + # if only a pair is provided... convert to a triplet with omegas == 0 + # so that behavior is preserved. + if angs.shape[-1] == 2: + angs = np.hstack((angs, np.zeros(angs.shape[:-1] + (1,)))) + + angs = np.ascontiguousarray(np.atleast_2d(angs)) + beam_vec = np.ascontiguousarray(beam_vec.flatten()) + eta_vec = np.ascontiguousarray(eta_vec.flatten()) + rmat_c = np.ascontiguousarray(rmat_c) + + return cpp_transforms.anglesToDVec(angs, beam_vec, eta_vec, chi, rmat_c) + + +def makeGVector( + hkl: NDArray[np.float64], bMat: NDArray[np.float64] +) -> NDArray[np.float64]: + """ + Take a crystal relative b matrix onto a list of hkls to output unit + reciprocal latice vectors (a.k.a. lattice plane normals) + + + Parameters + ---------- + hkl : ndarray + (3,n) ndarray of n hstacked reciprocal lattice vector component + triplets + bMat : ndarray + (3,3) ndarray of the change of basis matrix from the reciprocal + lattice to the crystal reference frame + + Returns + ------- + ndarray + (3,n) ndarray of n unit reciprocal lattice vectors (a.k.a. lattice + plane normals) + + """ + if hkl.shape[0] != 3: + raise ValueError('hkl input must be (3, n)') + return unit_vector(np.dot(bMat, hkl)) + + +def gvec_to_xy( + gvec_c: NDArray[np.float64], + rmat_d: NDArray[np.float64], + rmat_s: NDArray[np.float64], + rmat_c: NDArray[np.float64], + tvec_d: NDArray[np.float64], + tvec_s: NDArray[np.float64], + tvec_c: NDArray[np.float64], + beam_vec: Optional[NDArray[np.float64]] = None, + vmat_inv: Optional[NDArray[np.float64]] = None, + bmat: Optional[NDArray[np.float64]] = None, +) -> NDArray[np.float64]: + """ + Takes a concatenated list of reciprocal lattice vectors components in the + CRYSTAL FRAME to the specified detector-relative frame, subject to the + following: + + 1) it must be able to satisfy a bragg condition + 2) the associated diffracted beam must intersect the detector plane + + Parameters + ---------- + gvec_c : array_like + ([N,] 3) G-vector components in the CRYSTAL FRAME. + rmat_d : array_like + The (3, 3) COB matrix taking components in the + DETECTOR FRAME to the LAB FRAME + rmat_s : array_like + The ([N,] 3, 3) COB matrix taking components in the SAMPLE FRAME to the + LAB FRAME. It may be a single (3, 3) rotation matrix to use for all + gvec_c, or just one rotation matrix per gvec. + rmat_c : array_like + The (3, 3) COB matrix taking components in the + CRYSTAL FRAME to the SAMPLE FRAME + tvec_d : array_like + The (3, ) translation vector connecting LAB FRAME to DETECTOR FRAME + tvec_s : array_like + The (3, ) translation vector connecting LAB FRAME to SAMPLE FRAME + tvec_c : array_like + The ([M,] 3, ) translation vector(s) connecting SAMPLE FRAME to + CRYSTAL FRAME + beam_vec : array_like, optional + The (3, ) incident beam propagation vector components in the LAB FRAME; + the default is [0, 0, -1], which is the standard setting. + vmat_inv : array_like, optional + The (3, 3) matrix of inverse stretch tensor components in the + SAMPLE FRAME. The default is None, which implies a strain-free state + (i.e. V = I). + bmat : array_like, optional + The (3, 3) COB matrix taking components in the + RECIPROCAL LATTICE FRAME to the CRYSTAL FRAME; if supplied, it is + assumed that the input `gvecs` are G-vector components in the + RECIPROCL LATTICE FRAME (the default is None, which implies components + in the CRYSTAL FRAME) + + Returns + ------- + array_like + The ([M, ][N, ] 2) array of [x, y] diffracted beam intersections for + each of the N input G-vectors in the DETECTOR FRAME (all Z_d + coordinates are 0 and excluded) and for each of the M candidate + positions. For each input G-vector that cannot satisfy a Bragg + condition or intersect the detector plane, [NaN, Nan] is returned. + + Notes + ----- + Previously only a single candidate position was allowed. This is in + fact a vectored version of the previous API function. It is backwards + compatible, as passing single tvec_c is supported and has the same + result. + + """ + beam_vec = beam_vec if beam_vec is not None else cnst.beam_vec + + orig_ndim = gvec_c.ndim + gvec_c = np.ascontiguousarray(np.atleast_2d(gvec_c)) + rmat_d = np.ascontiguousarray(rmat_d) + rmat_s = np.ascontiguousarray(rmat_s) + rmat_c = np.ascontiguousarray(rmat_c) + tvec_d = np.ascontiguousarray(tvec_d.flatten()) + tvec_s = np.ascontiguousarray(tvec_s.flatten()) + tvec_c = np.ascontiguousarray(tvec_c.flatten()) + beam_vec = np.ascontiguousarray(beam_vec.flatten()) + + # depending on the number of dimensions of rmat_s use either the array + # version or the "scalar" (over rmat_s) version. Note that rmat_s is either + # a 3x3 matrix (ndim 2) or an nx3x4 array of matrices (ndim 3) + if rmat_s.ndim > 2: + result = transforms_c_api.gvecToDetectorXYArray( + gvec_c, rmat_d, rmat_s, rmat_c, tvec_d, tvec_s, tvec_c, beam_vec + ) + else: + result = transforms_c_api.gvecToDetectorXY( + gvec_c, rmat_d, rmat_s, rmat_c, tvec_d, tvec_s, tvec_c, beam_vec + ) + return result[0] if orig_ndim == 1 else result + + +def xy_to_gvec( + xy_d: NDArray[np.float64], + rmat_d: NDArray[np.float64], + rmat_s: NDArray[np.float64], + tvec_d: NDArray[np.float64], + tvec_s: NDArray[np.float64], + tvec_c: NDArray[np.float64], + rmat_b: Optional[NDArray[np.float64]] = None, + distortion: Optional[DistortionABC] = None, + output_ref: Optional[bool] = False, +) -> tuple[NDArray[np.float64], ...]: + """ + Takes a list cartesian (x, y) pairs in the DETECTOR FRAME and calculates + the associated reciprocal lattice (G) vectors and (bragg angle, azimuth) + pairs with respect to the specified beam and azimth (eta) reference + directions. + + Parameters + ---------- + xy_d : array_like + (n, 2) array of n (x, y) coordinates in DETECTOR FRAME + rmat_d : array_like + (3, 3) COB matrix taking components in the + DETECTOR FRAME to the LAB FRAME + rmat_s : array_like + (3, 3) COB matrix taking components in the + SAMPLE FRAME to the LAB FRAME + tvec_d : array_like + (3, ) translation vector connecting LAB FRAME to DETECTOR FRAME + tvec_s : array_like + (3, ) translation vector connecting LAB FRAME to SAMPLE FRAME + tvec_c : array_like + (3, ) translation vector connecting SAMPLE FRAME to CRYSTAL FRAME + rmat_b : array_like, optional + (3, 3) COB matrix taking components in the BEAM FRAME to the LAB FRAME; + defaults to None, which implies the standard setting of identity. + distortion : distortion class, optional + Default is None + output_ref : bool, optional + If True, prepends the apparent bragg angle and azimuth with respect to + the SAMPLE FRAME (ignoring effect of non-zero tvec_c) + + Returns + ------- + array_like + (n, 2) ndarray containing the (tth, eta) pairs associated with each + (x, y) associated with gVecs + array_like + (n, 3) ndarray containing the associated G vector directions in the + LAB FRAME + array_like, optional + if output_ref is True + """ + # It seems that the output_ref version is not present as the argument + # gets ignored + + rmat_b = rmat_b if rmat_b is not None else cnst.identity_3x3 + + # the code seems to ignore this argument, assume output_ref == True not + # implemented + assert not output_ref, 'output_ref not implemented' + + if distortion is not None: + xy_d = distortion.apply_inverse(xy_d) + + xy_d = np.ascontiguousarray(np.atleast_2d(xy_d)) + rmat_d = np.ascontiguousarray(rmat_d) + rmat_s = np.ascontiguousarray(rmat_s) + tvec_d = np.ascontiguousarray(tvec_d.flatten()) + tvec_s = np.ascontiguousarray(tvec_s.flatten()) + tvec_c = np.ascontiguousarray(tvec_c.flatten()) + beam_vec = np.ascontiguousarray((-rmat_b[:, 2]).flatten()) + eta_vec = np.ascontiguousarray(rmat_b[:, 0].flatten()) + return transforms_c_api.detectorXYToGvec( + xy_d, rmat_d, rmat_s, tvec_d, tvec_s, tvec_c, beam_vec, eta_vec + ) + + +def oscill_angles_of_hkls( + hkls: NDArray[np.float64], + chi: float, + rmat_c: NDArray[np.float64], + bmat: NDArray[np.float64], + wavelength: float, + v_inv: Optional[NDArray[np.float64]] = None, + beam_vec: NDArray[np.float64] = cnst.beam_vec, + eta_vec: NDArray[np.float64] = cnst.eta_vec, +) -> Tuple[NDArray[np.float64], NDArray[np.float64]]: + """ + Takes a list of unit reciprocal lattice vectors in crystal frame to the + specified detector-relative frame, subject to the conditions: + + 1) the reciprocal lattice vector must be able to satisfy a bragg condition + 2) the associated diffracted beam must intersect the detector plane + + Parameters + ---------- + hkls -- (n, 3) ndarray of n reciprocal lattice vectors + in the CRYSTAL FRAME + chi -- float representing the inclination angle of the + oscillation axis (std coords) + rmat_c -- (3, 3) ndarray, the COB taking CRYSTAL FRAME components + to SAMPLE FRAME + bmat -- (3, 3) ndarray, the COB taking RECIPROCAL LATTICE components + to CRYSTAL FRAME + wavelength -- float representing the x-ray wavelength in Angstroms + + Optional Keyword Arguments: + v_inv -- (6, ) ndarray, vec representation inverse stretch tensor + components in the SAMPLE FRAME. The default is None, which + implies a strain-free state (i.e. V = I). + beam_vec -- (3, 1) mdarray containing the incident beam direction + components in the LAB FRAME + eta_vec -- (3, 1) mdarray containing the reference azimuth direction + components in the LAB FRAME + + Outputs: + ome0 -- (n, 3) ndarray containing the feasible (tTh, eta, ome) triplets + for each input hkl (first solution) + ome1 -- (n, 3) ndarray containing the feasible (tTh, eta, ome) triplets + for each input hkl (second solution) + + Notes: + ------------------------------------------------------------------------ + The reciprocal lattice vector, G, will satisfy the the Bragg condition + when: + + b.T * G / ||G|| = -sin(theta) + + where b is the incident beam direction (k_i) and theta is the Bragg + angle consistent with G and the specified wavelength. The components of + G in the lab frame in this case are obtained using the crystal + orientation, Rc, and the single-parameter oscillation matrix, Rs(ome): + + Rs(ome) * Rc * G / ||G|| + + The equation above can be rearranged to yield an expression of the form: + + a*sin(ome) + b*cos(ome) = c + + which is solved using the relation: + + a*sin(x) + b*cos(x) = sqrt(a**2 + b**2) * sin(x + alpha) + + --> sin(x + alpha) = c / sqrt(a**2 + b**2) + + where: + + alpha = atan2(b, a) + + The solutions are: + + / + | arcsin(c / sqrt(a**2 + b**2)) - alpha + x = < + | pi - arcsin(c / sqrt(a**2 + b**2)) - alpha + \ + + There is a double root in the case the reflection is tangent to the + Debye-Scherrer cone (c**2 = a**2 + b**2), and no solution if the + Laue condition cannot be satisfied (filled with NaNs in the results + array here) + """ + # this was oscillAnglesOfHKLs + hkls = np.array(hkls, dtype=float, order='C') + if v_inv is None: + v_inv = np.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]) + else: + v_inv = np.ascontiguousarray(v_inv.flatten()) + rmat_c = np.ascontiguousarray(rmat_c) + beam_vec = np.ascontiguousarray(beam_vec.flatten()) + eta_vec = np.ascontiguousarray(eta_vec.flatten()) + bmat = np.ascontiguousarray(bmat) + return transforms_c_api.oscillAnglesOfHKLs( + hkls, chi, rmat_c, bmat, wavelength, v_inv, beam_vec, eta_vec + ) + + +def unit_vector(vec_in: NDArray[np.float64]) -> NDArray[np.float64]: + """ + Normalize the input vector(s) to unit length. + + Parameters + ---------- + vec_in : ndarray + The input vector(s) (n,3) to normalize. + + Returns + ------- + ndarray + The normalized vector(s) of the same shape as the input. + + Raises + ------ + ValueError + If the input vector(s) do not have the correct shape. + """ + vec_in = np.ascontiguousarray(vec_in) + if vec_in.ndim == 1: + return transforms_c_api.unitRowVector(vec_in) + elif vec_in.ndim == 2: + return transforms_c_api.unitRowVectors(vec_in) + else: + raise ValueError( + "incorrect arg shape; must be 1-d or 2-d, yours is %d-d" % (vec_in.ndim) + ) + + +def make_detector_rmat(tilt_angles: NDArray[np.float64]) -> NDArray[np.float64]: + """ + Form the (3, 3) tilt rotations from the tilt angle list: + + tilt_angles = [gamma_Xl, gamma_Yl, gamma_Zl] in radians + Returns the (3, 3) rotation matrix from the detector frame to the lab frame + """ + t_angles = np.ascontiguousarray(np.r_[tilt_angles].flatten()) + return transforms_c_api.makeDetectorRotMat(t_angles) + + +# make_sample_rmat in CAPI is split between makeOscillRotMat +# and makeOscillRotMatArray... + + +def make_oscill_rmat(oscillAngles): + chi, ome = oscillAngles + ome = np.atleast_1d(ome) + result = transforms_c_api.makeOscillRotMat(chi, ome) + return result.reshape((3, 3)) + + +def make_oscill_rmat_array(chi, omeArray): + arg = np.ascontiguousarray(omeArray) + return transforms_c_api.makeOscillRotMat(chi, arg) + + +def make_sample_rmat( + chi: float, ome: float | NDArray[np.float64] +) -> NDArray[np.float64]: + # TODO: Check this docstring + """ + Make a rotation matrix representing the COB from sample frame to the lab + frame. + + Parameters + ---------- + chi : float + The inclination angle of the sample frame about the lab frame Y. + ome : float or ndarray + The oscillation angle(s) about the sample frame Y. + + Returns + ------- + ndarray + A 3x3 rotation matrix representing the input oscillation angles. + + """ + ome_array = np.atleast_1d(ome) + if ome is ome_array: + ome_array = np.ascontiguousarray(ome_array) + result = cpp_transforms.makeOscillRotMat(chi, ome_array).reshape( + len(ome_array), 3, 3 + ) + else: + # converted to 1d array of 1 element, no need + # to call ascontiguousarray, but need to remove + # the outer dimension from the result + result = cpp_transforms.makeOscillRotMat(chi, ome_array) + + return result + + +def make_rmat_of_expmap( + exp_map: Sequence[float] | NDArray[np.float64], +) -> NDArray[np.float64]: + """ + Calculate the rotation matrix of an exponential map. + + Parameters + ---------- + exp_map : array_like + A 3-element sequence representing the exponential map n*phi. + + Returns + ------- + NDArray[np.float64] + A 3x3 rotation matrix representing the input exponential map + """ + arg = np.ascontiguousarray(exp_map).flatten() + return cpp_transforms.make_rot_mat_of_exp_map(arg) + + +def make_binary_rmat(axis: NDArray[np.float64]) -> NDArray[np.float64]: + """ + Create a 180 degree rotation matrix about the input axis. + + Parameters + ---------- + axis : ndarray + 3-element sequence representing the rotation axis, in radians. + + Returns + ------- + ndarray + A 3x3 rotation matrix representing a 180 degree rotation about the + input axis. + """ + + arg = np.ascontiguousarray(axis.flatten()) + return cpp_transforms.make_binary_rot_mat(arg) + + +def make_beam_rmat( + bvec_l: NDArray[np.float64], evec_l: NDArray[np.float64] +) -> NDArray[np.float64]: + """ + Creates a COB matrix from the beam frame to the lab frame + Note: beam and eta vectors must not be colinear + + Parameters + ---------- + bvec_l : ndarray + The beam vector in the lab frame, The (3, ) incident beam propagation + vector components in the lab frame. + the default is [0, 0, -1], which is the standard setting. + evec_l : ndarray + Vector defining eta=0 in the lab frame. Defaults to [1,0,0] + """ + arg1 = np.ascontiguousarray(bvec_l.flatten()) + arg2 = np.ascontiguousarray(evec_l.flatten()) + return transforms_c_api.makeEtaFrameRotMat(arg1, arg2) + + +def validate_angle_ranges( + ang_list: float | NDArray[np.float64], + start_angs: float | NDArray[np.float64], + stop_angs: float | NDArray[np.float64], + ccw: bool = True, +) -> NDArray[np.bool_]: + """ + Find out if angles are in the CCW or CW range from start to stop + + Parameters + ---------- + ang_list : ndarray + The list of angles to validate + start_angs : ndarray + The starting angles + stop_angs : ndarray + The stopping angles + ccw : bool, optional + If True, the angles are in the CCW range from start to stop. + Defaults to True. + + Returns + ------- + ndarray + List of bools indicating if the angles are in the correct range + + """ + ang_list = np.atleast_1d(ang_list).astype(np.double, order='C') + start_angs = np.atleast_1d(start_angs).astype(np.double, order='C') + stop_angs = np.atleast_1d(stop_angs).astype(np.double, order='C') + + return transforms_c_api.validateAngleRanges(ang_list, start_angs, stop_angs, ccw) + + +def rotate_vecs_about_axis( + angle: float | NDArray[np.float64], + axis: NDArray[np.float64], + vecs: NDArray[np.float64], +) -> NDArray[np.float64]: + """ + Rotate vectors about an axis + + Parameters + ---------- + angle : array_like + Array of angles (len==1 or n) + axis : ndarray + Array of unit vectors (shape == (3,1) or (3,n)) + vecs : ndarray + Array of vectors to rotate (shape == (3,n)) + + Returns + ------- + ndarray + Array of rotated vectors (shape == (3,n)) + """ + angle = np.asarray(angle) + axis = np.ascontiguousarray(axis.T) + vecs = np.ascontiguousarray(vecs.T) + result = transforms_c_api.rotate_vecs_about_axis(angle, axis, vecs) + return result.T + + +def quat_distance( + q1: NDArray[np.float64], q2: NDArray[np.float64], qsym: NDArray[np.float64] +) -> NDArray[np.float64]: + """Distance between two quaternions, taking quaternions of symmetry into account. + + Parameters + ---------- + q1 : arary_like + First quaternion. + q2 : arary_like + Second quaternion. + qsym : ndarray + List of symmetry quaternions. + + Returns + ------- + float + The distance between the two quaternions. + """ + q1 = np.ascontiguousarray(q1.flatten()) + q2 = np.ascontiguousarray(q2.flatten()) + # C module expects quaternions in row major, numpy code in column major. + qsym = np.ascontiguousarray(qsym.T) + return transforms_c_api.quat_distance(q1, q2, qsym) \ No newline at end of file diff --git a/hexrd/powder/wppf/texture.py b/hexrd/powder/wppf/texture.py index 9e3d2d0bc..032dca7e0 100644 --- a/hexrd/powder/wppf/texture.py +++ b/hexrd/powder/wppf/texture.py @@ -5,7 +5,7 @@ # hexrd imports # ------------- -from hexrd.core.transforms.xfcapi import anglesToGVec +from hexrd.core.transforms.xfcapi import angles_to_gvec from hexrd.core import constants from hexrd.powder.wppf.phase import Material_Rietveld @@ -907,7 +907,7 @@ def precompute_spherical_harmonics( ) ).T - pfgrid = anglesToGVec(angs, bHat_l=self.bvec, eHat_l=self.etavec) + pfgrid = angles_to_gvec(angs, beam_vec=self.bvec, eta_vec=self.etavec) # Next part runs in a numba function to accelerate it. if calc_type == 'texture_factor': diff --git a/setup.py b/setup.py index 71d5d9200..5bef58817 100644 --- a/setup.py +++ b/setup.py @@ -154,23 +154,9 @@ def get_cpp_extensions(): return [transforms_ext, inverse_distortion_ext] -def get_old_xfcapi_extension_modules(): - # for transforms - srclist = ['transforms_CAPI.c', 'transforms_CFUNC.c'] - srclist = [os.path.join('hexrd/core/transforms', f) for f in srclist] - transforms_mod = Extension( - 'hexrd.core.extensions._transforms_CAPI', - sources=srclist, - include_dirs=[np_include_dir], - extra_compile_args=compiler_flags, - ) - - return [transforms_mod] - - def get_new_xfcapi_extension_modules(): transforms_mod = Extension( - 'hexrd.core.extensions._new_transforms_capi', + 'hexrd.core.extensions.transforms_c_api', sources=['hexrd/core/transforms/new_capi/module.c'], include_dirs=[np_include_dir], extra_compile_args=compiler_flags, @@ -184,7 +170,6 @@ def get_extension_modules(): return [ item for sublist in ( - get_old_xfcapi_extension_modules(), get_new_xfcapi_extension_modules(), get_convolution_extensions(), get_cpp_extensions(), diff --git a/tests/transforms/common.py b/tests/transforms/common.py index 3afa7b7a0..997815355 100644 --- a/tests/transforms/common.py +++ b/tests/transforms/common.py @@ -3,7 +3,7 @@ import numpy as np import hexrd.core.constants as ct -from hexrd.core.transforms.new_capi.xf_new_capi import unit_vector +from hexrd.core.transforms.xfcapi import unit_vector def convert_axis_angle_to_rmat(axis, angle): diff --git a/tests/transforms/test_angles_to_dvec_from_file.py b/tests/transforms/test_angles_to_dvec_from_file.py index c4f8a7c0e..5ab67a585 100644 --- a/tests/transforms/test_angles_to_dvec_from_file.py +++ b/tests/transforms/test_angles_to_dvec_from_file.py @@ -4,7 +4,7 @@ from __future__ import absolute_import import numpy as np -from hexrd.core.transforms.new_capi.xf_new_capi import angles_to_dvec +from hexrd.core.transforms.xfcapi import angles_to_dvec # from common import random_rotation_matrix, random_unit_vectors diff --git a/tests/transforms/test_angles_to_gvec_from_file.py b/tests/transforms/test_angles_to_gvec_from_file.py index 0273d7a40..7ed9c66dc 100644 --- a/tests/transforms/test_angles_to_gvec_from_file.py +++ b/tests/transforms/test_angles_to_gvec_from_file.py @@ -5,7 +5,7 @@ from __future__ import absolute_import import numpy as np -from hexrd.core.transforms.new_capi.xf_new_capi import angles_to_gvec +from hexrd.core.transforms.xfcapi import angles_to_gvec # from common import random_rotation_matrix, random_unit_vectors diff --git a/tests/transforms/test_gvec_to_xy.py b/tests/transforms/test_gvec_to_xy.py index a681ca01d..a12db71e1 100644 --- a/tests/transforms/test_gvec_to_xy.py +++ b/tests/transforms/test_gvec_to_xy.py @@ -12,7 +12,7 @@ from common import convert_axis_angle_to_rmat -from hexrd.core.transforms.new_capi.xf_new_capi import gvec_to_xy +from hexrd.core.transforms.xfcapi import gvec_to_xy # gvec_to_xy intersects vectors from crystal position with the detector plane. diff --git a/tests/transforms/test_gvec_to_xy_from_file.py b/tests/transforms/test_gvec_to_xy_from_file.py index 9f9bd5ebc..6563d0abb 100644 --- a/tests/transforms/test_gvec_to_xy_from_file.py +++ b/tests/transforms/test_gvec_to_xy_from_file.py @@ -4,7 +4,7 @@ from __future__ import absolute_import import numpy as np -from hexrd.core.transforms.new_capi.xf_new_capi import gvec_to_xy +from hexrd.core.transforms.xfcapi import gvec_to_xy # from common import random_rotation_matrix, random_unit_vectors diff --git a/tests/transforms/test_make_beam_rmat_from_file.py b/tests/transforms/test_make_beam_rmat_from_file.py index c3c8ea7bd..851ccb6fc 100644 --- a/tests/transforms/test_make_beam_rmat_from_file.py +++ b/tests/transforms/test_make_beam_rmat_from_file.py @@ -5,7 +5,7 @@ from __future__ import absolute_import import numpy as np -from hexrd.core.transforms.new_capi.xf_new_capi import make_beam_rmat +from hexrd.core.transforms.xfcapi import make_beam_rmat # from common import random_unit_vectors diff --git a/tests/transforms/test_make_binary_rmat.py b/tests/transforms/test_make_binary_rmat.py index ff809cb12..2666b1145 100644 --- a/tests/transforms/test_make_binary_rmat.py +++ b/tests/transforms/test_make_binary_rmat.py @@ -5,7 +5,7 @@ from __future__ import absolute_import import numpy as np -from hexrd.transforms.new_capi.xf_new_capi import make_binary_rmat +from hexrd.transforms.xfcapi import make_binary_rmat from hexrd.rotations import quatOfAngleAxis, rotMatOfQuat from transforms.common import random_unit_vectors diff --git a/tests/transforms/test_make_detector_rmat_from_file.py b/tests/transforms/test_make_detector_rmat_from_file.py index 7d8343ee1..1cea16592 100644 --- a/tests/transforms/test_make_detector_rmat_from_file.py +++ b/tests/transforms/test_make_detector_rmat_from_file.py @@ -5,7 +5,7 @@ from __future__ import absolute_import import numpy as np -from hexrd.core.transforms.new_capi.xf_new_capi import make_detector_rmat +from hexrd.core.transforms.xfcapi import make_detector_rmat def test_make_detector_rmat_from_file(test_data_dir): diff --git a/tests/transforms/test_make_rmat_of_expmap_from_file.py b/tests/transforms/test_make_rmat_of_expmap_from_file.py index 678194a68..dcb9d7723 100644 --- a/tests/transforms/test_make_rmat_of_expmap_from_file.py +++ b/tests/transforms/test_make_rmat_of_expmap_from_file.py @@ -5,7 +5,7 @@ from __future__ import absolute_import import numpy as np -from hexrd.core.transforms.new_capi.xf_new_capi import make_rmat_of_expmap +from hexrd.core.transforms.xfcapi import make_rmat_of_expmap def test_make_rmat_of_expmap_from_file(test_data_dir): diff --git a/tests/transforms/test_make_sample_rmat_from_file.py b/tests/transforms/test_make_sample_rmat_from_file.py index 49657994e..2a1fecf2b 100644 --- a/tests/transforms/test_make_sample_rmat_from_file.py +++ b/tests/transforms/test_make_sample_rmat_from_file.py @@ -4,7 +4,7 @@ from __future__ import absolute_import import numpy as np -from hexrd.core.transforms.new_capi.xf_new_capi import make_sample_rmat +from hexrd.core.transforms.xfcapi import make_sample_rmat def test_make_sample_rmat_from_file(test_data_dir): diff --git a/tests/transforms/test_quat_distance_from_file.py b/tests/transforms/test_quat_distance_from_file.py index b7d35e655..55be3f0b5 100644 --- a/tests/transforms/test_quat_distance_from_file.py +++ b/tests/transforms/test_quat_distance_from_file.py @@ -4,7 +4,7 @@ from __future__ import absolute_import import numpy as np -from hexrd.core.transforms.new_capi.xf_new_capi import quat_distance +from hexrd.core.transforms.xfcapi import quat_distance # from common import random_unit_vectors # from hexrd.core.rotations import quatOfLaueGroup diff --git a/tests/transforms/test_rotate_vecs_about_axis.py b/tests/transforms/test_rotate_vecs_about_axis.py index 94200285f..d6f0df1e4 100644 --- a/tests/transforms/test_rotate_vecs_about_axis.py +++ b/tests/transforms/test_rotate_vecs_about_axis.py @@ -1,4 +1,4 @@ -from hexrd.core.transforms.new_capi.xf_new_capi import rotate_vecs_about_axis +from hexrd.core.transforms.xfcapi import rotate_vecs_about_axis import numpy as np diff --git a/tests/transforms/test_validate_angle_ranges_from_file.py b/tests/transforms/test_validate_angle_ranges_from_file.py index 6c43012c4..f2d3fb0b9 100644 --- a/tests/transforms/test_validate_angle_ranges_from_file.py +++ b/tests/transforms/test_validate_angle_ranges_from_file.py @@ -4,7 +4,7 @@ from __future__ import absolute_import import numpy as np -from hexrd.core.transforms.new_capi.xf_new_capi import validate_angle_ranges +from hexrd.core.transforms.xfcapi import validate_angle_ranges def test_validate_angle_ranges_from_file(test_data_dir): diff --git a/tests/transforms/test_xy_to_gvec.py b/tests/transforms/test_xy_to_gvec.py index 54cddf1cf..899a71afd 100644 --- a/tests/transforms/test_xy_to_gvec.py +++ b/tests/transforms/test_xy_to_gvec.py @@ -9,7 +9,7 @@ from collections import namedtuple import pytest import numpy as np -from hexrd.core.transforms.new_capi.xf_new_capi import xy_to_gvec +from hexrd.core.transforms.xfcapi import xy_to_gvec Experiment = namedtuple( From 2d1a2858f06277ee1ada7350dd2f4d9691ef419f Mon Sep 17 00:00:00 2001 From: Zack Singer Date: Tue, 17 Feb 2026 14:51:55 -0500 Subject: [PATCH 2/2] Add deprecation exceptions to old function calls. --- hexrd/core/transforms/xfcapi.py | 113 +++++++++++++++++++++++++++++++- 1 file changed, 112 insertions(+), 1 deletion(-) diff --git a/hexrd/core/transforms/xfcapi.py b/hexrd/core/transforms/xfcapi.py index fb7e1cc7f..0cc167c86 100644 --- a/hexrd/core/transforms/xfcapi.py +++ b/hexrd/core/transforms/xfcapi.py @@ -688,4 +688,115 @@ def quat_distance( q2 = np.ascontiguousarray(q2.flatten()) # C module expects quaternions in row major, numpy code in column major. qsym = np.ascontiguousarray(qsym.T) - return transforms_c_api.quat_distance(q1, q2, qsym) \ No newline at end of file + return transforms_c_api.quat_distance(q1, q2, qsym) + +#################### +# Sunset functions # +#################### + + +def anglesToGVec(*args, **kwargs): + """Deprecated alias for angles_to_gvec.""" + raise RuntimeError( + "anglesToGVec is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#anglestogvec" + ) + + +def anglesToDVec(*args, **kwargs): + """Deprecated alias for angles_to_dvec.""" + raise RuntimeError( + "anglesToDVec is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#anglestodvec" + ) + + +def detectorXYToGvec(*args, **kwargs): + """Deprecated alias for xy_to_gvec.""" + raise RuntimeError( + "detectorXYToGvec is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#detectorxytogvec" + ) + + +def gvecToDetectorXY(*args, **kwargs): + """Deprecated alias for gvec_to_xy.""" + raise RuntimeError( + "gvecToDetectorXY is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#gvectodetectorxy" + ) + + +def gvecToDetectorXYArray(*args, **kwargs): + """Deprecated alias for gvec_to_xy.""" + raise RuntimeError( + "gvecToDetectorXYArray is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#gvectodetectorxyarray" + ) + + +def oscillAnglesOfHKLs(*args, **kwargs): + """Deprecated alias for oscill_angles_of_hkls.""" + raise RuntimeError( + "oscillAnglesOfHKLs is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#oscillanglesofhkls" + ) + + +def angularDifference(*args, **kwargs): + """Deprecated alias for angular_difference.""" + raise RuntimeError( + "angularDifference is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#angulardifference" + ) + + +def makeEtaFrameRotMat(*args, **kwargs): + """Deprecated alias for make_beam_rmat.""" + raise RuntimeError( + "makeEtaFrameRotMat is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#makeetaframerotmat" + ) + + +def makeOscillRotMat(*args, **kwargs): + """Deprecated alias for make_oscill_rmat.""" + raise RuntimeError( + "makeOscillRotMat is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#makeoscillrotmat" + ) + + +def makeOscillRotMatArray(*args, **kwargs): + """Deprecated alias for make_oscill_rmat_array.""" + raise RuntimeError( + "makeOscillRotMatArray is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#makeoscillrotmatarray" + ) + + +def makeRotMatOfExpMap(*args, **kwargs): + """Deprecated alias for make_rmat_of_expmap.""" + raise RuntimeError( + "makeRotMatOfExpMap is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#makerotmatofexpmap" + ) + + +def makeRotMatOfQuat(*args, **kwargs): + """Deprecated alias for make_rot_mat_of_quat.""" + raise RuntimeError( + "makeRotMatOfQuat is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#makerotmatofquat" + ) + + +def mapAngle(*args, **kwargs): + """Deprecated alias for map_angle.""" + raise RuntimeError( + "mapAngle is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#mapangle" + ) + + +def rowNorm(*args, **kwargs): + """Deprecated alias for row_norm.""" + raise RuntimeError( + "rowNorm is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#rownorm" + ) + + +def unitRowVector(*args, **kwargs): + """Deprecated alias for unit_vector.""" + raise RuntimeError( + "unitRowVector is deprecated, refer to the transition instructions here: https://github.com/HEXRD/hexrd/wiki/old_xfcapi-Conversion-Guide#unitrowvector" + ) + +