SymOp

// ******************************************************************************
//
// Title:       Force Field X.
// Description: Force Field X - Software for Molecular Biophysics.
// Copyright:   Copyright (c) Michael J. Schnieders 2001-2023.
//
// This file is part of Force Field X.
//
// Force Field X is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License version 3 as published by
// the Free Software Foundation.
//
// Force Field X 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 GNU General Public License for more
// details.
//
// You should have received a copy of the GNU General Public License along with
// Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
// Place, Suite 330, Boston, MA 02111-1307 USA
//
// Linking this library statically or dynamically with other modules is making a
// combined work based on this library. Thus, the terms and conditions of the
// GNU General Public License cover the whole combination.
//
// As a special exception, the copyright holders of this library give you
// permission to link this library with independent modules to produce an
// executable, regardless of the license terms of these independent modules, and
// to copy and distribute the resulting executable under terms of your choice,
// provided that you also meet, for each linked independent module, the terms
// and conditions of the license of that module. An independent module is a
// module which is not derived from or based on this library. If you modify this
// library, you may extend this exception to your version of the library, but
// you are not obligated to do so. If you do not wish to do so, delete this
// exception statement from your version.
//
// ******************************************************************************
package ffx.crystal;

import static ffx.numerics.math.DoubleMath.dot;
import static ffx.numerics.math.MatrixMath.mat3Inverse;
import static ffx.numerics.math.MatrixMath.mat4Mat4;
import static ffx.numerics.math.ScalarMath.mod;
import static java.lang.Double.parseDouble;
import static java.lang.String.format;
import static org.apache.commons.math3.util.FastMath.PI;
import static org.apache.commons.math3.util.FastMath.cos;
import static org.apache.commons.math3.util.FastMath.random;
import static org.apache.commons.math3.util.FastMath.rint;
import static org.apache.commons.math3.util.FastMath.sin;
import static org.apache.commons.math3.util.FastMath.sqrt;

/**
 * The SymOp class defines the rotation and translation of a single symmetry operator.
 *
 * @author Michael J. Schnieders
 * @see SpaceGroup
 * @since 1.0
 */
public class SymOp {

  /** Constant <code>ZERO = 0.0</code> */
  private static final double ZERO = 0.0;
  /** Constant <code>ZERO_ROTATION = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}</code> */
  public static final double[][] ZERO_ROTATION = {
      {1.0, ZERO, ZERO},
      {ZERO, 1.0, ZERO},
      {ZERO, ZERO, 1.0}};
  /** Constant <code>Tr_0_0_0={ZERO, ZERO, ZERO}</code> */
  public static final double[] Tr_0_0_0 = {ZERO, ZERO, ZERO};
  /** Constant <code>f12 = 1.0 / 2.0</code> */
  private static final double f12 = 1.0 / 2.0;
  /** Constant <code>Tr_12_0_12={f12, ZERO, f12}</code> */
  static final double[] Tr_12_0_12 = {f12, ZERO, f12};
  /** Constant <code>Tr_0_12_12={ZERO, f12, f12}</code> */
  static final double[] Tr_0_12_12 = {ZERO, f12, f12};
  /** Constant <code>Tr_12_12_0={f12, f12, ZERO}</code> */
  static final double[] Tr_12_12_0 = {f12, f12, ZERO};
  /** Constant <code>Tr_12_12_12={f12, f12, f12}</code> */
  static final double[] Tr_12_12_12 = {f12, f12, f12};
  /** Constant <code>Tr_0_12_0={ZERO, f12, ZERO}</code> */
  static final double[] Tr_0_12_0 = {ZERO, f12, ZERO};
  /** Constant <code>Tr_12_0_0={f12, ZERO, ZERO}</code> */
  static final double[] Tr_12_0_0 = {f12, ZERO, ZERO};
  /** Constant <code>Tr_0_0_12={ZERO, ZERO, f12}</code> */
  static final double[] Tr_0_0_12 = {ZERO, ZERO, f12};
  /** Constant <code>f13 = 1.0 / 3.0</code> */
  private static final double f13 = 1.0 / 3.0;
  /** Constant <code>Tr_0_0_13={ZERO, ZERO, f13}</code> */
  static final double[] Tr_0_0_13 = {ZERO, ZERO, f13};
  /** Constant <code>f23 = 2.0 / 3.0</code> */
  private static final double f23 = 2.0 / 3.0;
  /** Constant <code>Tr_23_13_13={f23, f13, f13}</code> */
  static final double[] Tr_23_13_13 = {f23, f13, f13};
  /** Constant <code>Tr_13_23_23={f13, f23, f23}</code> */
  static final double[] Tr_13_23_23 = {f13, f23, f23};
  /** Constant <code>Tr_0_0_23={ZERO, ZERO, f23}</code> */
  static final double[] Tr_0_0_23 = {ZERO, ZERO, f23};
  /** Constant <code>f14 = 1.0 / 4.0</code> */
  private static final double f14 = 1.0 / 4.0;
  /** Constant <code>Tr_12_0_14={f12, ZERO, f14}</code> */
  static final double[] Tr_12_0_14 = {f12, ZERO, f14};
  /** Constant <code>Tr_0_12_14={ZERO, f12, f14}</code> */
  static final double[] Tr_0_12_14 = {ZERO, f12, f14};
  /** Constant <code>Tr_14_14_14={f14, f14, f14}</code> */
  static final double[] Tr_14_14_14 = {f14, f14, f14};
  /** Constant <code>Tr_12_12_14={f12, f12, f14}</code> */
  static final double[] Tr_12_12_14 = {f12, f12, f14};
  /** Constant <code>Tr_0_0_14={ZERO, ZERO, f14}</code> */
  static final double[] Tr_0_0_14 = {ZERO, ZERO, f14};
  /** Constant <code>f34 = 3.0 / 4.0</code> */
  private static final double f34 = 3.0 / 4.0;
  /** Constant <code>Tr_0_0_34={ZERO, ZERO, f34}</code> */
  static final double[] Tr_0_0_34 = {ZERO, ZERO, f34};
  /** Constant <code>Tr_12_0_34={f12, ZERO, f34}</code> */
  static final double[] Tr_12_0_34 = {f12, ZERO, f34};
  /** Constant <code>Tr_0_12_34={ZERO, f12, f34}</code> */
  static final double[] Tr_0_12_34 = {ZERO, f12, f34};
  /** Constant <code>Tr_34_14_14={f34, f14, f14}</code> */
  static final double[] Tr_34_14_14 = {f34, f14, f14};
  /** Constant <code>Tr_14_14_34={f14, f14, f34}</code> */
  static final double[] Tr_14_14_34 = {f14, f14, f34};
  /** Constant <code>Tr_14_34_14={f14, f34, f14}</code> */
  static final double[] Tr_14_34_14 = {f14, f34, f14};
  /** Constant <code>Tr_12_12_34={f12, f12, f34}</code> */
  static final double[] Tr_12_12_34 = {f12, f12, f34};
  /** Constant <code>Tr_14_34_34={f14, f34, f34}</code> */
  static final double[] Tr_14_34_34 = {f14, f34, f34};
  /** Constant <code>Tr_34_34_14={f34, f34, f14}</code> */
  static final double[] Tr_34_34_14 = {f34, f34, f14};
  /** Constant <code>Tr_34_34_34={f34, f34, f34}</code> */
  static final double[] Tr_34_34_34 = {f34, f34, f34};
  /** Constant <code>Tr_34_14_34={f34, f14, f34}</code> */
  static final double[] Tr_34_14_34 = {f34, f14, f34};
  /** Constant <code>f16 = 1.0 / 6.0</code> */
  private static final double f16 = 1.0 / 6.0;
  /** Constant <code>Tr_13_23_16={f13, f23, f16}</code> */
  static final double[] Tr_13_23_16 = {f13, f23, f16};
  /** Constant <code>Tr_0_0_16={ZERO, ZERO, f16}</code> */
  static final double[] Tr_0_0_16 = {ZERO, ZERO, f16};
  /** Constant <code>f56 = 5.0 / 6.0</code> */
  private static final double f56 = 5.0 / 6.0;
  /** Constant <code>Tr_0_0_56={ZERO, ZERO, f56}</code> */
  static final double[] Tr_0_0_56 = {ZERO, ZERO, f56};
  /** Constant <code>Tr_23_13_56={f23, f13, f56}</code> */
  static final double[] Tr_23_13_56 = {f23, f13, f56};
  /** Constant <code>X = {1.0, ZERO, ZERO}</code> */
  private static final double[] X = {1.0, ZERO, ZERO};
  /** Constant <code>Y = {ZERO, 1.0, ZERO}</code> */
  private static final double[] Y = {ZERO, 1.0, ZERO};
  /** Constant <code>Z = {ZERO, ZERO, 1.0}</code> */
  private static final double[] Z = {ZERO, ZERO, 1.0};
  /** Constant <code>Rot_Y_Z_X={Y, Z, X}</code> */
  static final double[][] Rot_Y_Z_X = {Y, Z, X};
  /** Constant <code>Rot_X_Y_Z={X, Y, Z}</code> */
  static final double[][] Rot_X_Y_Z = {X, Y, Z};
  /** Constant <code>Rot_Z_X_Y={Z, X, Y}</code> */
  static final double[][] Rot_Z_X_Y = {Z, X, Y};
  /** Constant <code>Rot_X_Z_Y={X, Z, Y}</code> */
  static final double[][] Rot_X_Z_Y = {X, Z, Y};
  /** Constant <code>Rot_Z_Y_X={Z, Y, X}</code> */
  static final double[][] Rot_Z_Y_X = {Z, Y, X};
  /** Constant <code>Rot_Y_X_Z={Y, X, Z}</code> */
  static final double[][] Rot_Y_X_Z = {Y, X, Z};
  /** Constant <code>mX = {-1.0, ZERO, ZERO}</code> */
  private static final double[] mX = {-1.0, ZERO, ZERO};
  /** Constant <code>Rot_Y_mX_Z={Y, mX, Z}</code> */
  static final double[][] Rot_Y_mX_Z = {Y, mX, Z};
  /** Constant <code>Rot_mX_Z_Y={mX, Z, Y}</code> */
  static final double[][] Rot_mX_Z_Y = {mX, Z, Y};
  /** Constant <code>Rot_Y_Z_mX={Y, Z, mX}</code> */
  static final double[][] Rot_Y_Z_mX = {Y, Z, mX};
  /** Constant <code>Rot_mX_Y_Z={mX, Y, Z}</code> */
  static final double[][] Rot_mX_Y_Z = {mX, Y, Z};
  /** Constant <code>Rot_Z_Y_mX={Z, Y, mX}</code> */
  static final double[][] Rot_Z_Y_mX = {Z, Y, mX};
  /** Constant <code>Rot_Z_mX_Y={Z, mX, Y}</code> */
  static final double[][] Rot_Z_mX_Y = {Z, mX, Y};
  /** Constant <code>mY = {ZERO, -1.0, ZERO}</code> */
  private static final double[] mY = {ZERO, -1.0, ZERO};
  /** Constant <code>Rot_Z_mY_X={Z, mY, X}</code> */
  static final double[][] Rot_Z_mY_X = {Z, mY, X};
  /** Constant <code>Rot_X_Z_mY={X, Z, mY}</code> */
  static final double[][] Rot_X_Z_mY = {X, Z, mY};
  /** Constant <code>Rot_mY_X_Z={mY, X, Z}</code> */
  static final double[][] Rot_mY_X_Z = {mY, X, Z};
  /** Constant <code>Rot_mY_Z_mX={mY, Z, mX}</code> */
  static final double[][] Rot_mY_Z_mX = {mY, Z, mX};
  /** Constant <code>Rot_mY_Z_X={mY, Z, X}</code> */
  static final double[][] Rot_mY_Z_X = {mY, Z, X};
  /** Constant <code>Rot_Z_X_mY={Z, X, mY}</code> */
  static final double[][] Rot_Z_X_mY = {Z, X, mY};
  /** Constant <code>Rot_Z_mX_mY={Z, mX, mY}</code> */
  static final double[][] Rot_Z_mX_mY = {Z, mX, mY};
  /** Constant <code>Rot_mX_Z_mY={mX, Z, mY}</code> */
  static final double[][] Rot_mX_Z_mY = {mX, Z, mY};
  /** Constant <code>Rot_X_mY_Z={X, mY, Z}</code> */
  static final double[][] Rot_X_mY_Z = {X, mY, Z};
  /** Constant <code>Rot_mY_mX_Z={mY, mX, Z}</code> */
  static final double[][] Rot_mY_mX_Z = {mY, mX, Z};
  /** Constant <code>Rot_Z_mY_mX={Z, mY, mX}</code> */
  static final double[][] Rot_Z_mY_mX = {Z, mY, mX};
  /** Constant <code>Rot_mX_mY_Z={mX, mY, Z}</code> */
  static final double[][] Rot_mX_mY_Z = {mX, mY, Z};
  /** Constant <code>mZ = {ZERO, ZERO, -1.0}</code> */
  private static final double[] mZ = {ZERO, ZERO, -1.0};
  /** Constant <code>Rot_Y_mX_mZ={Y, mX, mZ}</code> */
  static final double[][] Rot_Y_mX_mZ = {Y, mX, mZ};
  /** Constant <code>Rot_mX_Y_mZ={mX, Y, mZ}</code> */
  static final double[][] Rot_mX_Y_mZ = {mX, Y, mZ};
  /** Constant <code>Rot_X_mZ_Y={X, mZ, Y}</code> */
  static final double[][] Rot_X_mZ_Y = {X, mZ, Y};
  /** Constant <code>Rot_mY_mZ_X={mY, mZ, X}</code> */
  static final double[][] Rot_mY_mZ_X = {mY, mZ, X};
  /** Constant <code>Rot_Y_X_mZ={Y, X, mZ}</code> */
  static final double[][] Rot_Y_X_mZ = {Y, X, mZ};
  /** Constant <code>Rot_Y_mZ_X={Y, mZ, X}</code> */
  static final double[][] Rot_Y_mZ_X = {Y, mZ, X};
  /** Constant <code>Rot_mX_mY_mZ={mX, mY, mZ}</code> */
  static final double[][] Rot_mX_mY_mZ = {mX, mY, mZ};
  /** Constant <code>Rot_X_Y_mZ={X, Y, mZ}</code> */
  static final double[][] Rot_X_Y_mZ = {X, Y, mZ};
  /** Constant <code>Rot_mZ_mY_mX={mZ, mY, mX}</code> */
  static final double[][] Rot_mZ_mY_mX = {mZ, mY, mX};
  /** Constant <code>Rot_X_mZ_mY={X, mZ, mY}</code> */
  static final double[][] Rot_X_mZ_mY = {X, mZ, mY};
  /** Constant <code>Rot_mY_mX_mZ={mY, mX, mZ}</code> */
  static final double[][] Rot_mY_mX_mZ = {mY, mX, mZ};
  /** Constant <code>Rot_mY_X_mZ={mY, X, mZ}</code> */
  static final double[][] Rot_mY_X_mZ = {mY, X, mZ};
  /** Constant <code>Rot_mX_mZ_mY={mX, mZ, mY}</code> */
  static final double[][] Rot_mX_mZ_mY = {mX, mZ, mY};
  /** Constant <code>Rot_mZ_mX_mY={mZ, mX, mY}</code> */
  static final double[][] Rot_mZ_mX_mY = {mZ, mX, mY};
  /** Constant <code>Rot_mZ_mY_X={mZ, mY, X}</code> */
  static final double[][] Rot_mZ_mY_X = {mZ, mY, X};
  /** Constant <code>Rot_mZ_Y_mX={mZ, Y, mX}</code> */
  static final double[][] Rot_mZ_Y_mX = {mZ, Y, mX};
  /** Constant <code>Rot_mZ_mX_Y={mZ, mX, Y}</code> */
  static final double[][] Rot_mZ_mX_Y = {mZ, mX, Y};
  /** Constant <code>Rot_mX_mZ_Y={mX, mZ, Y}</code> */
  static final double[][] Rot_mX_mZ_Y = {mX, mZ, Y};
  /** Constant <code>Rot_X_mY_mZ={X, mY, mZ}</code> */
  static final double[][] Rot_X_mY_mZ = {X, mY, mZ};
  /** Constant <code>Rot_mZ_X_Y={mZ, X, Y}</code> */
  static final double[][] Rot_mZ_X_Y = {mZ, X, Y};
  /** Constant <code>Rot_Y_mZ_mX={Y, mZ, mX}</code> */
  static final double[][] Rot_Y_mZ_mX = {Y, mZ, mX};
  /** Constant <code>Rot_mY_mZ_mX={mY, mZ, mX}</code> */
  static final double[][] Rot_mY_mZ_mX = {mY, mZ, mX};
  /** Constant <code>Rot_mZ_Y_X={mZ, Y, X}</code> */
  static final double[][] Rot_mZ_Y_X = {mZ, Y, X};
  /** Constant <code>Rot_mZ_X_mY={mZ, X, mY}</code> */
  static final double[][] Rot_mZ_X_mY = {mZ, X, mY};
  /** Constant <code>XmY = {1.0, -1.0, ZERO}</code> */
  private static final double[] XmY = {1.0, -1.0, ZERO};
  /** Constant <code>Rot_XmY_X_mZ={XmY, X, mZ}</code> */
  static final double[][] Rot_XmY_X_mZ = {XmY, X, mZ};
  /** Constant <code>Rot_XmY_X_Z={XmY, X, Z}</code> */
  static final double[][] Rot_XmY_X_Z = {XmY, X, Z};
  /** Constant <code>Rot_XmY_mY_Z={XmY, mY, Z}</code> */
  static final double[][] Rot_XmY_mY_Z = {XmY, mY, Z};
  /** Constant <code>Rot_X_XmY_Z={X, XmY, Z}</code> */
  static final double[][] Rot_X_XmY_Z = {X, XmY, Z};
  /** Constant <code>Rot_X_XmY_mZ={X, XmY, mZ}</code> */
  static final double[][] Rot_X_XmY_mZ = {X, XmY, mZ};
  /** Constant <code>Rot_mY_XmY_mZ={mY, XmY, mZ}</code> */
  static final double[][] Rot_mY_XmY_mZ = {mY, XmY, mZ};
  /** Constant <code>Rot_mY_XmY_Z={mY, XmY, Z}</code> */
  static final double[][] Rot_mY_XmY_Z = {mY, XmY, Z};
  /** Constant <code>Rot_XmY_mY_mZ={XmY, mY, mZ}</code> */
  static final double[][] Rot_XmY_mY_mZ = {XmY, mY, mZ};
  /** Constant <code>mXY = {-1.0, 1.0, ZERO}</code> */
  private static final double[] mXY = {-1.0, 1.0, ZERO};
  /** Constant <code>Rot_Y_mXY_Z={Y, mXY, Z}</code> */
  static final double[][] Rot_Y_mXY_Z = {Y, mXY, Z};
  /** Constant <code>Rot_mX_mXY_mZ={mX, mXY, mZ}</code> */
  static final double[][] Rot_mX_mXY_mZ = {mX, mXY, mZ};
  /** Constant <code>Rot_mXY_Y_Z={mXY, Y, Z}</code> */
  static final double[][] Rot_mXY_Y_Z = {mXY, Y, Z};
  /** Constant <code>Rot_mXY_mX_Z={mXY, mX, Z}</code> */
  static final double[][] Rot_mXY_mX_Z = {mXY, mX, Z};
  /** Constant <code>Rot_mXY_Y_mZ={mXY, Y, mZ}</code> */
  static final double[][] Rot_mXY_Y_mZ = {mXY, Y, mZ};
  /** Constant <code>Rot_mXY_mX_mZ={mXY, mX, mZ}</code> */
  static final double[][] Rot_mXY_mX_mZ = {mXY, mX, mZ};
  /** Constant <code>Rot_mX_mXY_Z={mX, mXY, Z}</code> */
  static final double[][] Rot_mX_mXY_Z = {mX, mXY, Z};
  /** Constant <code>Rot_Y_mXY_mZ={Y, mXY, mZ}</code> */
  static final double[][] Rot_Y_mXY_mZ = {Y, mXY, mZ};
  /** The rotation matrix in fractional coordinates. */
  public final double[][] rot;
  /** The translation vector in fractional coordinates. */
  public final double[] tr;
  /** A mask equal to 0 for X-coordinates. */
  private static final int XX = 0;
  /** A mask equal to 1 for Y-coordinates. */
  private static final int YY = 1;
  /** A mask equal to 2 for Z-coordinates. */
  private static final int ZZ = 2;
  /** Replicates position for the given symmetry operator (LxMxN). */
  public final int[] replicatesVector;

  /**
   * The SymOp constructor using a rotation matrix and translation vector.
   *
   * @param rot The rotation matrix.
   * @param tr The translation vector.
   */
  public SymOp(double[][] rot, double[] tr) {
    this.rot = rot;
    this.tr = tr;
    replicatesVector = new int[]{0, 0, 0};
  }

  /**
   * The SymOp constructor using a rotation matrix and translation vector.
   *
   * @param rot The rotation matrix.
   * @param tr The translation vector.
   * @param replicatesVector Describes symmetry operators location within replicates crystal.
   */
  public SymOp(double[][] rot, double[] tr, int[] replicatesVector) {
    this.rot = rot;
    this.tr = tr;
    this.replicatesVector = replicatesVector;
  }


  /**
   * The SymOp constructor using a 4x4 matrix.
   *
   * @param m The rotation matrix and translation vector as a 4x4 matrix.
   */
  public SymOp(double[][] m) {
    this.rot = new double[3][3];
    rot[0][0] = m[0][0];
    rot[0][1] = m[0][1];
    rot[0][2] = m[0][2];
    rot[1][0] = m[1][0];
    rot[1][1] = m[1][1];
    rot[1][2] = m[1][2];
    rot[2][0] = m[2][0];
    rot[2][1] = m[2][1];
    rot[2][2] = m[2][2];

    this.tr = new double[3];
    tr[0] = m[0][3] / m[3][3];
    tr[1] = m[1][3] / m[3][3];
    tr[2] = m[2][3] / m[3][3];

    replicatesVector = new int[]{0, 0, 0};
  }

  /**
   * symPhaseShift
   *
   * @param hkl an array of double.
   * @return a double.
   */
  public double symPhaseShift(double[] hkl) {
    // Apply translation
    return -2.0 * PI * (hkl[0] * tr[0] + hkl[1] * tr[1] + hkl[2] * tr[2]);
  }

  /**
   * symPhaseShift
   *
   * @param hkl a {@link HKL} object.
   * @return a double.
   */
  public double symPhaseShift(HKL hkl) {
    // Apply translation
    return -2.0 * PI * (hkl.getH() * tr[0] + hkl.getK() * tr[1] + hkl.getL() * tr[2]);
  }

  /**
   * Return the SymOp as a 4x4 matrix.
   *
   * @return A 4x4 matrix representation of the SymOp.
   */
  public double[][] asMatrix() {
    return new double[][] {
        {rot[0][0], rot[0][1], rot[0][2], tr[0]},
        {rot[1][0], rot[1][1], rot[1][2], tr[1]},
        {rot[2][0], rot[2][1], rot[2][2], tr[2]},
        {0.0, 0.0, 0.0, 1.0}};
  }

  /**
   * Return the combined SymOp that is equivalent to first applying <code>this</code> SymOp and then
   * the argument. Note: Applied as rotation then translation.
   * <code>X' = S_arg(S_this(X))</code>
   * <code>X' = S_combined(X)</code>
   *
   * @param symOp The SymOp to append to <code>this</code> SymOp.
   * @return The combined SymOp.
   */
  public SymOp append(SymOp symOp) {
    return new SymOp(mat4Mat4(symOp.asMatrix(), asMatrix()));
  }

  /**
   * Return the combined SymOp that is equivalent to first applying the argument and then
   * <code>this</code> SymOp. Note: Applied as rotation then translation.
   * <code>X' = S_this(S_arg(X))</code>
   * <code>X' = S_combined(X)</code>
   *
   * @param symOp The SymOp to prepend to <code>this</code> SymOp.
   * @return The combined SymOp.
   */
  public SymOp prepend(SymOp symOp) {
    return new SymOp(mat4Mat4(asMatrix(), symOp.asMatrix()));
  }

  /**
   * Return the combined SymOp that is equivalent to first applying symOp1 and then SymOp2. Note:
   * Applied as rotation then translation.
   * <p>
   * <code>X' = S_2(S_1(X))</code>
   * <code>X' = S_combined(X)</code>
   *
   * @param symOp1 The fist SymOp.
   * @param symOp2 The second SymOp.
   * @return The combined SymOp.
   */
  public static SymOp combineSymOps(SymOp symOp1, SymOp symOp2) {
    return new SymOp(mat4Mat4(symOp2.asMatrix(), symOp1.asMatrix()));
  }

  /**
   * Apply a Cartesian symmetry operator to an array of Cartesian coordinates. If the arrays x, y or
   * z are null or not of length n, the method returns immediately. If mateX, mateY or mateZ are null
   * or not of length n, the method returns immediately.
   *
   * @param n Number of atoms.
   * @param x Input cartesian x-coordinates.
   * @param y Input cartesian y-coordinates.
   * @param z Input cartesian z-coordinates.
   * @param mateX Output cartesian x-coordinates.
   * @param mateY Output cartesian y-coordinates.
   * @param mateZ Output cartesian z-coordinates.
   * @param symOp The cartesian symmetry operator.
   */
  public static void applyCartSymOp(int n, double[] x, double[] y, double[] z,
      double[] mateX, double[] mateY, double[] mateZ, SymOp symOp) {
    if (x == null || y == null || z == null) {
      throw new IllegalArgumentException("The input arrays x, y and z must not be null.");
    }
    if (x.length < n || y.length < n || z.length < n) {
      throw new IllegalArgumentException("The input arrays x, y and z must be of length n: " + n);
    }
    if (mateX == null || mateY == null || mateZ == null) {
      throw new IllegalArgumentException("The output arrays mateX, mateY and mateZ must not be null.");
    }
    if (mateX.length < n || mateY.length < n || mateZ.length < n) {
      throw new IllegalArgumentException("The output arrays mateX, mateY and mateZ must be of length n: " + n);
    }

    final double[][] rot = symOp.rot;
    final double[] trans = symOp.tr;

    final double rot00 = rot[0][0];
    final double rot10 = rot[1][0];
    final double rot20 = rot[2][0];
    final double rot01 = rot[0][1];
    final double rot11 = rot[1][1];
    final double rot21 = rot[2][1];
    final double rot02 = rot[0][2];
    final double rot12 = rot[1][2];
    final double rot22 = rot[2][2];
    final double t0 = trans[0];
    final double t1 = trans[1];
    final double t2 = trans[2];
    for (int i = 0; i < n; i++) {
      double xc = x[i];
      double yc = y[i];
      double zc = z[i];
      // Apply Symmetry Operator.
      mateX[i] = rot00 * xc + rot01 * yc + rot02 * zc + t0;
      mateY[i] = rot10 * xc + rot11 * yc + rot12 * zc + t1;
      mateZ[i] = rot20 * xc + rot21 * yc + rot22 * zc + t2;
    }
  }

  /**
   * Apply a cartesian symmetry operator to an array of coordinates.
   *
   * @param xyz Input  cartesian coordinates.
   * @param mate Symmetry mate  cartesian coordinates.
   * @param symOp The cartesian symmetry operator.
   */
  public static void applyCartesianSymOp(double[] xyz, double[] mate, SymOp symOp) {
    applyCartesianSymOp(xyz, mate, symOp, null);
  }

  /**
   * Apply a cartesian symmetry operator to an array of coordinates.
   *
   * @param xyz Input  cartesian coordinates.
   * @param mate Symmetry mate  cartesian coordinates.
   * @param symOp The cartesian symmetry operator.
   * @param mask Only apply the SymOp if the per atom mask is true.
   */
  public static void applyCartesianSymOp(double[] xyz, double[] mate, SymOp symOp, boolean[] mask) {
    var rot = symOp.rot;
    var trans = symOp.tr;

    assert (xyz.length % 3 == 0);
    assert (xyz.length == mate.length);
    int len = xyz.length / 3;

    // Load the SymOp into local variables.
    var rot00 = rot[0][0];
    var rot10 = rot[1][0];
    var rot20 = rot[2][0];
    var rot01 = rot[0][1];
    var rot11 = rot[1][1];
    var rot21 = rot[2][1];
    var rot02 = rot[0][2];
    var rot12 = rot[1][2];
    var rot22 = rot[2][2];
    var tx = trans[0];
    var ty = trans[1];
    var tz = trans[2];

    if (mask == null) {
      for (int i = 0; i < len; i++) {
        int index = i * 3;
        var xc = xyz[index + XX];
        var yc = xyz[index + YY];
        var zc = xyz[index + ZZ];
        // Apply Symmetry Operator.
        mate[index + XX] = rot00 * xc + rot01 * yc + rot02 * zc + tx;
        mate[index + YY] = rot10 * xc + rot11 * yc + rot12 * zc + ty;
        mate[index + ZZ] = rot20 * xc + rot21 * yc + rot22 * zc + tz;
      }
    } else {
      for (int i = 0; i < len; i++) {
        int index = i * 3;
        var xc = xyz[index + XX];
        var yc = xyz[index + YY];
        var zc = xyz[index + ZZ];
        if (mask[i]) {
          // Apply Symmetry Operator.
          mate[index + XX] = rot00 * xc + rot01 * yc + rot02 * zc + tx;
          mate[index + YY] = rot10 * xc + rot11 * yc + rot12 * zc + ty;
          mate[index + ZZ] = rot20 * xc + rot21 * yc + rot22 * zc + tz;
        } else {
          mate[index + XX] = xc;
          mate[index + YY] = yc;
          mate[index + ZZ] = zc;
        }
      }
    }
  }

  /**
   * Apply a fractional symmetry operator to one set of coordinates.
   *
   * @param xyz Input fractional coordinates.
   * @param mate Symmetry mate fractional coordinates.
   * @param symOp The fractional symmetry operator.
   */
  public static void applyFracSymOp(double[] xyz, double[] mate, SymOp symOp) {
    var rot = symOp.rot;
    var trans = symOp.tr;
    var xf = xyz[0];
    var yf = xyz[1];
    var zf = xyz[2];
    // Apply Symmetry Operator.
    mate[0] = rot[0][0] * xf + rot[0][1] * yf + rot[0][2] * zf + trans[0];
    mate[1] = rot[1][0] * xf + rot[1][1] * yf + rot[1][2] * zf + trans[1];
    mate[2] = rot[2][0] * xf + rot[2][1] * yf + rot[2][2] * zf + trans[2];
  }

  /**
   * Apply a symmetry operator to one set of coordinates.
   *
   * @param h Input coordinates.
   * @param k Input coordinates.
   * @param l Input coordinates.
   * @param mate Symmetry mate coordinates.
   * @param symOp The symmetry operator.
   * @param nx number of unit cell translations
   * @param ny number of unit cell translations
   * @param nz number of unit cell translations
   */
  public static void applySymOp(int h, int k, int l, int[] mate, SymOp symOp, int nx, int ny,
      int nz) {
    var rot = symOp.rot;
    var trans = symOp.tr;
    // Apply Symmetry Operator.
    mate[0] =
        (int) rot[0][0] * h + (int) rot[0][1] * k + (int) rot[0][2] * l + (int) rint(nx * trans[0]);
    mate[1] =
        (int) rot[1][0] * h + (int) rot[1][1] * k + (int) rot[1][2] * l + (int) rint(ny * trans[1]);
    mate[2] =
        (int) rot[2][0] * h + (int) rot[2][1] * k + (int) rot[2][2] * l + (int) rint(nz * trans[2]);
    mate[0] = mod(mate[0], nx);
    mate[1] = mod(mate[1], ny);
    mate[2] = mod(mate[2], nz);
  }

  /**
   * Apply a symmetry operator to one HKL.
   *
   * @param hkl Input HKL.
   * @param mate Symmetry mate HKL.
   * @param symOp The symmetry operator.
   */
  public static void applySymRot(HKL hkl, HKL mate, SymOp symOp) {
    var rot = symOp.rot;
    double h = hkl.getH();
    double k = hkl.getK();
    double l = hkl.getL();
    double hs = rot[0][0] * h + rot[0][1] * k + rot[0][2] * l;
    double ks = rot[1][0] * h + rot[1][1] * k + rot[1][2] * l;
    double ls = rot[2][0] * h + rot[2][1] * k + rot[2][2] * l;
    // Convert back to HKL
    mate.setH((int) rint(hs));
    mate.setK((int) rint(ks));
    mate.setL((int) rint(ls));
  }

  /**
   * Apply a Cartesian symmetry rotation to an array of Cartesian coordinates. The length of xyz must
   * be divisible by 3 and mate must have the same length.
   *
   * @param xyz Input cartesian x, y, z-coordinates.
   * @param mate Output cartesian x, y, z-coordinates.
   * @param symOp The fractional symmetry operator.
   */
  public static void applyCartesianSymRot(double[] xyz, double[] mate, SymOp symOp) {
    applyCartesianSymRot(xyz, mate, symOp, null);
  }

  /**
   * Apply a Cartesian symmetry rotation to an array of Cartesian coordinates. The length of xyz must
   * be divisible by 3 and mate must have the same length.
   *
   * @param xyz Input cartesian x, y, z-coordinates.
   * @param mate Output cartesian x, y, z-coordinates.
   * @param symOp The fractional symmetry operator.
   * @param mask Only apply the SymOp if the per atom mask is true.
   */
  public static void applyCartesianSymRot(double[] xyz, double[] mate, SymOp symOp, boolean[] mask) {
    int l = xyz.length;
    int n = l / 3;
    assert (l % 3 == 0);
    assert (mate.length == l);

    // Load the rotation matrix
    var rot = symOp.rot;
    var rot00 = rot[0][0];
    var rot10 = rot[1][0];
    var rot20 = rot[2][0];
    var rot01 = rot[0][1];
    var rot11 = rot[1][1];
    var rot21 = rot[2][1];
    var rot02 = rot[0][2];
    var rot12 = rot[1][2];
    var rot22 = rot[2][2];

    if (mask == null) {
      for (int i = 0; i < n; i++) {
        int index = i * 3;
        var xi = xyz[index + XX];
        var yi = xyz[index + YY];
        var zi = xyz[index + ZZ];
        // Apply Symmetry Operator.
        mate[index + XX] = rot00 * xi + rot01 * yi + rot02 * zi;
        mate[index + YY] = rot10 * xi + rot11 * yi + rot12 * zi;
        mate[index + ZZ] = rot20 * xi + rot21 * yi + rot22 * zi;
      }
    } else {
      for (int i = 0; i < n; i++) {
        int index = i * 3;
        var xi = xyz[index + XX];
        var yi = xyz[index + YY];
        var zi = xyz[index + ZZ];
        if (mask[i]) {
          // Apply Symmetry Operator.
          mate[index + XX] = rot00 * xi + rot01 * yi + rot02 * zi;
          mate[index + YY] = rot10 * xi + rot11 * yi + rot12 * zi;
          mate[index + ZZ] = rot20 * xi + rot21 * yi + rot22 * zi;
        } else {
          mate[index + XX] = xi;
          mate[index + YY] = yi;
          mate[index + ZZ] = zi;
        }
      }
    }
  }

  /**
   * Apply a transpose rotation symmetry operator to one HKL.
   *
   * @param hkl Input HKL.
   * @param mate Symmetry mate HKL.
   * @param symOp The symmetry operator.
   */
  public static void applyTransSymRot(HKL hkl, HKL mate, SymOp symOp) {
    double[][] rot = symOp.rot;
    double h = hkl.getH();
    double k = hkl.getK();
    double l = hkl.getL();
    // Apply transpose Symmetry Operator.
    double hs = rot[0][0] * h + rot[1][0] * k + rot[2][0] * l;
    double ks = rot[0][1] * h + rot[1][1] * k + rot[2][1] * l;
    double ls = rot[0][2] * h + rot[1][2] * k + rot[2][2] * l;
    // Convert back to HKL
    mate.setH((int) rint(hs));
    mate.setK((int) rint(ks));
    mate.setL((int) rint(ls));
  }

  /**
   * Print a Sym Op matrix as a continued line string.
   *
   * @param symOp Symmetry operation to print.
   * @return Continued line string.
   */
  public static String asMatrixString(SymOp symOp) {
    double[][] values = symOp.asMatrix();
    return format("""
                     %14.8f %14.8f %14.8f \\
                               %14.8f %14.8f %14.8f \\
                               %14.8f %14.8f %14.8f \\
                               %14.8f %14.8f %14.8f \
                    """, values[0][0],
            values[0][1], values[0][2], values[1][0], values[1][1], values[1][2], values[2][0],
            values[2][1], values[2][2], values[0][3], values[1][3], values[2][3]);
  }

  /**
   * Invert a symmetry operator.
   *
   * @param symOp Original symmetry operator of which the inverse is desired.
   * @return SymOp The inverse symmetry operator of the one supplied.
   */
  public static SymOp invertSymOp(SymOp symOp) {
    var tr = symOp.tr;
    var rot = symOp.rot;
    var inv = mat3Inverse(rot);
    return new SymOp(inv,
        new double[] {-dot(inv[0], tr), -dot(inv[1], tr), -dot(inv[2], tr)});
  }

  /**
   * Create a SymOp from an input String.
   *
   * @param s Input <code>String</code> containing 12 white-space delimited double values.
   * @return The SymOp.
   */
  public static SymOp parse(String s) {
    String[] tokens = s.split(" +");
    if (tokens.length < 12) {
      return null;
    }
    return new SymOp(new double[][] {
        {parseDouble(tokens[0]), parseDouble(tokens[1]), parseDouble(tokens[2])},
        {parseDouble(tokens[3]), parseDouble(tokens[4]), parseDouble(tokens[5])},
        {parseDouble(tokens[6]), parseDouble(tokens[7]), parseDouble(tokens[8])}},
        new double[] {parseDouble(tokens[9]), parseDouble(tokens[10]), parseDouble(tokens[11])});
  }

  /**
   * Generate a random Cartesian Symmetry Operator.
   *
   * @param scalar The range of translations will be from -scalar/2 .. scalar/2.
   * @return A Cartesian SymOp with a random rotation and translation.
   */
  public static SymOp randomSymOpFactory(double scalar) {
    double[] tr = {scalar * (random() - 0.5), scalar * (random() - 0.5), scalar * (random() - 0.5)};
    return randomSymOpFactory(tr);
  }

  /**
   * Generate a random Cartesian Symmetry Operator.
   *
   * <p>The random rotation matrix is derived from: Arvo, James (1992), "Fast random rotation
   * matrices", in David Kirk, Graphics Gems III, San Diego: Academic Press Professional, pp.
   * 117–120, ISBN 978-0-12-409671-4
   *
   * @param tr The translations to apply.
   * @return A Cartesian SymOp with a random rotation and translation.
   */
  public static SymOp randomSymOpFactory(double[] tr) {
    double[][] rot = new double[3][3];
    double PI2 = 2.0 * PI;
    double[] x = new double[3];
    x[0] = random();
    x[1] = random();
    x[2] = random();
    /* Rotation about the pole (Z).      */
    double theta = x[0] * PI2;
    /* For direction of pole deflection. */
    double phi = x[1] * PI2;
    /* For magnitude of pole deflection. */
    double z = x[2] * 2.0;
    /*
     Compute a vector V used for distributing points over the sphere via
     the reflection I - V Transpose(V). This formulation of V will
     guarantee that if x[1] and x[2] are uniformly distributed, the
     reflected points will be uniform on the sphere. Note that V has
     length sqrt(2) to eliminate the 2 in the Householder matrix.
    */
    double r = sqrt(z);
    double Vx = sin(phi) * r;
    double Vy = cos(phi) * r;
    double Vz = sqrt(2.0 - z);
    /*
     Compute the row vector S = Transpose(V) * R, where R is a simple
     rotation by theta about the z-axis. No need to compute Sz since it's
     just Vz.
    */
    double st = sin(theta);
    double ct = cos(theta);
    double Sx = Vx * ct - Vy * st;
    double Sy = Vx * st + Vy * ct;

    // Construct the rotation matrix ( V Transpose(V) - I ) R, which is equivalent to V S - R.
    rot[0][0] = Vx * Sx - ct;
    rot[0][1] = Vx * Sy - st;
    rot[0][2] = Vx * Vz;
    rot[1][0] = Vy * Sx + st;
    rot[1][1] = Vy * Sy - ct;
    rot[1][2] = Vy * Vz;
    rot[2][0] = Vz * Sx;
    rot[2][1] = Vz * Sy;
    // This equals Vz * Vz - 1.0
    rot[2][2] = 1.0 - z;

    return new SymOp(rot, tr);
  }

  /**
   * trtoString
   *
   * @param tr a double.
   * @return a {@link java.lang.String} object.
   */
  private static String trtoString(double tr) {
    if (tr == f12) {
      return "+1/2";
    }
    if (tr == f13) {
      return "+1/3";
    }
    if (tr == f23) {
      return "+2/3";
    }
    if (tr == f14) {
      return "+1/4";
    }
    if (tr == f34) {
      return "+3/4";
    }
    if (tr == f16) {
      return "+1/6";
    }
    if (tr == f56) {
      return "+5/6";
    }

    return "";
  }

  /** {@inheritDoc} */
  @Override
  public String toString() {
    StringBuilder sb = new StringBuilder(" Rotation operator:\n");
    sb.append(format(" [[%4.1f,%4.1f,%4.1f]\n  [%4.1f,%4.1f,%4.1f]\n  [%4.1f,%4.1f,%4.1f]]\n",
        rot[0][0], rot[0][1], rot[0][2],
        rot[1][0], rot[1][1], rot[1][2],
        rot[2][0], rot[2][1], rot[2][2]));
    sb.append(" Translation:\n");
    sb.append(format(" [%4.2f,%4.2f,%4.2f]", tr[0], tr[1], tr[2]));
    return sb.toString();
  }

  /**
   * Print the symmetry operator with double precision.
   *
   * @return String of rotation/translation with double precision.
   */
  public String toStringPrecise() {
    StringBuilder sb = new StringBuilder(" Rotation operator:\n");
    sb.append(format(
        " [[%18.16e,%18.16e,%18.16e]\n  [%18.16e,%18.16e,%18.16e]\n  [%18.16e,%18.16e,%18.16e]]\n",
        rot[0][0], rot[0][1], rot[0][2],
        rot[1][0], rot[1][1], rot[1][2],
        rot[2][0], rot[2][1], rot[2][2]));
    sb.append(" Translation:\n");
    sb.append(format(" [%18.16e,%18.16e,%18.16e]", tr[0], tr[1], tr[2]));
    return sb.toString();
  }

  /**
   * toXYZString
   *
   * @return a {@link java.lang.String} object.
   */
  public String toXYZString() {
    StringBuilder sb = new StringBuilder();
    for (int i = 0; i < 3; i++) {
      boolean s = false;
      if (rot[i][0] < 0.0) {
        sb.append("-X");
        s = true;
      } else if (rot[i][0] > 0.0) {
        sb.append("X");
        s = true;
      }
      if (rot[i][1] < 0.0) {
        sb.append("-Y");
        s = true;
      } else if (rot[i][1] > 0.0) {
        sb.append(s ? "+Y" : "Y");
        s = true;
      }
      if (rot[i][2] < 0.0) {
        sb.append("-Z");
      } else if (rot[i][2] > 0.0) {
        sb.append(s ? "+Z" : "Z");
      }

      if (tr[i] > 0.0) {
        sb.append(trtoString(tr[i]));
      }

      if (i < 2) {
        sb.append(", ");
      }
    }
    return sb.toString();
  }
}