'use strict';
/**
* @fileoverview TensorData class to store tensors like NMR data and such.
* @module
*/
import _ from 'lodash';
import * as mjs from 'mathjs';
import { arraysAlmostEqual, rotate_matrix } from './utils.js';
import { RelativeTensorOrientation } from './relative-orientation.js';
const efg2hz = 234964.77815245767;
const isc2hz = 1.6784031762379067e-16; // hbar/(2*pi)*1e19
const PI = mjs.pi
function createRz(angle) {
return mjs.matrix([
[mjs.cos(angle), -mjs.sin(angle), 0],
[mjs.sin(angle), mjs.cos(angle), 0],
[0, 0, 1]
]);
}
function createRy(angle) {
return mjs.matrix([
[mjs.cos(angle), 0, mjs.sin(angle)],
[0, 1, 0],
[-mjs.sin(angle), 0, mjs.cos(angle)]
]);
}
function createRx(angle) {
return mjs.matrix([
[1, 0, 0],
[0, mjs.cos(angle), -mjs.sin(angle)],
[0, mjs.sin(angle), mjs.cos(angle)]
]);
}
/**
*
* @param {*} alpha in radians
* @param {*} beta in radians
* @param {*} gamma in radians
* @param {*} PAS (principal axis system)
* @param {*} mode (e.g. 'zyz')
* @param {*} active (true: active rotation, false: passive rotation)
*
* @returns {mathjs.Matrix} rotation matrix
*
*/
function eulerRotation(alpha, beta, gamma, mode = 'zyz', active = true) {
if (mode != 'zyz' && mode != 'zxz')
throw new Error('Only zyz and zxz modes are implemented so far');
// Initialize the rotation matrix as diagonal
let R = mjs.matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]
]);
let rz1, rx2, rz3;
if (mode === 'zyz') {
rz1 = createRz(alpha);
rx2 = createRy(beta);
rz3 = createRz(gamma);
}
else if (mode === 'zxz') {
rz1 = createRz(alpha);
rx2 = createRx(beta);
rz3 = createRz(gamma);
}
if (active) {
R = mjs.multiply(mjs.multiply(rz1, rx2), rz3);
} else {
rz1 = mjs.inv(rz1);
rx2 = mjs.inv(rx2);
rz3 = mjs.inv(rz3);
R = mjs.multiply(mjs.multiply(rz3, rx2), rz1);
}
return R;
}
function rotateTensor(alpha, beta, gamma, PAS, mode, active) {
const R = eulerRotation(alpha, beta, gamma, mode, active);
// R*PAS R^-1
return mjs.multiply(mjs.multiply(R, PAS), mjs.inv(R));
}
class TensorData {
/**
* Create a TensorData object, to store whole tensors,
* diagonalise their symmetric part and return derived quantities
*
* @param {Array | mathjs.Matrix | TensorData} M Tensor in 3x3 matrix form
*
*/
constructor(M) {
M = convertToMatrix(M);
this._M = M.clone();
var MT = mjs.transpose(M);
this._Msymm = mjs.divide(mjs.add(M, MT), 2.0);
this._Masymm = mjs.subtract(M, this._Msymm);
// Diagonalize
var eigs = mjs.eigs(this._Msymm);
// Sort by eigenvalue magnitude
var evecs = eigs.eigenvectors.map(e => e.vector._data);
eigs = _.zip(eigs.values._data, evecs);
eigs = _.sortBy(eigs, function(x) {
return x[0];
});
eigs = _.unzip(eigs);
this._evals = eigs[0];
// store original eigenvectors and eigenvalues (before making them right-handed)
this._orig_evals = Array.from(this._evals);
this._orig_evecs = mjs.transpose(eigs[1]);
// Make it right-handed
evecs = eigs[1];
evecs[2] = mjs.cross(evecs[0], evecs[1]);
this._evecs = mjs.transpose(evecs);
// Isotropy
this._iso = mjs.mean(this._evals);
// number of degeneracies
const EPS = 1e-6;
const degeneracies = _.countBy(this._evals, function(x) {
return mjs.round(x/EPS)*EPS;
});
// symmetry
// 2 if all eigenvalues are the same
// 1 if two eigenvalues are the same
// 0 if all eigenvalues are different
this._symmetry = this._evals.length - _.keys(degeneracies).length;
// Save evals and evecs in Haeberlen and NQR order for convenience
this._haeb_evals = this.sorted_eigenvalues('haeberlen');
this._haeb_evecs = this.sorted_eigenvectors('haeberlen');
this._nqr_evals = this.sorted_eigenvalues('nqr');
this._nqr_evecs = this.sorted_eigenvectors('nqr');
// Default convention is increasing
this._convention = 'increasing';
}
/**
* This method sorts the eigenvalues and eigenvectors in Haeberlen and NQR order.
*
* Haeberlen Order:
* The method first calculates the absolute difference between each eigenvalue and the isotropic value.
* These differences are then sorted in ascending order. The sorted eigenvalues are stored in `this._haeb_evals`.
* The corresponding eigenvectors are also sorted in the same order, but with the first and second elements swapped.
* The sorted eigenvectors are stored in `this._haeb_evecs`.
*
* NQR Order:
* The method then repeats the same process for the original eigenvalues, but without swapping the first and second elements.
* The sorted original eigenvalues are stored in `this._nqr_evals`.
*
* Note: The method uses the lodash (`_`) library for the `zip`, `range`, and `sortBy` functions.
* The `mjs.transpose` function is used to transpose the array of eigenvectors, and `mjs.cross` is used to calculate the cross product of two vectors.
*/
sort_eigs() {
var evecs = mjs.transpose(this._evecs);
// Haeberlen order
var iso = this._iso;
var haeb = _.zip(_.range(3), this._evals);
haeb = _.sortBy(haeb, function(x) {
return mjs.abs(x[1] - iso);
});
// x and y are swapped in this order
this._haeb_evals = [
this._evals[haeb[1][0]],
this._evals[haeb[0][0]],
this._evals[haeb[2][0]]
];
this._haeb_evecs = mjs.transpose([
evecs[haeb[1][0]],
evecs[haeb[0][0]],
mjs.cross(evecs[haeb[1][0]], evecs[haeb[0][0]])
]);
// NQR order
var nqr = _.zip(_.range(3), this._orig_evals);
// same sorting function as for Haeberlen
nqr = _.sortBy(nqr, function(x) {
return mjs.abs(x[1] - iso);
});
// but in this case we don't swap the x and y
this._nqr_evals = [
this._orig_evals[nqr[0][0]],
this._orig_evals[nqr[1][0]],
this._orig_evals[nqr[2][0]]
];
let nqr_evecs = mjs.transpose(this._orig_evecs);
this._nqr_evecs = mjs.transpose([
nqr_evecs[nqr[0][0]],
nqr_evecs[nqr[1][0]],
mjs.cross(nqr_evecs[nqr[0][0]], nqr_evecs[nqr[1][0]])
]);
}
get data() {
return JSON.parse(JSON.stringify(this._M._data));
}
get symmetric() {
return JSON.parse(JSON.stringify(this._Msymm._data));
}
get asymmetric() {
return JSON.parse(JSON.stringify(this._Masymm._data));
}
get convention() {
return this._convention;
}
set convention(convention) {
this._convention = convention;
}
get eigenvalues() {
return Array.from(this.sorted_eigenvalues(this.convention));
}
get eigenvectors() {
return JSON.parse(JSON.stringify(this.sorted_eigenvectors(this.convention)));
}
get haeberlen_eigenvalues() {
return Array.from(this._haeb_evals);
}
get haeberlen_eigenvectors() {
return JSON.parse(JSON.stringify(this._haeb_evecs));
}
get nqr_eigenvalues() {
return Array.from(this._nqr_evals);
}
get nqr_eigenvectors() {
return JSON.parse(JSON.stringify(this._nqr_evecs));
}
get isotropy() {
return this._iso;
}
get anisotropy() {
return this._haeb_evals[2] - (this._haeb_evals[0] + this._haeb_evals[1]) / 2.0;
}
get reduced_anisotropy() {
return this._haeb_evals[2] - this._iso;
}
get asymmetry() {
var ra = this.reduced_anisotropy;
return (this._haeb_evals[1] - this._haeb_evals[0]) / ra;
}
get span() {
return this._evals[2] - this._evals[0];
}
get skew() {
var s = this.span;
return 3 * (this._evals[1] - this._iso) / s;
}
get symmetry() {
// 2 if all eigenvalues are the same
// 1 if two eigenvalues are the same
// 0 if all eigenvalues are different
return this._symmetry;
}
/**
* Rotate the TensorData by a given basis, either as passive or active
* transformation. Returns the rotated TensorData (does not modify this in
* place). Default is passive. The convention is such that for a symmetric
* tensor,
*
* T.rotate(T.eigenvectors)
*
* returns the diagonalised tensor.
*
* @param {Array | mathjs.Matrix | TensorData} basis Basis to rotate into
* @param {Boolean} active If true, make it an active transformation (default is false)
*
* @return {TensorData} Rotated tensor
*/
rotate(basis, active = false) {
// Rotate the tensor by the given basis of vectors
if (basis instanceof mjs.Matrix)
basis = basis._data;
if (basis instanceof TensorData)
basis = basis._M._data;
var bR = basis;
var bL = mjs.transpose(basis);
if (active) {
bR = bL;
bL = basis;
}
var rdata = mjs.multiply(bL, this._M._data, bR);
return new TensorData(rdata);
}
/**
* Rotate the TensorData by a given angle around a given axis.
* Uses the rotate_matrix utility function.
*
* @param {number|Array} a - The angle that the matrix is rotated around the vector 'v'. 'a'
* can also be a vector and then 'a' is rotated into 'v'. Angles are in radians.
* @param {string|Array} v - Vector to rotate the matrix around. Vectors can be given as
* strings: 'x', '-x', 'y', ... .These correspond to the *local* coordinate system of the tensor if intrinsic == true.
* @param {boolean} intrinsic - If true, rotate the tensor intrinsically (default is true). i.e. the tensor is rotated in the local coordinate system.
*
* @return {TensorData} Rotated tensor
*
*/
rotateByAngleAxis(angle, axis, intrinsic = true) {
// if (intrinsic) {
// // convert axis to local coordinates
// axis = mjs.multiply(mjs.inv(this.eigenvectors), string2vector(axis));
// if (typeof angle != 'number') {
// angle = mjs.multiply(mjs.inv(this.eigenvectors), string2vector(angle));
// }
// }
if (typeof angle === 'number') {
angle = angle * 180 / PI; // convert to degrees for the rotate_matrix function
}
// Rotate the tensor by the given angle around the given axis
let M = rotate_matrix(this._M._data, angle, axis);
return new TensorData(M);
}
/**
* Calculate the Euler angles for the TensorData. Returns an array of
* [alpha, beta, gamma], in radians unless degrees is true.
*
* This function is based on the code from TensorView for MATLAB by
* Leo Svenningsson and Leonard J. Mueller
* https://doi.org/10.1016/j.ssnmr.2022.101849
*
* @param {String} mode Mode to use: either 'zyz' or 'zxz' (default is 'zyz')
* @param {Boolean} active If true, make it an active transformation (default is true)
* @param {String} convention Eval ordering convention to use (default is whatever is set in 'this.convention)
* @param {Boolean} degrees If true, return angles in degrees (default is false)
* @param {Number} EPS Tolerance for comparing eigenvalues etc. (default is 1e-6)
*
* @return {Array} Array of [alpha, beta, gamma] in radians
*
* TODO: check edge cases, esp. axial symmetry
* TODO: add in conversion from radians to degrees
*/
euler(mode = 'zyz', active = true, eval_convention = null, degrees = false, EPS = 1e-6) {
let convention = eval_convention != null ? eval_convention : this.convention;
// Eigenvalues and eigenvectors of the symmetric part of the tensor
const sorted_evals = this.sorted_eigenvalues(convention);
// deep copy of the eigenvectors
let U = JSON.parse(JSON.stringify(this.sorted_eigenvectors(convention)));
// invert U if not active
// Remove this for now as I've moved it into the eulerFromU function
// if (!active) {
// // I'm not sure why we need to multiply by -1 here
// // but I was getting the wrong alpha angles without it
// // (wrong wrt TensorView/Soprano)
// U = mjs.inv(mjs.multiply(U, -1));
// }
let alpha = 0
let beta = 0
let gamma = 0
// number of unique eigenvalues (within tolerance)
const symm = this.symmetry;
if (symm == 2) {
// all eigenvalues are the same -> spherical symmetry
// no need to calculate the Euler angles
return [0, 0, 0];
}
// PAS = principal axis system
// it has the eigenvalues on the diagonals, zeros elsewhere
const PAS = new mjs.diag(sorted_evals);
// Euler angles
// initialising the angles
[alpha, beta, gamma] = eulerFromU(U, mode, active, EPS);
let M = rotateTensor(alpha, beta, gamma, PAS, mode, active).toArray();
// U PAS U^-1
let CF = mjs.multiply(mjs.multiply(U, PAS), mjs.inv(U));
// now let's check to make sure M and Cf are nearly the same
// this is because the eigenvalue solver is not always perfect
// and we want to make sure that the angles are correct
// TODO! this doesn't seem to be working since the rotation doesn't give us
// the same tensor as the diagonalisation even after the U -> -U trick
if (!arraysAlmostEqual(M, CF, EPS)) {
console.log('M and CF are not the same. Trying again with -U');
// if they are not the same, we can simply flip all the signs of U
// and try again
U = mjs.multiply(-1, U);
[alpha, beta, gamma] = eulerFromU(U, mode, active, EPS);
M = rotateTensor(alpha, beta, gamma, PAS, mode, active);
CF = mjs.multiply(mjs.multiply(U, PAS), mjs.inv(U));
// if they are still not the same, then we have a problem
if (!arraysAlmostEqual(M, CF, EPS)) {
console.warn('M and CF are still not the same. Something is wrong.');
// throw new Error('M and CF are still not the same. Something is wrong.');
}
}
// standardise/normalise angle ranges is done inside the handleEulerEdgeCases function
[alpha, beta, gamma] = handleEulerEdgeCases([alpha, beta, gamma], sorted_evals, this.symmetric, mode, active, EPS);
// convert to degrees if necessary
if (degrees) {
alpha = mjs.multiply(alpha, 180 / PI);
beta = mjs.multiply(beta, 180 / PI);
gamma = mjs.multiply(gamma, 180 / PI);
}
return [alpha, beta, gamma];
}
// method wrapping around the equivalentEuler function
equivalentEuler(mode = 'zyz', active = true, eval_convention = null, degrees = false) {
let eulers = this.euler(mode, active, eval_convention, degrees);
return equivalentEuler(eulers, active, degrees);
}
relativeOrientationTo(otherTensor, specification = {}) {
return new RelativeTensorOrientation(this, otherTensor, specification);
}
/**
* Convert this TensorData to return a clone that has been converted from
* atomic units to Hertz, assuming it's an Electric Field Gradient tensor.
*
* @param {Number} Q Quadrupolar moment of the given nucleus
* (millibarn, as in the isotope data tables)
*
* @return {TensorData} Converted tensor
*/
efgAtomicToHz(Q) {
// Clone self, then multiply all the necessary quantities
var clone = _.clone(this);
var k = efg2hz*Q;
clone._M = mjs.multiply(this._M, k);
clone._Msymm = mjs.multiply(this._Msymm, k);
clone._Masymm = mjs.multiply(this._Masymm, k);
clone._iso = this._iso*k;
clone._evals = mjs.multiply(this._evals, k);
clone._haeb_evals = mjs.multiply(this._haeb_evals, k);
clone._nqr_evals = mjs.multiply(this._nqr_evals, k);
return clone;
}
/**
* Convert this TensorData to return a clone that has been converted from
* atomic units to Hertz, assuming it's an Indirect Spin-spin Coupling
* tensor.
*
* @param {Number} g1 Gyromagnetic ratio of the first atom
* @param {Number} g2 Gyromagnetic ratio of the second atom
*
* @return {TensorData} Converted tensor
*/
iscAtomicToHz(g1, g2) {
// Clone self, then multiply all the necessary quantities
var clone = _.clone(this);
var k = isc2hz*g1*g2;
clone._M = mjs.multiply(this._M, k);
clone._Msymm = mjs.multiply(this._Msymm, k);
clone._Masymm = mjs.multiply(this._Masymm, k);
clone._iso = this._iso*k;
clone._evals = mjs.multiply(this._evals, k);
clone._haeb_evals = mjs.multiply(this._haeb_evals, k);
clone._nqr_evals = mjs.multiply(this._nqr_evals, k);
return clone;
}
/**
* Get the eigenvalues of the tensor, sorted in the given convention.
*
* @param {String} convention Convention to sort the eigenvalues by.
* Can be "increasing", "decreasing",
* "haeberlen" or "nqr".
* @param {Boolean} return_indices If true, return the indices of the
* sorted eigenvalues, in addition to the
* sorted eigenvalues.
*
* @return {Array} Array of sorted eigenvalues, or array of sorted
* eigenvalues and indices, depending on the value of
* return_indices.
*
* @throws {Error} If the convention is not one of the above.
*
* @example
* // Get the eigenvalues of a tensor, sorted in increasing order
* var tensor = new TensorData([[1, 0, 0], [0, 2, 0], [0, 0, -6]]);
* var sorted = tensor.sorted_eigenvalues("increasing");
* // sorted = [-6, 1, 2]
*
*/
sorted_eigenvalues(convention, return_indices=false) {
const EPS = 1e-6;
let evals = this._orig_evals;
let iso = this._iso;
let sort_i = [];
if (mjs.norm(mjs.subtract(evals, iso)) < EPS) {
// isotropic tensor (all eigenvalues are the same)
sort_i = [0, 1, 2];
// TODO: what about the case where two eigenvalues are the same?
} else {
convention = convention.toLowerCase();
switch (convention) {
case "increasing":
case "i":
sort_i = _.sortBy(_.range(3), x => evals[x]);
break;
case "decreasing":
case "d":
sort_i = _.sortBy(_.range(3), x => evals[x]).reverse();
break;
case "haeberlen":
case "h":
sort_i = _.sortBy(_.range(3), x => mjs.abs(evals[x] - iso));
[sort_i[0], sort_i[1]] = [sort_i[1], sort_i[0]];
break;
case "nqr":
case "n":
sort_i = _.sortBy(_.range(3), x => mjs.abs(evals[x] - iso));
break;
default:
throw new Error("Unknown eigenvalue sorting convention: " + convention);
}
}
let sorted_evals = sort_i.map(i => this._evals[i]);
return return_indices ? [sorted_evals, sort_i] : sorted_evals;
}
/**
* Get the eigenvectors of the tensor, sorted in the given convention.
*
* @param {String} convention Convention to sort the eigenvalues by.
* Can be "increasing", "decreasing",
* "haeberlen" or "nqr".
*
* @return {Array} Array of sorted eigenvectors.
*/
sorted_eigenvectors(convention) {
// use this.sorted_eigenvalues to get the indices
let [, sort_i] = this.sorted_eigenvalues(convention, true);
let evecs = mjs.transpose(this._orig_evecs);
// make right-handed
let sorted_evecs = mjs.transpose([
evecs[sort_i[0]],
evecs[sort_i[1]],
mjs.cross(evecs[sort_i[0]], evecs[sort_i[1]])
]);
return sorted_evecs;
}
}
function eulerFromU(U, mode='zyz', active=true, EPS=1e-6) {
let alpha = 0;
let beta = 0;
let gamma = 0;
if (!active) {
U = mjs.inv(U);
}
let cos_beta = U[2][2];
// Fix for the occasional numerical error
cos_beta = Math.min(Math.max(cos_beta, -1), 1);
beta = mjs.acos(cos_beta);
if (mode.toLowerCase() === 'zyz') {
if (mjs.abs(cos_beta - 1) < EPS) {
console.warn()
// beta = 0
alpha = mjs.acos(U[0][0]);
gamma = 0;
} else {
alpha = mjs.atan2(U[1][2] / mjs.sin(beta), U[0][2] / mjs.sin(beta));
gamma = mjs.atan2(U[2][1] / mjs.sin(beta), -U[2][0] / mjs.sin(beta));
}
}
else if (mode.toLowerCase() === 'zxz') {
if (mjs.abs(cos_beta - 1) < EPS) {
// beta = 0
alpha = mjs.acos(U[0][0]);
gamma = 0;
} else {
alpha = mjs.atan2(U[0][2] / mjs.sin(beta), -U[1][2] / mjs.sin(beta));
gamma = mjs.atan2(U[2][0] / mjs.sin(beta), U[2][1] / mjs.sin(beta));
}
}
else {
throw new Error('Only zyz and zxz modes are implemented so far');
}
return [alpha, beta, gamma];
}
function normaliseEulerAngles(eulerAngles, active = true, eps = 1e-6) {
const passive = !active;
let [alpha, beta, gamma] = eulerAngles;
// wrap any negative angles
alpha = mjs.mod(alpha, 2 * PI);
beta = mjs.mod(beta, 2 * PI);
gamma = mjs.mod(gamma, 2 * PI);
if (passive) {
if (beta > PI) {
beta = 2 * PI - beta;
gamma = gamma - PI;
gamma = mjs.mod(gamma, 2 * PI);
}
if (beta >= PI / 2 - eps) {
alpha = PI - alpha;
alpha = mjs.mod(alpha, 2 * PI);
beta = PI - beta;
beta = mjs.mod(beta, 2 * PI);
gamma = PI + gamma;
gamma = mjs.mod(gamma, 2 * PI);
}
if (alpha >= PI - eps) alpha = alpha - PI;
} else {
if (beta > PI) {
beta = 2 * PI - beta;
alpha = alpha - PI;
alpha = mjs.mod(alpha, 2 * PI);
}
if (beta >= PI / 2 - eps) {
alpha = alpha + PI;
alpha = mjs.mod(alpha, 2 * PI);
beta = PI - beta;
beta = mjs.mod(beta, 2 * PI);
gamma = PI - gamma;
gamma = mjs.mod(gamma, 2 * PI);
}
if (gamma >= PI - eps) gamma = gamma - PI;
}
return [alpha, beta, gamma];
}
function handleEulerEdgeCases(eulerAngles, eigenvalues, originalTensor, mode = 'zyz', active = false, eps = 1e-6) {
const A = [...originalTensor];
const passive = !active;
const [e1, e2, e3] = eigenvalues;
const degeneracy = eigenvalues.filter(value => mjs.abs(value - e1) < eps).length;
if (degeneracy === 1) return normaliseEulerAngles(eulerAngles, active, eps);
if (degeneracy === 3) return [0, 0, 0];
if (degeneracy !== 2) throw new Error('Degeneracy must be 1, 2, or 3.');
if (mjs.abs(e1 - e2) < eps) {
// the rotation about the unique axis is free: last angle for active, first for passive
eulerAngles[passive ? 0 : 2] = 0;
} else if (mjs.abs(e2 - e3) < eps) {
// unreachable for sorted eigenvalues (a unique e1 gives degeneracy 1 above); kept as a safeguard
eulerAngles[0] = 0;
if (mode.toLowerCase() === 'zyz') {
const gamma = mjs.abs(mjs.asin(mjs.sqrt((A[1][1] - e2) / (e1 - e2))));
let beta;
if (mjs.abs(gamma - PI / 2) < eps) {
beta = 0;
} else if (mjs.abs(A[1][2]) < eps && mjs.abs(A[0][1]) < eps) {
beta = mjs.abs(mjs.asin(mjs.sqrt((A[2][2] - e3) / (e1 - e3 + (e2 - e1) * mjs.sin(gamma) ** 2))));
} else {
beta = mjs.atan2(-A[1][2] / (mjs.sin(gamma) * mjs.cos(gamma) * (e1 - e2)), A[0][1] / (mjs.sin(gamma) * mjs.cos(gamma) * (e1 - e2)));
}
eulerAngles[1] = beta;
eulerAngles[2] = gamma;
} else if (mode.toLowerCase() === 'zxz') {
const alpha = mjs.abs(mjs.asin(mjs.sqrt((A[0][0] - e2) / (e1 - e2))));
let beta;
if (mjs.abs(alpha - PI / 2) < eps) {
beta = 0;
} else if (mjs.abs(A[0][1]) < eps && mjs.abs(A[1][2]) < eps) {
beta = mjs.abs(mjs.asin(mjs.sqrt((A[2][2] - e3) / (e1 - e3 + (e2 - e1) * mjs.sin(alpha) ** 2))));
} else {
beta = mjs.atan2(-A[0][1] / (mjs.sin(alpha) * mjs.cos(alpha) * (e1 - e2)), A[1][2] / (mjs.sin(alpha) * mjs.cos(alpha) * (e1 - e2)));
}
eulerAngles[0] = alpha;
eulerAngles[1] = beta;
}
if (passive) {
const [alpha, beta, gamma] = eulerAngles;
eulerAngles = normaliseEulerAngles([-gamma, -beta, -alpha], active, eps);
}
} else {
throw new Error(`Unexpected degeneracy when computing Euler angles. Eigenvalues are ordered: ${eigenvalues}`);
}
return normaliseEulerAngles(eulerAngles, active, eps);
}
/**
* Calculates equivalent Euler angles for a given set of Euler angles.
*
* @param {Array} eulerAngles - An array of three Euler angles [alpha, beta, gamma] in **radians**.
* @param {boolean} [active=true] - Whether the rotation is active or passive.
* @param {boolean} [degrees=false] - Whether to return the angles in degrees. If false, the angles are returned in radians.
*
* @returns {Array} An array of four sets of equivalent Euler angles. Each set of angles is an array of three elements.
*
* The first set of angles is the original Euler angles. The remaining sets are calculated based on the convention and whether the rotation is active or passive.
*
* If degrees is true, all angles are converted to degrees.
*/
function equivalentEuler(eulerAngles, active = true, degrees = false) {
// the order of these doesn't really matter, but has been chosen to match
// that in the TensorView for MATLAB code
// (which is different to that in the corresponding paper)
const passive = !active;
let equivAngles = mjs.zeros([4, 3]);
let [alpha, beta, gamma] = eulerAngles;
// set the first row of the array to the original Euler angles
equivAngles[0] = [alpha, beta, gamma];
if (passive) {
equivAngles[1] = [PI + alpha, beta, gamma];
equivAngles[2] = [PI - alpha, PI - beta, PI + gamma];
equivAngles[3] = [2 * PI - alpha, PI - beta, PI + gamma];
} else {
equivAngles[1] = [alpha, beta, PI + gamma];
equivAngles[2] = [PI + alpha, PI - beta, PI - gamma];
equivAngles[3] = [PI + alpha, PI - beta, 2 * PI - gamma];
}
// now wrap any negative angles in equivAngles or > 2pi
equivAngles = equivAngles.map(value => value.map(value => mjs.mod(value, 2 * PI)));
// convert to degrees if necessary
if (degrees) {
equivAngles = equivAngles.map(value => value.map(value => mjs.multiply(value, 180 / PI)));
}
return equivAngles;
}
function convertToMatrix(M) {
if (M instanceof TensorData) {
return M._M;
} else if (Array.isArray(M)) {
if (M.length === 3 && M[0].length === 3) {
// Convert 3x3 array to mathjs matrix
return mjs.matrix(M);
} else if (M.length === 9) {
// Convert 9-element array to 3x3 mathjs matrix
return mjs.matrix([M.slice(0, 3), M.slice(3, 6), M.slice(6, 9)]);
} else if (M.length === 2) {
// Assume it's the eigenvalues and eigenvectors
const eigs = M[0];
const evecs = M[1];
M = mjs.multiply(mjs.multiply(evecs, mjs.diag(eigs)), mjs.transpose(evecs));
return mjs.matrix(M);
} else {
throw new Error('Array must be 3x3 or 9-element or [eigenvalues, eigenvectors]');
}
} else if (M instanceof mjs.Matrix) {
if (M.size()[0] !== 3 || M.size()[1] !== 3) {
throw new Error('Matrix must be 3x3');
}
return M;
} else if (typeof M === 'object' && M.eigenvalues && M.eigenvectors) {
const eigs = M.eigenvalues;
const evecs = M.eigenvectors;
M = mjs.multiply(mjs.multiply(evecs, mjs.diag(eigs)), mjs.transpose(evecs));
return mjs.matrix(M);
} else {
throw new Error('Invalid input type');
}
}
export {
equivalentEuler,
eulerFromU,
rotateTensor,
TensorData,
convertToMatrix,
}