UnitCell.cpp

#include "MantidGeometry/Crystal/UnitCell.h"
#include "MantidKernel/Matrix.h"
#include "MantidKernel/StringTokenizer.h"
#include "MantidKernel/System.h"
#include "MantidKernel/V3D.h"
#include <cfloat>
#include <iomanip>
#include <ios>
#include <stdexcept>

#include <boost/lexical_cast.hpp>

namespace Mantid {
namespace Geometry {
using Mantid::Kernel::DblMatrix;
using Mantid::Kernel::V3D;

/** Default constructor.
 \f$ a = b = c =  1 \mbox{\AA, } \alpha = \beta = \gamma = 90^\circ \f$ */
UnitCell::UnitCell()
    : da(6), ra(6), errorda(6), G(3, 3), Gstar(3, 3), B(3, 3), ModHKL(3, 3),
      errorModHKL(3, 3) {
  da[0] = da[1] = da[2] = 1.;
  da[3] = da[4] = da[5] = deg2rad * 90.0;
  errorda[0] = errorda[1] = errorda[2] = errorda[3] = errorda[4] = errorda[5] =
      0.0;
  MaxOrder = 0;
  CrossTerm = false;
  recalculate();
}

/** Constructor
 @param _a, _b, _c :: lattice parameters \f$ a, b, c \f$ \n
 with \f$\alpha = \beta = \gamma = 90^\circ \f$*/
UnitCell::UnitCell(double _a, double _b, double _c)
    : da(6), ra(6), errorda(6), G(3, 3), Gstar(3, 3), B(3, 3), ModHKL(3, 3),
      errorModHKL(3, 3) {
  da[0] = _a;
  da[1] = _b;
  da[2] = _c;
  // Angles are 90 degrees in radians ->Pi/2
  da[3] = da[4] = da[5] = 0.5 * M_PI;
  errorda[0] = errorda[1] = errorda[2] = errorda[3] = errorda[4] = errorda[5] =
      0.0;
  MaxOrder = 0;
  CrossTerm = false;
  recalculate();
}

/** Constructor
 @param _a, _b, _c, _alpha, _beta, _gamma :: lattice parameters\n
 @param angleunit :: units for angle, of type #AngleUnits. Default is degrees.
 */
UnitCell::UnitCell(double _a, double _b, double _c, double _alpha, double _beta,
                   double _gamma, const int angleunit)
    : da(6), ra(6), errorda(6), G(3, 3), Gstar(3, 3), B(3, 3), ModHKL(3, 3),
      errorModHKL(3, 3) {
  da[0] = _a;
  da[1] = _b;
  da[2] = _c;
  // Angle transformed in radians
  if (angleunit == angDegrees) {
    da[3] = deg2rad * _alpha;
    da[4] = deg2rad * _beta;
    da[5] = deg2rad * _gamma;
  } else {
    da[3] = _alpha;
    da[4] = _beta;
    da[5] = _gamma;
  }
  errorda[0] = errorda[1] = errorda[2] = errorda[3] = errorda[4] = errorda[5] =
      0.0;
  MaxOrder = 0;
  CrossTerm = false;
  recalculate();
}

/** Get lattice parameter
 @return a1 :: lattice parameter \f$ a \f$ (in \f$ \mbox{\AA} \f$ )
 @see a()*/
double UnitCell::a1() const { return da[0]; }

/** Get lattice parameter
 @return a2 :: lattice parameter \f$ b \f$ (in \f$ \mbox{\AA} \f$ )
 @see b()*/
double UnitCell::a2() const { return da[1]; }

/** Get lattice parameter
 @return a3 :: lattice parameter \f$ c \f$ (in \f$ \mbox{\AA} \f$ )
 @see c()*/

double UnitCell::a3() const { return da[2]; }
/** Get lattice parameter a1-a3 as function of index (0-2)
 @return a_n :: lattice parameter \f$ a,b or c \f$ (in \f$ \mbox{\AA} \f$ )
 */
double UnitCell::a(int nd) const {
  if (nd < 0 || nd > 2)
    throw(std::invalid_argument(
        "lattice parameter index can change from 0 to 2 "));
  return da[nd];
}

/** Get lattice parameter
 @return alpha1 :: lattice parameter \f$ \alpha \f$ (in radians)
 @see alpha()*/
double UnitCell::alpha1() const { return da[3]; }

/** Get lattice parameter
 @return alpha2 :: lattice parameter \f$ \beta \f$ (in radians)
 @see beta()*/
double UnitCell::alpha2() const { return da[4]; }

/** Get lattice parameter
 @return alpha3 :: lattice parameter \f$ \gamma \f$ (in radians)
 @see gamma()*/
double UnitCell::alpha3() const { return da[5]; }

/** Get lattice parameter
 @return a :: lattice parameter \f$ a \f$ (in \f$ \mbox{\AA} \f$ )
 @see a1()*/
double UnitCell::a() const { return da[0]; }

/** Get lattice parameter
 @return b :: lattice parameter \f$ b \f$ (in \f$ \mbox{\AA} \f$ )
 @see a2()*/
double UnitCell::b() const { return da[1]; }

/** Get lattice parameter
 @return c :: lattice parameter \f$ c \f$ (in \f$ \mbox{\AA} \f$ )
 @see a3()*/
double UnitCell::c() const { return da[2]; }

/** Get lattice parameter
 @return alpha :: lattice parameter \f$ \alpha \f$ (in degrees)
 @see alpha1()*/
double UnitCell::alpha() const { return da[3] * rad2deg; }

/** Get lattice parameter
 @return beta :: lattice parameter \f$ \beta \f$ (in degrees)
 @see alpha2()*/
double UnitCell::beta() const { return da[4] * rad2deg; }

/** Get lattice parameter
 @return gamma :: lattice parameter \f$ \gamma \f$ (in degrees)
 @see alpha3()*/
double UnitCell::gamma() const { return da[5] * rad2deg; }

/** Get reciprocal lattice parameter
 @return b1 :: lattice parameter \f$ a^{*} \f$ (in \f$ \mbox{\AA}^{-1} \f$ )
 @see astar()*/
double UnitCell::b1() const { return ra[0]; }

/** Get reciprocal lattice parameter
 @return b2 :: lattice parameter \f$ b^{*} \f$ (in \f$ \mbox{\AA}^{-1} \f$ )
 @see bstar()*/
double UnitCell::b2() const { return ra[1]; }

/** Get reciprocal lattice parameter
 @return b3 :: lattice parameter \f$ c^{*} \f$ (in \f$ \mbox{\AA}^{-1} \f$ )
 @see cstar()*/
double UnitCell::b3() const { return ra[2]; }

/** Get reciprocal lattice parameter
 @return beta1 :: lattice parameter \f$ \alpha^{*} \f$ (in radians)
 @see alphastar()*/
double UnitCell::beta1() const { return ra[3]; }

/** Get reciprocal lattice parameter
 @return beta2 :: lattice parameter \f$ \beta^{*} \f$ (in radians)
 @see betastar()*/
double UnitCell::beta2() const { return ra[4]; }

/** Get reciprocal lattice parameter
 @return beta3 :: lattice parameter \f$ \gamma^{*} \f$ (in radians)
 @see gammastar()*/
double UnitCell::beta3() const { return ra[5]; }

/** Get reciprocal lattice parameter
 @return astar :: lattice parameter \f$ a^{*} \f$ (in \f$ \mbox{\AA}^{-1} \f$ )
 @see b1()*/
double UnitCell::astar() const { return ra[0]; }

/** Get reciprocal lattice parameter
 @return bstar :: lattice parameter \f$ b^{*} \f$ (in \f$ \mbox{\AA}^{-1} \f$ )
 @see b2()*/
double UnitCell::bstar() const { return ra[1]; }

/** Get reciprocal lattice parameter
 @return cstar :: lattice parameter \f$ c^{*} \f$ (in \f$ \mbox{\AA}^{-1} \f$ )
 @see b3()*/
double UnitCell::cstar() const { return ra[2]; }

/** Get reciprocal lattice parameter
 @return alphastar :: lattice parameter \f$ \alpha^{*} \f$ (in degrees)
 @see beta1()*/
double UnitCell::alphastar() const { return ra[3] * rad2deg; }

/** Get reciprocal lattice parameter
 @return  betastar:: lattice parameter \f$ \beta^{*} \f$ (in degrees)
 @see beta2()*/
double UnitCell::betastar() const { return ra[4] * rad2deg; }

/** Get reciprocal lattice parameter
 @return  gammastar:: lattice parameter \f$ \gamma^{*} \f$ (in degrees)
 @see beta3()*/
double UnitCell::gammastar() const { return ra[5] * rad2deg; }

/** Get lattice parameter error
 @return errora :: errorlattice parameter \f$ a \f$ (in \f$ \mbox{\AA} \f$ )
 */
double UnitCell::errora() const { return errorda[0]; }

/** Get lattice parameter error
 @return errorb :: errorlattice parameter \f$ b \f$ (in \f$ \mbox{\AA} \f$ )
 */
double UnitCell::errorb() const { return errorda[1]; }

/** Get lattice parameter error
 @return errorc :: errorlattice parameter \f$ c \f$ (in \f$ \mbox{\AA} \f$ )
 */
double UnitCell::errorc() const { return errorda[2]; }

/** Get lattice parameter error
 @param angleunit :: units for angle, of type #AngleUnits . Default is degrees.
 @return erroralpha :: errorlattice parameter \f$ alpha \f$ (in degrees or
 radians )
 */
double UnitCell::erroralpha(const int angleunit) const {
  if (angleunit == angDegrees) {
    return errorda[3] * rad2deg;
  } else {
    return errorda[3];
  }
}

/** Get lattice parameter error
 @param angleunit :: units for angle, of type #AngleUnits . Default is degrees.
 @return erroralpha :: errorlattice parameter \f$ beta \f$ (in degrees or
 radians
 )
 */
double UnitCell::errorbeta(const int angleunit) const {
  if (angleunit == angDegrees) {
    return errorda[4] * rad2deg;
  } else {
    return errorda[4];
  }
}

/** Get lattice parameter error
 @param angleunit :: units for angle, of type #AngleUnits . Default is degrees.
 @return erroralpha :: errorlattice parameter \f$ gamma \f$ (in degrees or
 radians )
 */
double UnitCell::errorgamma(const int angleunit) const {
  if (angleunit == angDegrees) {
    return errorda[5] * rad2deg;
  } else {
    return errorda[5];
  }
}

/** Get lattice parameter error
 @return errorc :: errorlattice parameter \f$ volume \f$ (in \f$ \mbox{\AA} \f$
 )
 */
double UnitCell::errorvolume() const {
  // From latcon.py by Art Schultz
  double V = volume();
  double delta_V_alphaV = 0.0;
  if (erroralpha() > 0.0) {
    double alpha1 = alpha() - 0.5 * erroralpha();
    double Va1 = UnitCell(a(), b(), c(), alpha1, beta(), gamma()).volume();
    double alpha2 = alpha() + 0.5 * erroralpha();
    double Va2 = UnitCell(a(), b(), c(), alpha2, beta(), gamma()).volume();
    delta_V_alphaV = (Va2 - Va1) / V;
  }

  double delta_V_betaV = 0.0;
  if (errorbeta() > 0.0) {
    double beta1 = beta() - 0.5 * errorbeta();
    double Va1 = UnitCell(a(), b(), c(), alpha(), beta1, gamma()).volume();
    double beta2 = beta() + 0.5 * errorbeta();
    double Va2 = UnitCell(a(), b(), c(), alpha(), beta2, gamma()).volume();
    delta_V_betaV = (Va2 - Va1) / V;
  }

  double delta_V_gammaV = 0.0;
  if (errorgamma() > 0.0) {
    double gamma1 = gamma() - 0.5 * errorgamma();
    double Va1 = UnitCell(a(), b(), c(), alpha(), beta(), gamma1).volume();
    double gamma2 = gamma() + 0.5 * errorgamma();
    double Va2 = UnitCell(a(), b(), c(), alpha(), beta(), gamma2).volume();
    delta_V_gammaV = (Va2 - Va1) / V;
  }

  return V * sqrt(std::pow(errora() / a(), 2) + std::pow(errorb() / b(), 2) +
                  std::pow(errorc() / c(), 2) + std::pow(delta_V_alphaV, 2) +
                  std::pow(delta_V_betaV, 2) + std::pow(delta_V_gammaV, 2));
}

/** Set lattice parameters
 @param _a, _b, _c, _alpha, _beta, _gamma :: lattice parameters\n
 @param angleunit :: units for angle, of type #AngleUnits . Default is degrees.
 */

void UnitCell::set(double _a, double _b, double _c, double _alpha, double _beta,
                   double _gamma, const int angleunit) {
  da[0] = _a;
  da[1] = _b;
  da[2] = _c;
  if (angleunit == angDegrees) {
    da[3] = deg2rad * _alpha;
    da[4] = deg2rad * _beta;
    da[5] = deg2rad * _gamma;
  } else {
    da[3] = _alpha;
    da[4] = _beta;
    da[5] = _gamma;
  }
  recalculate();
}

/** Set lattice parameter errors
 @param _aerr, _berr, _cerr, _alphaerr, _betaerr, _gammaerr :: lattice
 parameter errors\n
 @param angleunit :: units for angle, of type #AngleUnits . Default is degrees.
 */

void UnitCell::setError(double _aerr, double _berr, double _cerr,
                        double _alphaerr, double _betaerr, double _gammaerr,
                        const int angleunit) {
  errorda[0] = _aerr;
  errorda[1] = _berr;
  errorda[2] = _cerr;
  if (angleunit == angDegrees) {
    errorda[3] = deg2rad * _alphaerr;
    errorda[4] = deg2rad * _betaerr;
    errorda[5] = deg2rad * _gammaerr;
  } else {
    errorda[3] = _alphaerr;
    errorda[4] = _betaerr;
    errorda[5] = _gammaerr;
  }
}

void UnitCell::setModHKL(double _dh1, double _dk1, double _dl1, double _dh2,
                         double _dk2, double _dl2, double _dh3, double _dk3,
                         double _dl3) {
  ModHKL[0][0] = _dh1;
  ModHKL[1][0] = _dk1;
  ModHKL[2][0] = _dl1;
  ModHKL[0][1] = _dh2;
  ModHKL[1][1] = _dk2;
  ModHKL[2][1] = _dl2;
  ModHKL[0][2] = _dh3;
  ModHKL[1][2] = _dk3;
  ModHKL[2][2] = _dl3;
}

void UnitCell::setModHKL(const DblMatrix &newModHKL) { ModHKL = newModHKL; }

void UnitCell::setErrorModHKL(const DblMatrix &newErrorModHKL) {
  errorModHKL = newErrorModHKL;
}

void UnitCell::setErrorModHKL(double _dh1err, double _dk1err, double _dl1err,
                              double _dh2err, double _dk2err, double _dl2err,
                              double _dh3err, double _dk3err, double _dl3err) {
  errorModHKL[0][0] = _dh1err;
  errorModHKL[1][0] = _dk1err;
  errorModHKL[2][0] = _dl1err;
  errorModHKL[0][1] = _dh2err;
  errorModHKL[1][1] = _dk2err;
  errorModHKL[2][1] = _dl2err;
  errorModHKL[0][2] = _dh3err;
  errorModHKL[1][2] = _dk3err;
  errorModHKL[2][2] = _dl3err;
}

void UnitCell::setModVec1(double _dh1, double _dk1, double _dl1) {
  ModHKL[0][0] = _dh1;
  ModHKL[1][0] = _dk1;
  ModHKL[2][0] = _dl1;
}

void UnitCell::setModVec2(double _dh2, double _dk2, double _dl2) {
  ModHKL[0][1] = _dh2;
  ModHKL[1][1] = _dk2;
  ModHKL[2][1] = _dl2;
}

void UnitCell::setModVec3(double _dh3, double _dk3, double _dl3) {
  ModHKL[0][2] = _dh3;
  ModHKL[1][2] = _dk3;
  ModHKL[2][2] = _dl3;
}

void UnitCell::setModVec1(const V3D &newModVec) {
  ModHKL[0][0] = newModVec[0];
  ModHKL[1][0] = newModVec[1];
  ModHKL[2][0] = newModVec[2];
}

void UnitCell::setModVec2(const V3D &newModVec) {
  ModHKL[0][1] = newModVec[0];
  ModHKL[1][1] = newModVec[1];
  ModHKL[2][1] = newModVec[2];
}

void UnitCell::setModVec3(const V3D &newModVec) {
  ModHKL[0][2] = newModVec[0];
  ModHKL[1][2] = newModVec[1];
  ModHKL[2][2] = newModVec[2];
}

void UnitCell::setModerr(int i, double _dherr, double _dkerr, double _dlerr) {
  errorModHKL[0][i] = _dherr;
  errorModHKL[1][i] = _dkerr;
  errorModHKL[2][i] = _dlerr;
}

void UnitCell::setModerr1(double _dh1err, double _dk1err, double _dl1err) {
  errorModHKL[0][0] = _dh1err;
  errorModHKL[1][0] = _dk1err;
  errorModHKL[2][0] = _dl1err;
}

void UnitCell::setModerr2(double _dh2err, double _dk2err, double _dl2err) {
  errorModHKL[0][1] = _dh2err;
  errorModHKL[1][1] = _dk2err;
  errorModHKL[2][1] = _dl2err;
}

void UnitCell::setModerr3(double _dh3err, double _dk3err, double _dl3err) {
  errorModHKL[0][2] = _dh3err;
  errorModHKL[1][2] = _dk3err;
  errorModHKL[2][2] = _dl3err;
}
void UnitCell::setMaxOrder(int MaxO) { MaxOrder = MaxO; }
void UnitCell::setCrossTerm(bool CT) { CrossTerm = CT; }

const Kernel::V3D UnitCell::getModVec(int j) const {
  return V3D(getdh(j), getdk(j), getdl(j));
}

const Kernel::V3D UnitCell::getVecErr(int j) const {
  return V3D(getdherr(j), getdkerr(j), getdlerr(j));
}

/*
const Kernel::V3D &UnitCell::getModVec2() const
{
    return V3D (getdh2(),getdk2(),getdl2());
}

const Kernel::V3D &UnitCell::getModVec3() const
{
    return V3D (getdh3(),getdk3(),getdl3());
}
*/

const Kernel::DblMatrix &UnitCell::getModHKL() const { return ModHKL; }

double UnitCell::getdh(int j) const { return ModHKL[0][j - 1]; }

double UnitCell::getdk(int j) const { return ModHKL[1][j - 1]; }

double UnitCell::getdl(int j) const { return ModHKL[2][j - 1]; }

double UnitCell::getdherr(int j) const { return errorModHKL[0][j - 1]; }

double UnitCell::getdkerr(int j) const { return errorModHKL[1][j - 1]; }

double UnitCell::getdlerr(int j) const { return errorModHKL[2][j - 1]; }

int UnitCell::getMaxOrder() const { return MaxOrder; }

bool UnitCell::getCrossTerm() const { return CrossTerm; }

/*
int UnitCell::getModDim()
{
    if (getdh(1)==0 && getdk(1)==0 && getdl(1)==0)
}
*/

/*
double UnitCell::getdh2() const { return ModHKL[0][1]; }

double UnitCell::getdk2() const { return ModHKL[1][1]; }

double UnitCell::getdl2() const { return ModHKL[2][1]; }

double UnitCell::getdh3() const { return ModHKL[0][2]; }

double UnitCell::getdk3() const { return ModHKL[1][2]; }

double UnitCell::getdl3() const { return ModHKL[2][2]; }
*/

/** Set lattice parameter
 @param _a :: lattice parameter \f$ a \f$ (in \f$ \mbox{\AA} \f$ )*/
void UnitCell::seta(double _a) {
  da[0] = _a;
  recalculate();
}
/** Set lattice parameter error
 @param _aerr :: lattice parameter \f$ a \f$ error (in \f$ \mbox{\AA} \f$ )*/
void UnitCell::setErrora(double _aerr) { errorda[0] = _aerr; }

/** Set lattice parameter
 @param _b :: lattice parameter \f$ b \f$ (in \f$ \mbox{\AA} \f$ )*/
void UnitCell::setb(double _b) {
  da[1] = _b;
  recalculate();
}
/** Set lattice parameter error
 @param _berr :: lattice parameter \f$ b \f$ error (in \f$ \mbox{\AA} \f$ )*/
void UnitCell::setErrorb(double _berr) { errorda[1] = _berr; }
/** Set lattice parameter
 @param _c :: lattice parameter \f$ c \f$ (in \f$ \mbox{\AA} \f$ )*/
void UnitCell::setc(double _c) {
  da[2] = _c;
  recalculate();
}
/** Set lattice parameter error
 @param _cerr :: lattice parameter \f$ c \f$ error (in \f$ \mbox{\AA} \f$ )*/
void UnitCell::setErrorc(double _cerr) { errorda[2] = _cerr; }
/** Set lattice parameter
 @param _alpha :: lattice parameter \f$ \alpha \f$
 @param angleunit :: units for angle, of type #AngleUnits. Default is degrees.
 */
void UnitCell::setalpha(double _alpha, const int angleunit) {
  if (angleunit == angDegrees)
    da[3] = deg2rad * _alpha;
  else
    da[3] = _alpha;
  recalculate();
}
/** Set lattice parameter error
 @param _alphaerr :: lattice parameter \f$ \alpha \f$ error
 @param angleunit :: units for angle, of type #AngleUnits. Default is degrees.
 */
void UnitCell::setErroralpha(double _alphaerr, const int angleunit) {
  if (angleunit == angDegrees)
    errorda[3] = deg2rad * _alphaerr;
  else
    errorda[3] = _alphaerr;
}
/** Set lattice parameter
 @param _beta :: lattice parameter \f$ \beta \f$
 @param angleunit :: units for angle, of type #AngleUnits. Default is degrees.
 */
void UnitCell::setbeta(double _beta, const int angleunit) {
  if (angleunit == angDegrees)
    da[4] = deg2rad * _beta;
  else
    da[4] = _beta;
  recalculate();
}

/** Set lattice parameter error
 @param _betaerr :: lattice parameter \f$ \beta \f$ error
 @param angleunit :: units for angle, of type #AngleUnits. Default is degrees.
 */
void UnitCell::setErrorbeta(double _betaerr, const int angleunit) {
  if (angleunit == angDegrees)
    errorda[4] = deg2rad * _betaerr;
  else
    errorda[4] = _betaerr;
}

/** Set lattice parameter
 @param _gamma :: lattice parameter \f$ \gamma \f$
 @param angleunit :: units for angle, of type #AngleUnits. Default is degrees.
 */
void UnitCell::setgamma(double _gamma, const int angleunit) {
  if (angleunit == angDegrees)
    da[5] = deg2rad * _gamma;
  else
    da[5] = _gamma;
  recalculate();
}

/** Set lattice parameter error
 @param _gammaerr :: lattice parameter \f$ \gamma \f$ error
 @param angleunit :: units for angle, of type #AngleUnits. Default is degrees.
 */
void UnitCell::setErrorgamma(double _gammaerr, const int angleunit) {
  if (angleunit == angDegrees)
    errorda[5] = deg2rad * _gammaerr;
  else
    errorda[5] = _gammaerr;
}

/// Return d-spacing (\f$ \mbox{ \AA } \f$) for a given h,k,l coordinate
double UnitCell::d(double h, double k, double l) const {
  return 1.0 / dstar(V3D(h, k, l));
}

/// Return d-spacing (\f$ \mbox{ \AA } \f$) for a given h,k,l coordinate
double UnitCell::d(const V3D &hkl) const { return 1.0 / dstar(hkl); }

/// Return d*=1/d (\f$ \mbox{ \AA }^{-1} \f$) for a given h,k,l coordinate
double UnitCell::dstar(double h, double k, double l) const {
  return dstar(V3D(h, k, l)); // create a V3D vector h,k,l
}

/// Return d*=1/d (\f$ \mbox{ \AA }^{-1} \f$) for a given h,k,l coordinate
double UnitCell::dstar(const V3D &hkl) const {
  V3D Q = B * hkl; // transform into $AA^-1$
  return Q.norm();
}

/// Calculate the angle in degrees or radians between two reciprocal vectors
/// (h1,k1,l1) and (h2,k2,l2)
double UnitCell::recAngle(double h1, double k1, double l1, double h2, double k2,
                          double l2, const int angleunit) const {
  V3D Q1(h1, k1, l1), Q2(h2, k2, l2);
  double E, ang;
  Q1 = Gstar * Q1;
  E = Q1.scalar_prod(Q2);
  double temp = E / dstar(h1, k1, l1) / dstar(h2, k2, l2);
  if (temp > 1)
    ang = 0.;
  else if (temp < -1)
    ang = M_PI;
  else
    ang = acos(temp);
  if (angleunit == angDegrees)
    return rad2deg * ang;
  else
    return ang;
}

/// Volume of the direct unit-cell
double UnitCell::volume() const {
  double volume = G.determinant();
  return sqrt(volume);
}

/// Volume of the reciprocal lattice
double UnitCell::recVolume() const {
  double recvolume = Gstar.determinant();
  return sqrt(recvolume);
}

/// Get the metric tensor
/// @return G :: metric tensor
const Kernel::DblMatrix &UnitCell::getG() const { return G; }

/// Get the reciprocal metric tensor
/// @return Gstar :: metric tensor of the reciprocal lattice
const Kernel::DblMatrix &UnitCell::getGstar() const { return Gstar; }

/// Get the B-matrix
/// @return B :: B matrix in Busing-Levy convention
const Kernel::DblMatrix &UnitCell::getB() const { return B; }

/// Get the inverse of the B-matrix
/// @return Binv :: inverse of the B matrix in Busing-Levy convention
const Kernel::DblMatrix &UnitCell::getBinv() const { return Binv; }

/// Private function, called at initialization or whenever lattice parameters
/// are changed
void UnitCell::recalculate() {
  if ((da[3] > da[4] + da[5]) || (da[4] > da[3] + da[5]) ||
      (da[5] > da[4] + da[3])) {
    throw std::invalid_argument("Invalid angles");
  }
  calculateG();
  calculateGstar();
  calculateReciprocalLattice();
  calculateB();
}

/// Private function to calculate #G matrix

/*
        void UnitCell::calculateModVec()
        {
            ModVec1(ModHKL[0][0],ModHKL[1][0],ModHKL[2][0]);
            ModVec1(ModHKL[0][1],ModHKL[1][1],ModHKL[2][1]);
            ModVec1(ModHKL[0][2],ModHKL[1][2],ModHKL[2][2]);
        }
 */

void UnitCell::calculateG() {
  G[0][0] = da[0] * da[0];
  G[1][1] = da[1] * da[1];
  G[2][2] = da[2] * da[2];
  G[0][1] = da[0] * da[1] * cos(da[5]);
  G[0][2] = da[0] * da[2] * cos(da[4]);
  G[1][2] = da[1] * da[2] * cos(da[3]);
  G[1][0] = G[0][1];
  G[2][0] = G[0][2];
  G[2][1] = G[1][2];
}

/// Private function to calculate #Gstar matrix
void UnitCell::calculateGstar() {
  // Reciprocal metrix tensor is simply the inverse of the direct one
  double det = G.determinant();
  if (det == 0) {
    throw std::range_error("UnitCell not properly initialized");
  }
  Gstar = G;
  if (Gstar.Invert() == 0) {
    throw std::range_error("UnitCell not properly initialized");
  }
}

/// Private function to calculate reciprocal lattice parameters
void UnitCell::calculateReciprocalLattice() {
  ra[0] = sqrt(Gstar[0][0]);                 // a*
  ra[1] = sqrt(Gstar[1][1]);                 // b*
  ra[2] = sqrt(Gstar[2][2]);                 // c*
  ra[3] = acos(Gstar[1][2] / ra[1] / ra[2]); // alpha*
  ra[4] = acos(Gstar[0][2] / ra[0] / ra[2]); // beta*
  ra[5] = acos(Gstar[0][1] / ra[0] / ra[1]); // gamma*
}

/// Private function to calculate #B matrix
void UnitCell::calculateB() {
  // B matrix using a right handed coordinate system with b1 along x and y in
  // the (b1,b2) plane.
  // This is the convention in Busing and Levy.
  // | b1 b2cos(beta3)      b3cos(beta2)        |
  // | 0  b2sin(beta3) -b3sin(beta2)cos(alpha1) |
  // | 0       0                  1/a3          |
  B[0][0] = ra[0];
  B[0][1] = ra[1] * cos(ra[5]);
  B[0][2] = ra[2] * cos(ra[4]);
  B[1][0] = 0.;
  B[1][1] = ra[1] * sin(ra[5]);
  B[1][2] = -ra[2] * sin(ra[4]) * cos(da[3]);
  B[2][0] = 0.;
  B[2][1] = 0.;
  B[2][2] = 1. / da[2];

  /// Now let's cache the inverse B
  Binv = B;
  Binv.Invert();
}

/// Recalculate lattice from reciprocal metric tensor (Gstar=transpose(UB)*UB)
void UnitCell::recalculateFromGstar(const DblMatrix &NewGstar) {
  if (NewGstar.numRows() != 3 || NewGstar.numCols() != 3) {
    std::ostringstream msg;
    msg << "UnitCell::recalculateFromGstar - Expected a 3x3 matrix but was "
           "given a "
        << NewGstar.numRows() << "x" << NewGstar.numCols();
    throw std::invalid_argument(msg.str());
  }

  if (NewGstar[0][0] * NewGstar[1][1] * NewGstar[2][2] <= 0.)
    throw std::invalid_argument("NewGstar");
  Gstar = NewGstar;
  calculateReciprocalLattice();
  G = Gstar;
  G.Invert();
  da[0] = sqrt(G[0][0]);                 // a
  da[1] = sqrt(G[1][1]);                 // b
  da[2] = sqrt(G[2][2]);                 // c
  da[3] = acos(G[1][2] / da[1] / da[2]); // alpha
  da[4] = acos(G[0][2] / da[0] / da[2]); // beta
  da[5] = acos(G[0][1] / da[0] / da[1]); // gamma
  calculateB();
}

std::ostream &operator<<(std::ostream &out, const UnitCell &unitCell) {
  // always show the lattice constants
  out << "Lattice Parameters:" << std::fixed << std::setprecision(6)
      << std::setw(12) << unitCell.a() << std::fixed << std::setprecision(6)
      << std::setw(12) << unitCell.b() << std::fixed << std::setprecision(6)
      << std::setw(12) << unitCell.c() << std::fixed << std::setprecision(6)
      << std::setw(12) << unitCell.alpha() << std::fixed << std::setprecision(6)
      << std::setw(12) << unitCell.beta() << std::fixed << std::setprecision(6)
      << std::setw(12) << unitCell.gamma() << std::fixed << std::setprecision(6)
      << " " << std::setw(12) << unitCell.volume();

  // write out the uncertainty if there is a positive one somewhere
  if ((unitCell.errora() > 0) || (unitCell.errorb() > 0) ||
      (unitCell.errorc() > 0) || (unitCell.erroralpha() > 0) ||
      (unitCell.errorbeta() > 0) || (unitCell.errorgamma() > 0))
    out << "\nParameter Errors  :" << std::fixed << std::setprecision(6)
        << std::setw(12) << unitCell.errora() << std::fixed
        << std::setprecision(6) << std::setw(12) << unitCell.errorb()
        << std::fixed << std::setprecision(6) << std::setw(12)
        << unitCell.errorc() << std::fixed << std::setprecision(6)
        << std::setw(12) << unitCell.erroralpha() << std::fixed
        << std::setprecision(6) << std::setw(12) << unitCell.errorbeta()
        << std::fixed << std::setprecision(6) << std::setw(12)
        << unitCell.errorgamma() << std::fixed << std::setprecision(6)
        << std::setw(12) << unitCell.errorvolume();

  return out;
}

std::string unitCellToStr(const UnitCell &unitCell) {
  std::ostringstream stream;
  stream << std::setprecision(9);

  stream << unitCell.a() << " " << unitCell.b() << " " << unitCell.c() << " "
         << unitCell.alpha() << " " << unitCell.beta() << " "
         << unitCell.gamma();

  return stream.str();
}

UnitCell strToUnitCell(const std::string &unitCellString) {

  Mantid::Kernel::StringTokenizer cellTokens(
      unitCellString, " ", Mantid::Kernel::StringTokenizer::TOK_IGNORE_EMPTY);

  std::vector<double> components;

  for (const auto &token : cellTokens) {
    components.push_back(boost::lexical_cast<double>(token));
  }

  switch (components.size()) {
  case 3:
    return UnitCell(components[0], components[1], components[2]);
  case 6:
    return UnitCell(components[0], components[1], components[2], components[3],
                    components[4], components[5]);
  default:
    throw std::runtime_error("Failed to parse unit cell input string: " +
                             unitCellString);
  }
}

} // namespace Geometry
} // namespace Mantid