ProSHADE  0.7.5.3 (FEB 2021)
Protein Shape Detection
ProSHADE_internal_maths Namespace Reference

This namespace contains the internal functions for common mathematical operations. More...

Classes

class  BicubicInterpolator
 

Functions

void complexMultiplication (proshade_double *r1, proshade_double *i1, proshade_double *r2, proshade_double *i2, proshade_double *retReal, proshade_double *retImag)
 Function to multiply two complex numbers. More...
 
void complexMultiplicationConjug (proshade_double *r1, proshade_double *i1, proshade_double *r2, proshade_double *i2, proshade_double *retReal, proshade_double *retImag)
 Function to multiply two complex numbers by using the second number's conjugate. More...
 
proshade_double complexMultiplicationRealOnly (proshade_double *r1, proshade_double *i1, proshade_double *r2, proshade_double *i2)
 Function to multiply two complex numbers and return the real part only. More...
 
proshade_double complexMultiplicationConjugRealOnly (proshade_double *r1, proshade_double *i1, proshade_double *r2, proshade_double *i2)
 Function to conjuggate multiply two complex numbers and return the real part only. More...
 
void vectorMeanAndSD (std::vector< proshade_double > *vec, proshade_double *&ret)
 Function to get vector mean and standard deviation. More...
 
void vectorMedianAndIQR (std::vector< proshade_double > *vec, proshade_double *&ret)
 Function to get vector median and inter-quartile range. More...
 
void arrayMedianAndIQR (proshade_double *vec, proshade_unsign vecSize, proshade_double *&ret)
 Function to get array median and inter-quartile range. More...
 
proshade_double pearsonCorrCoeff (proshade_double *valSet1, proshade_double *valSet2, proshade_unsign length)
 Function for computing the Pearson's correlation coefficient. More...
 
void getLegendreAbscAndWeights (proshade_unsign order, proshade_double *abscissas, proshade_double *weights, proshade_unsign taylorSeriesCap)
 Function to prepare abscissas and weights for Gauss-Legendre integration. More...
 
void getGLPolyAtZero (proshade_unsign order, proshade_double *polyValue, proshade_double *deriValue)
 This function obtains the Legendre polynomial values and its derivative at zero for any positive integer order polynomial. More...
 
void getGLFirstEvenRoot (proshade_double polyAtZero, proshade_unsign order, proshade_double *abscAtZero, proshade_double *weighAtZero, proshade_unsign taylorSeriesCap)
 This function finds the first root for Legendre polynomials of odd order. More...
 
proshade_double evaluateGLSeries (proshade_double *series, proshade_double target, proshade_unsign terms)
 This function evaluates the Taylor expansion. More...
 
proshade_double advanceGLPolyValue (proshade_double from, proshade_double to, proshade_double valAtFrom, proshade_unsign noSteps, proshade_unsign taylorSeriesCap)
 This function finds the next value of the polynomial. More...
 
void completeLegendreSeries (proshade_unsign order, proshade_double *abscissa, proshade_double *weights, proshade_unsign taylorSeriesCap)
 This function completes the Legendre polynomial series assuming you have obtained the first values. More...
 
proshade_double gaussLegendreIntegrationReal (proshade_double *vals, proshade_unsign valsSize, proshade_unsign order, proshade_double *abscissas, proshade_double *weights, proshade_double integralOverRange, proshade_double maxSphereDists)
 Function to compute real part of the Gauss-Legendre integration over spherical harmonic values in different shells. More...
 
void gaussLegendreIntegration (proshade_complex *vals, proshade_unsign valsSize, proshade_unsign order, proshade_double *abscissas, proshade_double *weights, proshade_double integralOverRange, proshade_double maxSphereDists, proshade_double *retReal, proshade_double *retImag)
 Function to compute the complete complex Gauss-Legendre integration over spherical harmonic values in different shells. More...
 
void complexMatrixSVDSigmasOnly (proshade_complex **mat, int dim, double *&singularValues)
 Function to compute the complete complex matrix SVD and return only the sigmas. More...
 
void complexMatrixSVDUandVOnly (proshade_double *mat, int dim, proshade_double *uAndV, bool fail=true)
 Function to compute the real matrix SVD and return the U and V matrices. More...
 
void getEulerZXZFromSOFTPosition (proshade_signed band, proshade_signed x, proshade_signed y, proshade_signed z, proshade_double *eulerAlpha, proshade_double *eulerBeta, proshade_double *eulerGamma)
 Function to find Euler angles (ZXZ convention) from index position in the inverse SOFT map. More...
 
void getSOFTPositionFromEulerZXZ (proshade_signed band, proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma, proshade_double *x, proshade_double *y, proshade_double *z)
 Function to find the index position in the inverse SOFT map from given Euler angles (ZXZ convention). More...
 
void getRotationMatrixFromEulerZXZAngles (proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma, proshade_double *matrix)
 Function to find the rotation matrix from Euler angles (ZXZ convention). More...
 
void getAxisAngleFromRotationMatrix (proshade_double *rotMat, proshade_double *x, proshade_double *y, proshade_double *z, proshade_double *ang)
 This function converts rotation matrix to the axis-angle representation. More...
 
void getAxisAngleFromRotationMatrix (std::vector< proshade_double > *rotMat, proshade_double *x, proshade_double *y, proshade_double *z, proshade_double *ang)
 This function converts rotation matrix to the axis-angle representation. More...
 
void getRotationMatrixFromAngleAxis (proshade_double *rotMat, proshade_double x, proshade_double y, proshade_double z, proshade_double ang)
 This function converts the axis-angle representation to the rotation matrix representation. More...
 
void getEulerZXZFromRotMatrix (proshade_double *rotMat, proshade_double *eA, proshade_double *eB, proshade_double *eG)
 This function converts rotation matrix to the Euler ZXZ angles representation. More...
 
void getEulerZXZFromAngleAxis (proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axAng, proshade_double *eA, proshade_double *eB, proshade_double *eG, proshade_unsign angDim)
 This function converts angle-axis representation to the Euler ZXZ angles representation. More...
 
void getEulerZXZFromAngleAxisFullSearch (proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axAng, proshade_double *eA, proshade_double *eB, proshade_double *eG, proshade_signed angDim)
 This function converts angle-axis representation to the Euler ZXZ angles representation using full search. More...
 
void multiplyTwoSquareMatrices (proshade_double *A, proshade_double *B, proshade_double *res, proshade_unsign dim=3)
 Function to compute matrix multiplication. More...
 
std::vector< proshade_signed > primeFactorsDecomp (proshade_signed number)
 Function to find prime factors of an integer. More...
 
proshade_double normalDistributionValue (proshade_double mean, proshade_double standardDev, proshade_double value)
 Function to the heiht of normal distribution given by mean and standard deviation for a given value. More...
 
proshade_double computeDotProduct (proshade_double *x1, proshade_double *y1, proshade_double *z1, proshade_double *x2, proshade_double *y2, proshade_double *z2)
 Simple 3D vector dot product computation. More...
 
proshade_double computeDotProduct (proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2)
 Simple 3D vector dot product computation. More...
 
proshade_double * computeCrossProduct (proshade_double *x1, proshade_double *y1, proshade_double *z1, proshade_double *x2, proshade_double *y2, proshade_double *z2)
 Simple 3D vector cross product computation. More...
 
proshade_double * compute3x3MatrixMultiplication (proshade_double *mat1, proshade_double *mat2)
 Function for computing a 3x3 matrix multiplication. More...
 
proshade_double * compute3x3MatrixVectorMultiplication (proshade_double *mat, proshade_double x, proshade_double y, proshade_double z)
 Function for computing a 3x3 matrix to 3x1 vector multiplication. More...
 
proshade_double * compute3x3MatrixInverse (proshade_double *mat)
 Function for computing a 3x3 matrix inverse. More...
 
void transpose3x3MatrixInPlace (proshade_double *mat)
 Transposes 3x3 matrix in place. More...
 
proshade_double * findRotMatMatchingVectors (proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2)
 Computation of rotation matrix rotating one vector onto the other. More...
 
std::vector< proshade_double > findVectorFromTwoVAndTwoD (proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2, proshade_double dot1, proshade_double dot2)
 Function for finding a vector which would have a given two dot products to two other vectors. More...
 
std::vector< proshade_double > findVectorFromThreeVAndThreeD (proshade_double x1, proshade_double y1, proshade_double z1, proshade_double x2, proshade_double y2, proshade_double z2, proshade_double x3, proshade_double y3, proshade_double z3, proshade_double dot1, proshade_double dot2, proshade_double dot3)
 Function for finding a vector which would have a given three dot products to three other vectors. More...
 
std::vector< proshade_double > multiplyGroupElementMatrices (std::vector< proshade_double > *el1, std::vector< proshade_double > *el2)
 This function computes matrix multiplication using the ProSHADE group element matrix format as input and output. More...
 
bool rotationMatrixSimilarity (std::vector< proshade_double > *mat1, std::vector< proshade_double > *mat2, proshade_double tolerance=0.1)
 This function compares the distance between two rotation matrices and decides if they are similar using tolerance. More...
 
bool vectorOrientationSimilarity (proshade_double a1, proshade_double a2, proshade_double a3, proshade_double b1, proshade_double b2, proshade_double b3, proshade_double tolerance=0.1)
 This function compares two vectors using cosine distance and decides if they are similar using tolerance. More...
 
bool vectorOrientationSimilaritySameDirection (proshade_double a1, proshade_double a2, proshade_double a3, proshade_double b1, proshade_double b2, proshade_double b3, proshade_double tolerance=0.1)
 This function compares two vectors using cosine distance and decides if they are similar using tolerance. More...
 
void optimiseAxisBiCubicInterpolation (proshade_double *bestLattitude, proshade_double *bestLongitude, proshade_double *bestSum, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun, proshade_double step=0.05)
 This function provides axis optimisation given starting lattitude and longitude indices. More...
 
void prepareBiCubicInterpolatorsMinusMinus (proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
 This function prepares the interpolation objects for the bi-cubic interpolation. More...
 
void prepareBiCubicInterpolatorsMinusPlus (proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
 This function prepares the interpolation objects for the bi-cubic interpolation. More...
 
void prepareBiCubicInterpolatorsPlusMinus (proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
 This function prepares the interpolation objects for the bi-cubic interpolation. More...
 
void prepareBiCubicInterpolatorsPlusPlus (proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
 This function prepares the interpolation objects for the bi-cubic interpolation. More...
 
bool isAxisUnique (std::vector< proshade_double * > *CSymList, proshade_double *axis, proshade_double tolerance=0.1, bool improve=false)
 This function checks if new axis is unique, or already detected. More...
 
bool isAxisUnique (std::vector< proshade_double * > *CSymList, proshade_double X, proshade_double Y, proshade_double Z, proshade_double fold, proshade_double tolerance)
 This function checks if new axis is unique, or already detected. More...
 
std::vector< proshade_unsign > findAllPrimes (proshade_unsign upTo)
 This function finds all prime numbers up to the supplied limit. More...
 

Detailed Description

This namespace contains the internal functions for common mathematical operations.

The ProSHADE_internal_maths namespace contains a set of common mathematical operations used in many places by ProSHADE. These typically include complex number operations and angle conversions.

Function Documentation

◆ advanceGLPolyValue()

proshade_double ProSHADE_internal_maths::advanceGLPolyValue ( proshade_double  from,
proshade_double  to,
proshade_double  valAtFrom,
proshade_unsign  noSteps,
proshade_unsign  taylorSeriesCap 
)

This function finds the next value of the polynomial.

Given the previous value of the polynomial, the distance to proceed and the number of steps to take, this function finds the next value of the polynomial using the Taylor series.

Parameters
[in]fromCurrent polynomial position.
[in]toPolynomial position to move to.
[in]valAtFromThe current value of the polynomial at the <from> position.
[in]noStepsNumber of steps in which to reach the <to> position.
[in]taylorSeriesCapThe limit on the Taylor series.
[out]XThe polynomial value at the <to> position.

Definition at line 477 of file ProSHADE_maths.cpp.

478 {
479  //================================================ Initialise variables
480  proshade_double hlpVal = 0.0;
481  proshade_double stepSize = 0.0;
482  proshade_double valChange = 0.0;
483  proshade_double valSecChange = 0.0;
484  proshade_double squareSteps = 0.0;
485  proshade_double curVal = 0.0;
486 
487  //================================================ Set initial values
488  stepSize = ( to - from ) / static_cast<proshade_double> ( taylorSeriesCap );
489  squareSteps = sqrt ( static_cast<proshade_double> ( noSteps * ( noSteps + 1 ) ) );
490  curVal = from;
491 
492  //================================================ Go through the series and iteratively improve the estimate
493  for ( proshade_unsign iter = 0; iter < taylorSeriesCap; iter++ )
494  {
495  hlpVal = ( 1.0 - valAtFrom ) * ( 1.0 + valAtFrom );
496  valChange = - stepSize * hlpVal / ( squareSteps * sqrt ( hlpVal ) - 0.5 * valAtFrom * sin ( 2.0 * curVal ) );
497  valAtFrom = valAtFrom + valChange;
498 
499  curVal = curVal + stepSize;
500 
501  hlpVal = ( 1.0 - valAtFrom ) * ( 1.0 + valAtFrom );
502  valSecChange = - stepSize * hlpVal / ( squareSteps * sqrt ( hlpVal ) - 0.5 * valAtFrom * sin ( 2.0 * curVal ) );
503  valAtFrom = valAtFrom + 0.5 * ( valSecChange - valChange );
504  }
505 
506  //================================================ Done
507  return valAtFrom;
508 
509 }

◆ arrayMedianAndIQR()

void ProSHADE_internal_maths::arrayMedianAndIQR ( proshade_double *  vec,
proshade_unsign  vecSize,
proshade_double *&  ret 
)

Function to get array median and inter-quartile range.

This function takes a pointer to a array of proshade_double's and returns the median and the inter-quartile range of such vector.

Parameters
[in]vecPointer to an array of proshade_double's for which median and IQR should be obtained.
[in]vecSizeThe length of the array.
[in]retPointer to array of 2 proshade_double's, which will be the return values - first median and second IQR.

Definition at line 198 of file ProSHADE_maths.cpp.

199 {
200  //================================================ Sort the vector
201  std::sort ( vec, vec + vecSize );
202 
203  //================================================ Get median
204  if ( vecSize % 2 == 0)
205  {
206  ret[0] = ( vec[ ( vecSize / 2 ) - 1 ] + vec[ vecSize / 2 ] ) / 2.0;
207  }
208  else
209  {
210  ret[0] = vec[ vecSize / 2 ];
211  }
212 
213  //================================================ Get first and third quartile
214  proshade_double Q1, Q3;
215  if ( vecSize % 2 == 0)
216  {
217  Q1 = ( vec[ ( vecSize / 4 ) - 1 ] + vec[ vecSize / 4 ] ) / 2.0;
218  Q3 = ( vec[ ( ( vecSize / 4 ) * 3 ) - 1 ] + vec[ ( vecSize / 4 ) * 3 ] ) / 2.0;
219  }
220  else
221  {
222  Q1 = vec[ vecSize / 4 ];
223  Q3 = vec[ ( vecSize / 4 ) * 3 ];
224  }
225 
226  //================================================ And now save the IQR
227  ret[1] = Q3 - Q1;
228 
229  //================================================ Return
230  return ;
231 
232 }

◆ completeLegendreSeries()

void ProSHADE_internal_maths::completeLegendreSeries ( proshade_unsign  order,
proshade_double *  abscissas,
proshade_double *  weights,
proshade_unsign  taylorSeriesCap 
)

This function completes the Legendre polynomial series assuming you have obtained the first values.

Given that the polynomial value at zero is known, this function will complete the Legendre polynomial and with it the absicassas and weights for the Gauss-Legendre integration using the other functions defined above.

Parameters
[in]orderThe positive integer value of the polynomial order.
[in]abscissasPointer to an array of abscissas containing the first value.
[in]weightsPointer to an array of weights containing the first value.
[in]taylorSeriesCapThe limit on the Taylor series.

Definition at line 521 of file ProSHADE_maths.cpp.

522 {
523  //================================================ Initialise internal variables
524  proshade_double hlpTaylorVal = 0.0;
525  proshade_double hlpOrderVal = static_cast<proshade_double> ( order );
526  proshade_double abscValueChange = 0.0;
527  proshade_double prevAbsc = 0.0;
528  proshade_double *hlpAbscSeries;
529  proshade_double *hlpWeightSeries;
530  proshade_unsign noSeriesElems = 0;
531  proshade_unsign oddEvenSwitch = 0;
532 
533  //================================================ Pre-set internal values
534  if ( order % 2 == 1 )
535  {
536  noSeriesElems = ( order - 1 ) / 2 - 1;
537  oddEvenSwitch = 1;
538  }
539  else
540  {
541  noSeriesElems = order / 2 - 1;
542  oddEvenSwitch = 0;
543  }
544 
545  //================================================ Allocate memory
546  hlpAbscSeries = new proshade_double[taylorSeriesCap+2];
547  hlpWeightSeries = new proshade_double[taylorSeriesCap+1];
548 
549  //================================================ For each series element
550  for ( proshade_unsign serIt = noSeriesElems + 1; serIt < order - 1; serIt++ )
551  {
552  //============================================ Init loop
553  prevAbsc = abscissas[serIt];
554  abscValueChange = advanceGLPolyValue ( M_PI/2.0, -M_PI/2.0, prevAbsc, order, taylorSeriesCap ) - prevAbsc;
555 
556  //============================================ Init abscissas
557  hlpAbscSeries[0] = 0.0;
558  hlpAbscSeries[1] = 0.0;
559  hlpAbscSeries[2] = weights[serIt];
560 
561  //============================================ Init weights
562  hlpWeightSeries[0] = 0.0;
563  hlpWeightSeries[1] = hlpAbscSeries[2];
564 
565  //============================================ Taylor expansion
566  for ( proshade_unsign tayIt = 0; tayIt <= taylorSeriesCap - 2; tayIt++ )
567  {
568  hlpTaylorVal = static_cast<proshade_double> ( tayIt );
569 
570  hlpAbscSeries[tayIt+3] = ( 2.0 * prevAbsc * ( hlpTaylorVal + 1.0 ) * hlpAbscSeries[tayIt+2] + ( hlpTaylorVal * ( hlpTaylorVal + 1.0 ) - hlpOrderVal *
571  ( hlpOrderVal + 1.0 ) ) * hlpAbscSeries[tayIt+1] / ( hlpTaylorVal + 1.0 ) ) / ( 1.0 - prevAbsc ) / ( 1.0 + prevAbsc ) /
572  ( hlpTaylorVal + 2.0 );
573 
574  hlpWeightSeries[tayIt+2] = ( hlpTaylorVal + 2.0 ) * hlpAbscSeries[tayIt+3];
575  }
576 
577  //============================================ Sum over results
578  for ( proshade_unsign iter = 0; iter < 5; iter++ )
579  {
580  abscValueChange = abscValueChange - evaluateGLSeries ( hlpAbscSeries, abscValueChange, taylorSeriesCap ) /
581  evaluateGLSeries ( hlpWeightSeries, abscValueChange, taylorSeriesCap-1 );
582  }
583 
584  //============================================ Save results
585  abscissas[serIt+1] = prevAbsc + abscValueChange;
586  weights[serIt+1] = evaluateGLSeries ( hlpWeightSeries, abscValueChange, taylorSeriesCap - 1 );
587  }
588 
589  for ( proshade_unsign serIt = 0; serIt <= noSeriesElems + oddEvenSwitch; serIt++ )
590  {
591  abscissas[serIt] = -abscissas[order-serIt-1];
592  weights[serIt] = weights[order-serIt-1];
593  }
594 
595  //================================================ Free memory
596  delete hlpAbscSeries;
597  delete hlpWeightSeries;
598 
599  //================================================ Done
600  return ;
601 
602 }

◆ complexMatrixSVDSigmasOnly()

void ProSHADE_internal_maths::complexMatrixSVDSigmasOnly ( proshade_complex **  mat,
int  dim,
double *&  singularValues 
)

Function to compute the complete complex matrix SVD and return only the sigmas.

This function converts the input proshade_complex matrix of dimensions dim onto the LAPACK compatible std::complex<double> matrix. It then proceeds to create a dummy variables for the U and V matrices for saving the SVD results as well as other required variables. It finally proceeds to call LAPACK ZGESDD function to compute the SVD of the complex matrix input, checks the results and terminates. Note that this function does not make use of most of the LAPACK capabilities and is limitted onto square matrices.

Parameters
[in]matPointer to a complex square matrix with dimensions dim * dim.
[in]dimThe dimension of the complex matrix.
[in]singularValuesEmpty array of size dim where the singular values will be saved.

Definition at line 802 of file ProSHADE_maths.cpp.

803 {
804  //================================================ Initialise local variables
805  char job = 'N'; // Save computation of parts of U and V matrices, they are not needed here
806  std::complex<double> *rotMatU = new std::complex<double> [dim*dim]; // The U matrix space
807  std::complex<double> *rotMatV = new std::complex<double> [dim*dim]; // The V^T matrix space
808  std::complex<double> *work = new std::complex<double> [( 4 * dim)]; // Workspace, minimum required is 3*dim, using more for performance
809  int workDim = ( 4 * dim); // Formalism stating just that
810  double* rwork = new double[(7 * dim)]; // Required by LAPACK, from 3.7 requires 7 * dim
811  int* iwork = new int[(8 * dim)]; // Required by LAPACK
812  int returnValue = 0; // This will tell if operation succeeded
813 
814  //================================================ Check memory allocation
815  ProSHADE_internal_misc::checkMemoryAllocation ( rotMatU, __FILE__, __LINE__, __func__ );
816  ProSHADE_internal_misc::checkMemoryAllocation ( rotMatV, __FILE__, __LINE__, __func__ );
817  ProSHADE_internal_misc::checkMemoryAllocation ( work, __FILE__, __LINE__, __func__ );
818  ProSHADE_internal_misc::checkMemoryAllocation ( rwork, __FILE__, __LINE__, __func__ );
819  ProSHADE_internal_misc::checkMemoryAllocation ( iwork, __FILE__, __LINE__, __func__ );
820 
821  //================================================ Load input data into array in column-major order
822  std::complex<double> *matrixToDecompose = new std::complex<double>[dim*dim];
823  ProSHADE_internal_misc::checkMemoryAllocation ( matrixToDecompose, __FILE__, __LINE__, __func__ );
824  for ( int rowIt = 0; rowIt < dim; rowIt++ )
825  {
826  for ( int colIt = 0; colIt < dim; colIt++ )
827  {
828  matrixToDecompose[(colIt*dim)+rowIt] = std::complex<double> ( mat[rowIt][colIt][0], mat[rowIt][colIt][1] );
829  }
830  }
831 
832  //================================================ Run LAPACK ZGESDD
833  zgesdd_ ( &job, &dim, &dim, matrixToDecompose, &dim, singularValues, rotMatU, &dim, rotMatV, &dim,
834  work, &workDim, rwork, iwork, &returnValue );
835 
836  //================================================ Free memory
837  delete[] rotMatU;
838  delete[] rotMatV;
839  delete[] work;
840  delete[] rwork;
841  delete[] iwork;
842  delete[] matrixToDecompose;
843 
844  //================================================ Check result
845  if ( returnValue != 0 )
846  {
847  throw ProSHADE_exception ( "The LAPACK complex SVD algorithm did not converge!", "EL00021", __FILE__, __LINE__, __func__, "LAPACK algorithm for computing the singular value\n : decomposition of complex matrices did not converge and\n : therefore it was not possible to combine SH coefficients\n : from multiple shells. Changing the resolution may help,\n : contact me if this error persists." );
848  }
849 
850  //================================================ Done
851  return ;
852 
853 }

◆ complexMatrixSVDUandVOnly()

void ProSHADE_internal_maths::complexMatrixSVDUandVOnly ( proshade_double *  mat,
int  dim,
proshade_double *  uAndV,
bool  fail = true 
)

Function to compute the real matrix SVD and return the U and V matrices.

This function converts the input proshade_double array of dimensions dim*dim onto the LAPACK compatible std::complex<double> matrix. It then proceeds to create a dummy variables for the U and V matrices for saving the SVD results as well as other required variables. It finally proceeds to call LAPACK ZGESDD function to compute the SVD of the real matrix input, checks the results and saves the U and V matrices for output. Note that this function does not make use of most of the LAPACK capabilities and is limitted onto square matrices.

Parameters
[in]matPointer to a real square matrix with dimensions dim * dim.
[in]dimThe dimension of the real matrix.
[in]uAndVEmpty and allocated array of size dim*6 where the U and V matrices will be saved.
[in]failIf true and an error is encountered (typically algorithm not converging), this function will stop the program (useful for distances computations). However, if false, the function will simply return -777 as the first matrix element and not fail.

Definition at line 869 of file ProSHADE_maths.cpp.

870 {
871  //================================================ Initialise local variables
872  char job = 'A'; // Save computation of parts of U and V matrices, they are not needed here
873  double* singularValues = new double[dim]; // The array of singular values
874  std::complex<double> *rotMatU = new std::complex<double> [dim*dim]; // The U matrix space
875  std::complex<double> *rotMatV = new std::complex<double> [dim*dim]; // The V^T matrix space
876  std::complex<double> *work = new std::complex<double> [static_cast<proshade_unsign>( ( 3 * dim) + pow( dim, 2 ) * dim)]; // Workspace, minimum required is 3*dim, using more for performance
877  int workDim = ( 3 * dim) + pow( dim, 2 ); // Formalism stating just that
878  double* rwork = new double[static_cast<proshade_unsign>((5 * dim) + 5 * pow(dim,2))]; // Required by LAPACK
879  int* iwork = new int[(8 * dim)]; // Required by LAPACK
880  int returnValue = 0; // This will tell if operation succeeded
881  ProSHADE_internal_misc::checkMemoryAllocation ( singularValues, __FILE__, __LINE__, __func__ );
882  ProSHADE_internal_misc::checkMemoryAllocation ( rotMatU, __FILE__, __LINE__, __func__ );
883  ProSHADE_internal_misc::checkMemoryAllocation ( rotMatV, __FILE__, __LINE__, __func__ );
884  ProSHADE_internal_misc::checkMemoryAllocation ( work, __FILE__, __LINE__, __func__ );
885  ProSHADE_internal_misc::checkMemoryAllocation ( rwork, __FILE__, __LINE__, __func__ );
886  ProSHADE_internal_misc::checkMemoryAllocation ( iwork, __FILE__, __LINE__, __func__ );
887 
888  //================================================ Load input data into array in column-major order
889  std::complex<double> *matrixToDecompose = new std::complex<double>[dim*dim];
890  ProSHADE_internal_misc::checkMemoryAllocation ( matrixToDecompose, __FILE__, __LINE__, __func__ );
891  for ( int rowIt = 0; rowIt < dim; rowIt++ )
892  {
893  for ( int colIt = 0; colIt < dim; colIt++ )
894  {
895  matrixToDecompose[(colIt*dim)+rowIt] = std::complex<double> ( mat[(rowIt*dim)+colIt], 0.0 );
896  }
897  }
898 
899  //================================================ Run LAPACK ZGESDD
900  zgesdd_ ( &job, &dim, &dim, matrixToDecompose, &dim, singularValues, rotMatU, &dim, rotMatV, &dim,
901  work, &workDim, rwork, iwork, &returnValue );
902 
903  //================================================ Free memory
904  delete[] work;
905  delete[] rwork;
906  delete[] iwork;
907  delete[] matrixToDecompose;
908  delete[] singularValues;
909 
910  //================================================ Check result
911  if ( ( returnValue != 0 ) && ( fail ) )
912  {
913  throw ProSHADE_exception ( "The LAPACK complex SVD algorithm did not converge!", "EL00022", __FILE__, __LINE__, __func__, "LAPACK algorithm for computing the singular value\n : decomposition of complex matrices did not converge and\n : therefore it was not possible to optimise the peak\n : positions in the (self-)rotation function. Changing the\n : resolution may help, contact me if this error persists." );
914  }
915  if ( ( returnValue != 0 ) && ( !fail ) )
916  {
917  uAndV[0] = -777.7;
918  return ;
919  }
920 
921  //================================================ Save U
922  for ( proshade_signed rowIt = 0; rowIt < dim; rowIt++ )
923  {
924  for ( proshade_signed colIt = 0; colIt < dim; colIt++ )
925  {
926  uAndV[(rowIt*3)+colIt] = rotMatU[( rowIt * 3 ) + colIt].real();
927  }
928  }
929 
930  //================================================ Save V
931  for ( proshade_signed rowIt = 0; rowIt < dim; rowIt++ )
932  {
933  for ( proshade_signed colIt = 0; colIt < dim; colIt++ )
934  {
935  uAndV[(rowIt*3)+colIt+9] = rotMatV[( rowIt * 3 ) + colIt].real();
936  }
937  }
938 
939  //================================================ Release the rest of the memory
940  delete[] rotMatU;
941  delete[] rotMatV;
942 
943  //================================================ Done
944  return ;
945 
946 }

◆ complexMultiplication()

void ProSHADE_internal_maths::complexMultiplication ( proshade_double *  r1,
proshade_double *  i1,
proshade_double *  r2,
proshade_double *  i2,
proshade_double *  retReal,
proshade_double *  retImag 
)

Function to multiply two complex numbers.

This function takes pointers to the real and imaginary parts of two complex numbers and returns the result of their multiplication.

Parameters
[in]r1Pointer to the real value of number 1.
[in]i1Pointer to the imaginary value of number 1.
[in]r2Pointer to the real value of number 2.
[in]i2Pointer to the imaginary value of number 2.
[in]retRealPointer to the real part of the complex variable to which the result will be saved.
[in]retImagPointer to the imaginary part of the complex variable to which the result will be saved.

Definition at line 38 of file ProSHADE_maths.cpp.

39 {
40  //================================================ Multiplication
41  *retReal = (*r1)*(*r2) - (*i1)*(*i2);
42  *retImag = (*r1)*(*i2) + (*i1)*(*r2);
43 
44  //================================================ Return
45  return ;
46 
47 }

◆ complexMultiplicationConjug()

void ProSHADE_internal_maths::complexMultiplicationConjug ( proshade_double *  r1,
proshade_double *  i1,
proshade_double *  r2,
proshade_double *  i2,
proshade_double *  retReal,
proshade_double *  retImag 
)

Function to multiply two complex numbers by using the second number's conjugate.

This function takes pointers to the real and imaginary parts of two complex numbers and returns the result of their multiplication, while using the conjugate of the second complex number.

Parameters
[in]r1Pointer to the real value of number 1.
[in]i1Pointer to the imaginary value of number 1.
[in]r2Pointer to the real value of number 2.
[in]i2Pointer to the imaginary value of number 2.
[in]retRealPointer to the real part of the complex variable to which the result will be saved.
[in]retImagPointer to the imaginary part of the complex variable to which the result will be saved.

Definition at line 62 of file ProSHADE_maths.cpp.

63 {
64  //================================================ Multiplication
65  *retReal = (*r1)*(*r2) + (*i1)*(*i2);
66  *retImag = -(*r1)*(*i2) + (*i1)*(*r2);
67 
68  //================================================ Return
69  return ;
70 
71 }

◆ complexMultiplicationConjugRealOnly()

proshade_double ProSHADE_internal_maths::complexMultiplicationConjugRealOnly ( proshade_double *  r1,
proshade_double *  i1,
proshade_double *  r2,
proshade_double *  i2 
)

Function to conjuggate multiply two complex numbers and return the real part only.

This function takes pointers to the real and imaginary parts of two complex numbers and returns the real part of the result of their conjugate multiplication.

Parameters
[in]r1Pointer to the real value of number 1.
[in]i1Pointer to the imaginary value of number 1.
[in]r2Pointer to the real value of number 2.
[in]i2Pointer to the imaginary value of number 2.

Definition at line 103 of file ProSHADE_maths.cpp.

104 {
105  //================================================ Multiplication
106  proshade_double ret = (*r1)*(*r2) + (*i1)*(*i2);
107 
108  //================================================ Return
109  return ( ret );
110 
111 }

◆ complexMultiplicationRealOnly()

proshade_double ProSHADE_internal_maths::complexMultiplicationRealOnly ( proshade_double *  r1,
proshade_double *  i1,
proshade_double *  r2,
proshade_double *  i2 
)

Function to multiply two complex numbers and return the real part only.

This function takes pointers to the real and imaginary parts of two complex numbers and returns the real part of the result of their multiplication.

Parameters
[in]r1Pointer to the real value of number 1.
[in]i1Pointer to the imaginary value of number 1.
[in]r2Pointer to the real value of number 2.
[in]i2Pointer to the imaginary value of number 2.

Definition at line 83 of file ProSHADE_maths.cpp.

84 {
85  //================================================ Multiplication
86  proshade_double ret = (*r1)*(*r2) - (*i1)*(*i2);
87 
88  //================================================ Return
89  return ( ret );
90 
91 }

◆ compute3x3MatrixInverse()

proshade_double * ProSHADE_internal_maths::compute3x3MatrixInverse ( proshade_double *  mat)

Function for computing a 3x3 matrix inverse.

Parameters
[in]matThe matrix to be inverted.
[out]inverseThe inverse of matrix mat.

Definition at line 1810 of file ProSHADE_maths.cpp.

1811 {
1812  //================================================ Allocate memory
1813  proshade_double* inverse = new proshade_double[9];
1814  ProSHADE_internal_misc::checkMemoryAllocation ( inverse, __FILE__, __LINE__, __func__ );
1815 
1816  //================================================ Compute determinant
1817  proshade_double matDet = ( mat[0] * mat[4] * mat[8] ) +
1818  ( mat[1] * mat[5] * mat[6] ) +
1819  ( mat[2] * mat[3] * mat[7] ) -
1820  ( mat[0] * mat[5] * mat[7] ) -
1821  ( mat[1] * mat[3] * mat[8] ) -
1822  ( mat[2] * mat[4] * mat[6] );
1823 
1824  //================================================ Compute inverse matrix
1825  inverse[0] = ( mat[4] * mat[8] - mat[5] * mat[7] ) / matDet;
1826  inverse[1] = ( mat[2] * mat[7] - mat[1] * mat[8] ) / matDet;
1827  inverse[2] = ( mat[1] * mat[5] - mat[2] * mat[4] ) / matDet;
1828  inverse[3] = ( mat[5] * mat[6] - mat[3] * mat[8] ) / matDet;
1829  inverse[4] = ( mat[0] * mat[8] - mat[2] * mat[6] ) / matDet;
1830  inverse[5] = ( mat[2] * mat[3] - mat[0] * mat[5] ) / matDet;
1831  inverse[6] = ( mat[3] * mat[7] - mat[4] * mat[6] ) / matDet;
1832  inverse[7] = ( mat[1] * mat[6] - mat[0] * mat[7] ) / matDet;
1833  inverse[8] = ( mat[0] * mat[4] - mat[1] * mat[3] ) / matDet;
1834 
1835  //================================================ Done
1836  return ( inverse );
1837 }

◆ compute3x3MatrixMultiplication()

proshade_double * ProSHADE_internal_maths::compute3x3MatrixMultiplication ( proshade_double *  mat1,
proshade_double *  mat2 
)

Function for computing a 3x3 matrix multiplication.

Parameters
[in]mat1The matrix to multiply mat2.
[in]mat2The matrix to be multiplied by mat1.
[out]retThe matrix resulting from matrix multiplication of mat1 and mat2 in this order.

Definition at line 1759 of file ProSHADE_maths.cpp.

1760 {
1761  //================================================ Allocate memory
1762  proshade_double* ret = new proshade_double[9];
1763  ProSHADE_internal_misc::checkMemoryAllocation ( ret, __FILE__, __LINE__, __func__ );
1764 
1765  //================================================ Multiply
1766  ret[0] = ( mat1[0] * mat2[0] ) + ( mat1[1] * mat2[3] ) + ( mat1[2] * mat2[6] );
1767  ret[1] = ( mat1[0] * mat2[1] ) + ( mat1[1] * mat2[4] ) + ( mat1[2] * mat2[7] );
1768  ret[2] = ( mat1[0] * mat2[2] ) + ( mat1[1] * mat2[5] ) + ( mat1[2] * mat2[8] );
1769  ret[3] = ( mat1[3] * mat2[0] ) + ( mat1[4] * mat2[3] ) + ( mat1[5] * mat2[6] );
1770  ret[4] = ( mat1[3] * mat2[1] ) + ( mat1[4] * mat2[4] ) + ( mat1[5] * mat2[7] );
1771  ret[5] = ( mat1[3] * mat2[2] ) + ( mat1[4] * mat2[5] ) + ( mat1[5] * mat2[8] );
1772  ret[6] = ( mat1[6] * mat2[0] ) + ( mat1[7] * mat2[3] ) + ( mat1[8] * mat2[6] );
1773  ret[7] = ( mat1[6] * mat2[1] ) + ( mat1[7] * mat2[4] ) + ( mat1[8] * mat2[7] );
1774  ret[8] = ( mat1[6] * mat2[2] ) + ( mat1[7] * mat2[5] ) + ( mat1[8] * mat2[8] );
1775 
1776  //================================================ Done
1777  return ( ret );
1778 
1779 }

◆ compute3x3MatrixVectorMultiplication()

proshade_double * ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( proshade_double *  mat,
proshade_double  x,
proshade_double  y,
proshade_double  z 
)

Function for computing a 3x3 matrix to 3x1 vector multiplication.

Parameters
[in]matThe matrix to multiply the vector with..
[in]xThe x-axis element of the vector which is to be multiplied by the matrix.
[in]yThe x-axis element of the vector which is to be multiplied by the matrix.
[in]zThe x-axis element of the vector which is to be multiplied by the matrix.
[out]retThe vector resulting from matrix multiplication of mat and the vector in this order.

Definition at line 1789 of file ProSHADE_maths.cpp.

1790 {
1791  //================================================ Allocate memory
1792  proshade_double* ret = new proshade_double[3];
1793  ProSHADE_internal_misc::checkMemoryAllocation ( ret, __FILE__, __LINE__, __func__ );
1794 
1795  //================================================ Compute the multiplication
1796  ret[0] = ( x * mat[0] ) + ( y * mat[1] ) + ( z * mat[2] );
1797  ret[1] = ( x * mat[3] ) + ( y * mat[4] ) + ( z * mat[5] );
1798  ret[2] = ( x * mat[6] ) + ( y * mat[7] ) + ( z * mat[8] );
1799 
1800  //================================================ Done
1801  return ( ret );
1802 
1803 }

◆ computeCrossProduct()

proshade_double * ProSHADE_internal_maths::computeCrossProduct ( proshade_double *  x1,
proshade_double *  y1,
proshade_double *  z1,
proshade_double *  x2,
proshade_double *  y2,
proshade_double *  z2 
)

Simple 3D vector cross product computation.

Parameters
[in]x1The x-axis element of the first vector.
[in]y1The y-axis element of the first vector.
[in]z1The z-axis element of the first vector.
[in]x2The x-axis element of the second vector.
[in]y2The y-axis element of the second vector.
[in]z2The z-axis element of the second vector.
[out]crossProdThe vector representing the cross product of the two input vectors.

Definition at line 1737 of file ProSHADE_maths.cpp.

1738 {
1739  //================================================ Allocate memory
1740  proshade_double* crossProd = new proshade_double[3];
1741  ProSHADE_internal_misc::checkMemoryAllocation ( crossProd, __FILE__, __LINE__, __func__ );
1742 
1743  //================================================ Compute
1744  crossProd[0] = ( (*y1) * (*z2) ) - ( (*z1) * (*y2) );
1745  crossProd[1] = ( (*z1) * (*x2) ) - ( (*x1) * (*z2) );
1746  crossProd[2] = ( (*x1) * (*y2) ) - ( (*y1) * (*x2) );
1747 
1748  //================================================ Done
1749  return ( crossProd );
1750 
1751 }

◆ computeDotProduct() [1/2]

proshade_double ProSHADE_internal_maths::computeDotProduct ( proshade_double *  x1,
proshade_double *  y1,
proshade_double *  z1,
proshade_double *  x2,
proshade_double *  y2,
proshade_double *  z2 
)

Simple 3D vector dot product computation.

Parameters
[in]x1The x-axis element of the first vector.
[in]y1The y-axis element of the first vector.
[in]z1The z-axis element of the first vector.
[in]x2The x-axis element of the second vector.
[in]y2The y-axis element of the second vector.
[in]z2The z-axis element of the second vector.
[out]XThe dot product of the two input vectors.

Definition at line 1705 of file ProSHADE_maths.cpp.

1706 {
1707  //================================================ Compute and return
1708  return ( (*x1 * *x2) + (*y1 * *y2) + (*z1 * *z2) );
1709 }

◆ computeDotProduct() [2/2]

proshade_double ProSHADE_internal_maths::computeDotProduct ( proshade_double  x1,
proshade_double  y1,
proshade_double  z1,
proshade_double  x2,
proshade_double  y2,
proshade_double  z2 
)

Simple 3D vector dot product computation.

Parameters
[in]x1The x-axis element of the first vector.
[in]y1The y-axis element of the first vector.
[in]z1The z-axis element of the first vector.
[in]x2The x-axis element of the second vector.
[in]y2The y-axis element of the second vector.
[in]z2The z-axis element of the second vector.
[out]XThe dot product of the two input vectors.

Definition at line 1721 of file ProSHADE_maths.cpp.

1722 {
1723  //================================================ Compute and return
1724  return ( (x1 * x2) + (y1 * y2) + (z1 * z2) );
1725 }

◆ evaluateGLSeries()

proshade_double ProSHADE_internal_maths::evaluateGLSeries ( proshade_double *  series,
proshade_double  target,
proshade_unsign  terms 
)

This function evaluates the Taylor expansion.

This function takes the series array, the target value and the cap on Taylor expansion and proceeds to evaluate the series. The main use of this is to evaluate the series twice, one where the series evaluation 'overshoots' and once where it 'undershoots' and taking value in between those, thus adding accuracy.

Parameters
[in]seriesPointer to array with the series values.
[in]targetThe target location on the series value.
[in]termsThe Taylor expansion cap.
[out]XThe value of the series at the target location.

Definition at line 447 of file ProSHADE_maths.cpp.

448 {
449  //================================================ Initalise
450  proshade_double factorialValue = 1.0;
451  proshade_double value = 0.0;
452 
453  //================================================ Compute
454  for ( proshade_unsign iter = 1; iter <= terms; iter++ )
455  {
456  value = value + series[iter] * factorialValue;
457  factorialValue = factorialValue * target;
458  }
459 
460  //================================================ Done
461  return ( value );
462 
463 }

◆ findAllPrimes()

std::vector< proshade_unsign > ProSHADE_internal_maths::findAllPrimes ( proshade_unsign  upTo)

This function finds all prime numbers up to the supplied limit.

This function uses the sieve of Eratosthenes algorithm to find all prime numbers from 2 to the supplied limit. This is not the fastest algorithm and it may become slow when the limit is high, but it is fine for small numbers and given that we will use it for symmetry folds, which should not got much over 20, this should be more than fast enough.

Parameters
[in]upToThe limit to which prime numbers should be sought.

Definition at line 2693 of file ProSHADE_maths.cpp.

2694 {
2695  //================================================ Initialise variables
2696  std::vector< proshade_unsign > ret;
2697  std::vector< std::pair< proshade_unsign, bool > > sieveOfEratosthenesArray;
2698 
2699  //================================================ Sanity check
2700  if ( upTo < 2 ) { return ( ret ); }
2701 
2702  //================================================ Initialise the sieve array up to the required number
2703  for ( proshade_unsign iter = 2; iter <= upTo; iter++ ) { sieveOfEratosthenesArray.emplace_back ( std::pair< proshade_unsign, bool > ( iter, true ) ); }
2704 
2705  //================================================ For each entry in the array
2706  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); iter++ )
2707  {
2708  //============================================ If this entry is still true
2709  if ( sieveOfEratosthenesArray.at(iter).second )
2710  {
2711  //======================================== Set all entries with the position x * [this entry value] to false
2712  for ( proshade_unsign it = iter + sieveOfEratosthenesArray.at(iter).first; it < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); it += sieveOfEratosthenesArray.at(iter).first )
2713  {
2714  sieveOfEratosthenesArray.at(it).second = false;
2715  }
2716  }
2717  }
2718 
2719  //================================================ Copy passing results to return vector
2720  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); iter++ )
2721  {
2722  if ( sieveOfEratosthenesArray.at(iter).second ) { ProSHADE_internal_misc::addToUnsignVector ( &ret, sieveOfEratosthenesArray.at(iter).first ); }
2723  }
2724 
2725  //================================================ Done
2726  return ( ret );
2727 
2728 }

◆ findRotMatMatchingVectors()

proshade_double * ProSHADE_internal_maths::findRotMatMatchingVectors ( proshade_double  x1,
proshade_double  y1,
proshade_double  z1,
proshade_double  x2,
proshade_double  y2,
proshade_double  z2 
)

Computation of rotation matrix rotating one vector onto the other.

This function starts by normalising both input vectors to have magnitude of 1.0. Then, it computes the cosine and sine of the angle between the two vectors (using the magnitude of the cross product and the dot product); these are then sufficient to build the rotation matrix for rotation in plane on which both of the vectors lie.

It then proceeds to compute the change of basis matrix and its inverse, which are in turn sufficient to to compute the rotation matrix in the original basis. This rotation matrix is then returned.

Parameters
[in]x1The x-axis element of the first vector.
[in]y1The y-axis element of the first vector.
[in]z1The z-axis element of the first vector.
[in]x2The x-axis element of the second vector.
[in]y2The y-axis element of the second vector.
[in]z2The z-axis element of the second vector.
[out]rotMatRotation matrix optimally rotating x1 ; y1 ; z1 to match x2 ; y2 ; z2.

Definition at line 1883 of file ProSHADE_maths.cpp.

1884 {
1885  //================================================ Allocate required memory
1886  proshade_double* inPlaneRotation = new proshade_double[9];
1887  proshade_double* basisChangeMat = new proshade_double[9];
1888  ProSHADE_internal_misc::checkMemoryAllocation ( inPlaneRotation, __FILE__, __LINE__, __func__ );
1889  ProSHADE_internal_misc::checkMemoryAllocation ( basisChangeMat, __FILE__, __LINE__, __func__ );
1890 
1891  //================================================ Normalise inputs
1892  proshade_double normF = std::sqrt( std::pow( x1, 2.0 ) + std::pow ( y1, 2.0 ) + std::pow ( z1, 2.0 ) );
1893  x1 /= normF; y1 /= normF; z1 /= normF;
1894 
1895  normF = std::sqrt( std::pow( x2, 2.0 ) + std::pow ( y2, 2.0 ) + std::pow ( z2, 2.0 ) );
1896  x2 /= normF; y2 /= normF; z2 /= normF;
1897 
1898  //================================================ Compute cross product's magnitude
1899  proshade_double* crossProd = ProSHADE_internal_maths::computeCrossProduct( &x1, &y1, &z1, &x2, &y2, &z2 );
1900  proshade_double crossProdMag = std::sqrt( std::pow( crossProd[0], 2.0 ) + std::pow ( crossProd[1], 2.0 ) + std::pow ( crossProd[2], 2.0 ) );
1901  delete[] crossProd;
1902 
1903  //================================================ Compute dot product
1904  proshade_double dotProd = ProSHADE_internal_maths::computeDotProduct ( &x1, &y1, &z1, &x2, &y2, &z2 );
1905 
1906  //================================================ Construct the in-plane rotation matrix
1907  inPlaneRotation[0] = dotProd; inPlaneRotation[1] = -crossProdMag; inPlaneRotation[2] = 0.0;
1908  inPlaneRotation[3] = crossProdMag; inPlaneRotation[4] = dotProd; inPlaneRotation[5] = 0.0;
1909  inPlaneRotation[6] = 0.0; inPlaneRotation[7] = 0.0; inPlaneRotation[8] = 1.0;
1910 
1911  //================================================ Construct change of basis matrix
1912  crossProd = ProSHADE_internal_maths::computeCrossProduct( &x2, &y2, &z2, &x1, &y1, &z1 );
1913  normF = std::sqrt ( std::pow ( x2 - ( dotProd * x1 ), 2.0 ) + std::pow ( y2 - ( dotProd * y1 ), 2.0 ) + std::pow ( z2 - ( dotProd * z1 ), 2.0 ) );
1914 
1915  basisChangeMat[0] = x1; basisChangeMat[1] = ( x2 - ( dotProd * x1 ) ) / normF; basisChangeMat[2] = crossProd[0];
1916  basisChangeMat[3] = y1; basisChangeMat[4] = ( y2 - ( dotProd * y1 ) ) / normF; basisChangeMat[5] = crossProd[1];
1917  basisChangeMat[6] = z1; basisChangeMat[7] = ( z2 - ( dotProd * z1 ) ) / normF; basisChangeMat[8] = crossProd[2];
1918 
1919  //================================================ Invert the change of basis matrix
1920  proshade_double* basisChangeMatInverse = ProSHADE_internal_maths::compute3x3MatrixInverse ( basisChangeMat );
1921 
1922  //================================================ Multiply inverse of change of basis matrix with the in plane rotation matrix, then multiply the result with the inverse
1923  proshade_double* tmpMat = ProSHADE_internal_maths::compute3x3MatrixMultiplication ( basisChangeMat, inPlaneRotation );
1924  proshade_double* rotMat = ProSHADE_internal_maths::compute3x3MatrixMultiplication ( tmpMat, basisChangeMatInverse );
1925 
1926  //================================================ Release memory
1927  delete[] crossProd;
1928  delete[] inPlaneRotation;
1929  delete[] basisChangeMat;
1930  delete[] basisChangeMatInverse;
1931  delete[] tmpMat;
1932 
1933  //================================================ Done
1934  return ( rotMat );
1935 
1936 }

◆ findVectorFromThreeVAndThreeD()

std::vector< proshade_double > ProSHADE_internal_maths::findVectorFromThreeVAndThreeD ( proshade_double  x1,
proshade_double  y1,
proshade_double  z1,
proshade_double  x2,
proshade_double  y2,
proshade_double  z2,
proshade_double  x3,
proshade_double  y3,
proshade_double  z3,
proshade_double  dot1,
proshade_double  dot2,
proshade_double  dot3 
)

Function for finding a vector which would have a given three dot products to three other vectors.

This function takes three vectors and three dot product values. It then basically solves the following set of equations for x, y and z:

solX*x1 + solY*y1 + solZ*z1 = dot1
solX*x2 + solY*y2 + solZ*z2 = dot2
solX*x3 + solY*y3 + solZ*z3 = dot3

This should result in a vector, which has the required angles to all input vectors and is normalised (this is done later as part of this function). The equations are courtesy of https://www.wolframalpha.com/input/?i=Solve%5B%7Ba+x+%2B+b+y+%2B+c+z+%3D%3D+f%2C+u+x+%2B+v+y+%2B+w+z+%3D%3D+g%2C+k+x+%2B+l+y+%2B+m+z+%3D%3D+h%7D%2C+%7Bx%2C+y%2C+z%7D%5D webpage of Wolfram Alpha. If in doubt, do not fear to derive yourself :-).

Parameters
[in]x1The x-axis element of the first vector.
[in]y1The y-axis element of the first vector.
[in]z1The z-axis element of the first vector.
[in]x2The x-axis element of the second vector.
[in]y2The y-axis element of the second vector.
[in]z2The z-axis element of the second vector.
[in]dot1The dot product specifying the angle between the sought vector and the first input vector.
[in]dot2The dot product specifying the angle between the sought vector and the second input vectors.
[out]vecA std::vector containing the three elements of the sought vector.

Definition at line 2101 of file ProSHADE_maths.cpp.

2102 {
2103  //================================================ Initialise variables
2104  std::vector < proshade_double > ret;
2105 
2106  //================================================ Solution
2107  proshade_double solX = - ( y1 * dot2 * z3 - y1 * dot3 * z2 - z1 * dot2 * y3 + z1 * dot3 * y2 + dot1 * y3 * z2 - dot1 * z3 * y2 ) /
2108  ( -x1 * y3 * z2 + x1 * z3 * y2 + y1 * x3 * z2 - y1 * z3 * x2 - z1 * x3 * y2 + z1 * y3 * x2 );
2109  proshade_double solY = - ( x1 * dot2 * z3 - x1 * dot3 * z2 - z1 * dot2 * x3 + z1 * dot3 * x2 + dot1 * x3 * z2 - dot1 * z3 * x2 ) /
2110  ( x1 * y3 * z2 - x1 * z3 * y2 - y1 * x3 * z2 + y1 * z3 * x2 + z1 * x3 * y2 - z1 * y3 * x2 );
2111  proshade_double solZ = - ( x1 * dot2 * y3 - x1 * dot3 * y2 - y1 * dot2 * x3 + y1 * dot3 * x2 + dot1 * x3 * y2 - dot1 * y3 * x2 ) /
2112  ( -x1 * y3 * z2 + x1 * z3 * y2 + y1 * x3 * z2 - y1 * z3 * x2 - z1 * x3 * y2 + z1 * y3 * x2 );
2113 
2114  //================================================ Normalise the axis to magnitude 1
2115  proshade_double normFactor = sqrt ( pow ( solX, 2.0 ) + pow ( solY, 2.0 ) + pow ( solZ, 2.0 ) );
2116  solX /= normFactor;
2117  solY /= normFactor;
2118  solZ /= normFactor;
2119 
2120  //================================================ Set largest axis element to positive (ProSHADE standard)
2121  if ( ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solX ) ) && ( solX < 0.0 ) ) ||
2122  ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solY ) ) && ( solY < 0.0 ) ) ||
2123  ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solZ ) ) && ( solZ < 0.0 ) ) ) { solX *= -1.0; solY *= -1.0; solZ *= -1.0; }
2124 
2125  //================================================ Save solutions
2129 
2130  //================================================ Done
2131  return ( ret );
2132 
2133 }

◆ findVectorFromTwoVAndTwoD()

std::vector< proshade_double > ProSHADE_internal_maths::findVectorFromTwoVAndTwoD ( proshade_double  x1,
proshade_double  y1,
proshade_double  z1,
proshade_double  x2,
proshade_double  y2,
proshade_double  z2,
proshade_double  dot1,
proshade_double  dot2 
)

Function for finding a vector which would have a given two dot products to two other vectors.

This function takes two vectors and two dot product values. It then basically solves the following set of equations for x, y and z:

solX*x1 + solY*y1 + solZ*z1 = dot1
solX*x2 + solY*y2 + solZ*z2 = dot2
sqrt ( solX^2 + solY^2 + solZ^2 ) = 1.0

This should result in a vector, which has the required angles to both input vectors and is normalised (this is done by the third equation, which is required to obtain system of three equations with three unkonwns). The equations are courtesy of https://www.wolframalpha.com/input/?i=Solve%5B%7Ba+x+%2B+b+y+%2B+c+z+%3D%3D+f%2C+k+x+%2B+l+y+%2B+m+z+%3D%3D+g%2C+Sqrt%5Bx%5E2+%2B+y%5E2+%2B+z%5E2%5D+%3D%3D+1%7D%2C+%7Bx%2C+y%2C+z%7D%5D webpage of Wolfram Alpha. If in doubt, do not fear to derive yourself :-).

Parameters
[in]x1The x-axis element of the first vector.
[in]y1The y-axis element of the first vector.
[in]z1The z-axis element of the first vector.
[in]x2The x-axis element of the second vector.
[in]y2The y-axis element of the second vector.
[in]z2The z-axis element of the second vector.
[in]dot1The dot product specifying the angle between the sought vector and the first input vector.
[in]dot2The dot product specifying the angle between the sought vector and the second input vectors.
[out]vecA std::vector containing the three elements of the sought vector.

Definition at line 1959 of file ProSHADE_maths.cpp.

1960 {
1961  //================================================ Initialise variables
1962  std::vector < proshade_double > ret;
1963 
1964  //================================================ Solution
1965  proshade_double solX = ( -sqrt ( pow ( 2.0 * x1 * y1 * dot2 * y2 + 2.0 * x1 * z1 * dot2 * z2 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( y1, 2.0 ) * dot2 * x2 + 2.0 * y1 * dot1 * x2 * y2 - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
1966  4.0 * ( pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * x1 * y1 * x2 * y2 - 2.0 * x1 * z1 * x2 * z2 + pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) ) *
1967  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 - 2.0 * y1 * dot1 * dot2 * y2 + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) - 2.0 * z1 * dot1 * dot2 * z2 + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) ) ) -
1968  2.0 * x1 * y1 * dot2 * y2 - 2.0 * x1 * z1 * dot2 * z2 + 2.0 * x1 * dot1 * pow ( y2, 2.0 ) + 2.0 * x1 * dot1 * pow ( z2, 2.0 ) + 2.0 * pow ( y1, 2.0 ) * dot2 * x2 - 2.0 * y1 * dot1 * x2 * y2 + 2.0 * pow ( z1, 2.0 ) * dot2 * x2 - 2.0 * z1 * dot1 * x2 * z2 ) /
1969  ( 2.0 * ( pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * x1 * y1 * x2 * y2 - 2.0 * x1 * z1 * x2 * z2 + pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) ) );
1970  proshade_double solY = ( ( dot2 * pow ( x2, 2.0 ) * pow ( z1, 3.0 ) ) /
1971  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1972  ( dot1 * pow ( x2, 2.0 ) * z2 * pow ( z1, 2.0 ) ) /
1973  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1974  ( 2.0 * x1 * dot2 * x2 * z2 * pow ( z1, 2.0 ) ) /
1975  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) - dot2 * z1 -
1976  ( x2 * sqrt ( pow ( - 2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
1977  4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
1978  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) * z1 ) /
1979  ( 2.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
1980  ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) * z1 ) /
1981  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
1982  ( x1 * dot1 * x2 * pow ( y2, 2.0 ) * z1 ) /
1983  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
1984  ( pow ( x1, 2.0 ) * dot2 * pow ( z2, 2.0 ) * z1 ) /
1985  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
1986  ( 2.0 * x1 * dot1 * x2 * pow ( z2, 2.0 ) * z1 ) /
1987  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1988  ( y1 * dot1 * pow ( x2, 2.0 ) * y2 * z1 ) /
1989  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1990  ( x1 * y1 * dot2 * x2 * y2 * z1 ) /
1991  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) + dot1 * z2 +
1992  ( x1 * z2 * sqrt ( pow ( -2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
1993  4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
1994  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) ) /
1995  ( 2.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
1996  ( pow ( x1, 2.0 ) * dot1 * pow ( z2, 3.0 ) ) /
1997  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1998  ( pow ( x1, 2.0 ) * dot1 * pow ( y2, 2.0 ) * z2 ) /
1999  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2000  ( x1 * pow ( y1, 2.0 ) * dot2 * x2 * z2 ) /
2001  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2002  ( pow ( x1, 2.0 ) * y1 * dot2 * y2 * z2 ) /
2003  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2004  ( x1 * y1 * dot1 * x2 * y2 * z2 ) /
2005  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) / ( y1 * z2 - z1 * y2 );
2006  proshade_double solZ = ( - ( dot2 * pow ( x2, 2.0 ) * y2 * pow ( z1, 3.0 ) ) /
2007  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2008  ( dot2 * pow ( x2, 2.0 ) * pow ( z1, 2.0 ) ) /
2009  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2010  ( dot1 * pow ( x2, 2.0 ) * y2 * z2 * pow ( z1, 2.0 ) ) /
2011  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2012  ( 2.0 * x1 * dot2 * x2 * y2 * z2 * pow ( z1, 2.0 ) ) /
2013  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2014  ( x2 * y2 * sqrt ( pow ( -2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
2015  4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
2016  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) * z1 ) /
2017  ( 2.0 * ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2018  ( dot2 * y2 * z1 ) / ( y1 * z2 - z1 * y2 ) +
2019  ( dot1 * pow ( x2, 2.0 ) * z2 * z1 ) /
2020  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2021  ( x1 * dot2 * x2 * z2 * z1 ) /
2022  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2023  ( x1 * dot1 * x2 * pow ( y2, 3.0 ) * z1 ) /
2024  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2025  ( y1 * dot1 * pow ( x2, 2.0 ) * pow ( y2, 2.0 ) * z1 ) /
2026  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2027  ( x1 * y1 * dot2 * x2 * pow ( y2, 2.0 ) * z1 ) /
2028  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2029  ( pow ( x1, 2.0 ) * dot2 * y2 * pow ( z2, 2.0 ) * z1 ) /
2030  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2031  ( 2.0 * x1 * dot1 * x2 * y2 * pow ( z2, 2.0 ) * z1 ) /
2032  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2033  ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) * y2 * z1 ) /
2034  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) + dot2 +
2035  ( x2 * sqrt ( pow ( - 2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
2036  4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
2037  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) ) /
2038  ( 2.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2039  ( x1 * y2 * z2 * sqrt ( pow ( - 2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
2040  4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
2041  ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) ) /
2042  ( 2.0 * ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2043  ( dot1 * y2 * z2 ) / ( y1 * z2 - z1 * y2 ) -
2044  ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) ) /
2045  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2046  ( x1 * dot1 * x2 * pow ( y2, 2.0 ) ) /
2047  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2048  ( x1 * dot1 * x2 * pow ( z2, 2.0 ) ) /
2049  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2050  ( y1 * dot1 * pow ( x2, 2.0 ) * y2 ) /
2051  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2052  ( x1 * y1 * dot2 * x2 * y2 ) /
2053  ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2054  ( pow ( x1, 2.0 ) * dot1 * y2 * pow ( z2, 3.0 ) ) /
2055  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2056  ( pow ( x1, 2.0 ) * dot1 * pow ( y2, 3.0 ) * z2 ) /
2057  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2058  ( pow ( x1, 2.0 ) * y1 * dot2 * pow ( y2, 2.0 ) * z2 ) /
2059  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2060  ( x1 * y1 * dot1 * x2 * pow ( y2, 2.0 ) * z2 ) /
2061  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2062  ( x1 * pow ( y1, 2.0 ) * dot2 * x2 * y2 * z2 ) /
2063  ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) ) / z2;
2064 
2065  //================================================ Set largest axis element to positive (ProSHADE standard)
2066  if ( ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solX ) ) && ( solX < 0.0 ) ) ||
2067  ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solY ) ) && ( solY < 0.0 ) ) ||
2068  ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solZ ) ) && ( solZ < 0.0 ) ) ) { solX *= -1.0; solY *= -1.0; solZ *= -1.0; }
2069 
2070  //================================================ Save solutions
2074 
2075  //================================================ Done
2076  return ( ret );
2077 
2078 }

◆ gaussLegendreIntegration()

void ProSHADE_internal_maths::gaussLegendreIntegration ( proshade_complex *  vals,
proshade_unsign  valsSize,
proshade_unsign  order,
proshade_double *  abscissas,
proshade_double *  weights,
proshade_double  integralOverRange,
proshade_double  maxSphereDists,
proshade_double *  retReal,
proshade_double *  retImag 
)

Function to compute the complete complex Gauss-Legendre integration over spherical harmonic values in different shells.

This function takes the real parts of the spherical harmonics value in different shells and proceeds to compute the Gauss-Legendre integration over them. It uses the shell positions to appropriately place abscissas and their weights, which it assumes were pre-computed by the getLegendreAbscAndWeights() function.

Parameters
[in]valsPointer to a complex array of values over which the integration to be done.
[in]valsSizeThe length of the input array.
[in]orderThe integration order value.
[in]abscissasThe allocated array for holding the abscissa values.
[in]weightsThe allocated array for holding the weight values.
[in]integralOverRangeThe range of the intgral. If progressive shell mapping is used, this will not be max shell radius.
[in]maxSphereDistsDistance between two shells.
[in]retRealThe real part of the complex result of Gauss-Legendre integration over the shperical harmonics values.
[in]retImagThe imaginary part of the complex result of Gauss-Legendre integration over the shperical harmonics values.

Definition at line 709 of file ProSHADE_maths.cpp.

710 {
711  //================================================ Initialise local variables
712  proshade_triplet* intData = new proshade_triplet [order];
713  ProSHADE_internal_misc::checkMemoryAllocation ( intData, __FILE__, __LINE__, __func__ );
714  proshade_triplet posVals;
715  proshade_unsign lesserPos = 0;
716  proshade_unsign upperPos = 0;
717  proshade_double lesserWeight = 0.0;
718  proshade_double upperWeight = 0.0;
719 
720  //================================================ Rescale to <order> points
721  for ( proshade_unsign absIter = 0; absIter < order; absIter++ )
722  {
723  //============================================ Init loop
724  posVals[0] = 0.0;
725  posVals[1] = 0.0;
726  posVals[2] = 0.0;
727 
728  //============================================ Find real position of abscissas
729  posVals[0] = ( ( abscissas[absIter] + 1.0 ) / 2.0 ) * integralOverRange;
730 
731 
732  //============================================ Find lesser and upper bounds
733  for ( proshade_unsign valIt = 0; valIt < valsSize; valIt++ )
734  {
735  if ( ( ( valIt * maxSphereDists ) <= posVals[0] ) && ( ( ( valIt + 1 ) * maxSphereDists ) > posVals[0] ) )
736  {
737  lesserPos = static_cast<proshade_unsign> ( valIt );
738  upperPos = static_cast<proshade_unsign> ( valIt + 1 );
739  break;
740  }
741  }
742 
743  //============================================ Linear Interpolation
744  lesserWeight = 0.0;
745  upperWeight = 0.0;
746  if ( lesserPos != 0 )
747  {
748  //======================================== Here we realise that the lesser and upper bounds were determined on scale 1 ... N, while our values are on scale 0 ... N-1 and therefore after determining the linear interpolation weights, we subtract 1 from both lesserPos and upperPos; however ...
749  lesserWeight = upperPos - ( posVals[0] / maxSphereDists );
750  upperWeight = 1.0 - lesserWeight;
751 
752  posVals[1] = ( lesserWeight * vals[lesserPos-1][0] ) + ( upperWeight * vals[upperPos-1][0] );
753  posVals[2] = ( lesserWeight * vals[lesserPos-1][1] ) + ( upperWeight * vals[upperPos-1][1] );
754  }
755  else
756  {
757  //======================================== ... this then means that we would require position -1 for when the integration value is between 0 and the first shell. To resolve this, we assume that the values are 0 below the first shell and proceed as follows:
758  upperWeight = 1.0 - ( upperPos - ( posVals[0] / maxSphereDists ) );
759 
760  posVals[1] = ( upperWeight * vals[upperPos-1][0] );
761  posVals[2] = ( upperWeight * vals[upperPos-1][1] );
762  }
763 
764  intData[absIter][0] = posVals[0];
765  intData[absIter][1] = posVals[1];
766  intData[absIter][2] = posVals[2];
767  }
768 
769  //================================================ Integrate
770  *retReal = 0.0;
771  *retImag = 0.0;
772  for ( proshade_unsign absPoint = 0; absPoint < order; absPoint++ )
773  {
774  *retReal += ( weights[absPoint] * intData[absPoint][1] );
775  *retImag += ( weights[absPoint] * intData[absPoint][2] );
776  }
777 
778  //================================================ Normalise
779  *retReal *= ( integralOverRange / 2.0 );
780  *retImag *= ( integralOverRange / 2.0 );
781 
782  //================================================ Release memory
783  delete[] intData;
784 
785  //================================================ Done
786  return ;
787 
788 }

◆ gaussLegendreIntegrationReal()

proshade_double ProSHADE_internal_maths::gaussLegendreIntegrationReal ( proshade_double *  vals,
proshade_unsign  valsSize,
proshade_unsign  order,
proshade_double *  abscissas,
proshade_double *  weights,
proshade_double  integralOverRange,
proshade_double  maxSphereDists 
)

Function to compute real part of the Gauss-Legendre integration over spherical harmonic values in different shells.

This function takes the real parts of the spherical harmonics value in different shells and proceeds to compute the Gauss-Legendre integration over them. It uses the shell positions to appropriately place abscissas and their weights, which it assumes were pre-computed by the getLegendreAbscAndWeights() function.

Parameters
[in]valsPointer to an array of values over which the integration to be done.
[in]valsSizeThe length of the input array.
[in]orderThe integration order value.
[in]abscissasThe allocated array for holding the abscissa values.
[in]weightsThe allocated array for holding the weight values.
[in]integralOverRangeThe range of the intgral. If progressive shell mapping is used, this will not be max shell radius.
[in]maxSphereDistsDistance between two shells.
[out]XThe real part of Gauss-Legendre integration over the shperical harmonics values.

Definition at line 619 of file ProSHADE_maths.cpp.

620 {
621  //================================================ Initialise local variables
622  proshade_double ret = 0.0;
623  proshade_complex* intData = new proshade_complex[order];
624  ProSHADE_internal_misc::checkMemoryAllocation ( intData, __FILE__, __LINE__, __func__ );
625  proshade_complex posVals;
626  proshade_unsign lesserPos = 0;
627  proshade_unsign upperPos = 0;
628  proshade_double lesserWeight = 0.0;
629  proshade_double upperWeight = 0.0;
630 
631  //================================================ Rescale to <order> points
632  for ( proshade_unsign absIter = 0; absIter < order; absIter++ )
633  {
634  //============================================ Init loop
635  posVals[0] = 0.0;
636  posVals[1] = 0.0;
637 
638  //============================================ Find real position of abscissas
639  posVals[0] = ( ( abscissas[absIter] + 1.0 ) / 2.0 ) * integralOverRange;
640 
641 
642  //============================================ Find lesser and upper bounds
643  for ( proshade_unsign valIt = 0; valIt < valsSize; valIt++ )
644  {
645  if ( ( ( valIt * maxSphereDists ) <= posVals[0] ) && ( ( ( valIt + 1 ) * maxSphereDists ) > posVals[0] ) )
646  {
647  lesserPos = static_cast<proshade_unsign> ( valIt );
648  upperPos = static_cast<proshade_unsign> ( valIt + 1 );
649  break;
650  }
651  }
652 
653  //============================================ Linear Interpolation
654  lesserWeight = 0.0;
655  upperWeight = 0.0;
656  if ( lesserPos != 0 )
657  {
658  //======================================== Here we realise that the lesser and upper bounds were determined on scale 1 ... N, while our values are on scale 0 ... N-1 and therefore after determining the linear interpolation weights, we subtract 1 from both lesserPos and upperPos; however ...
659  lesserWeight = upperPos - ( posVals[0] / maxSphereDists );
660  upperWeight = 1.0 - lesserWeight;
661 
662  posVals[1] = ( lesserWeight * vals[lesserPos-1] ) + ( upperWeight * vals[upperPos-1] );
663  }
664  else
665  {
666  //======================================== ... this then means that we would require position -1 for when the integration value is between 0 and the first shell. To resolve this, we assume that the values are 0 below the first shell and proceed as follows:
667  upperWeight = 1.0 - ( upperPos - ( posVals[0] / maxSphereDists ) );
668 
669  posVals[1] = ( upperWeight * vals[upperPos-1] );
670  }
671 
672  intData[absIter][0] = posVals[0];
673  intData[absIter][1] = posVals[1];
674  }
675 
676  //================================================ Integrate
677  for ( proshade_unsign absPoint = 0; absPoint < order; absPoint++ )
678  {
679  ret += ( weights[absPoint] * intData[absPoint][1] );
680  }
681 
682  //================================================ Normalise
683  ret *= ( integralOverRange / 2.0 );
684 
685  //================================================ Release memory
686  delete[] intData;
687 
688  //================================================ Done
689  return ( ret );
690 
691 }

◆ getAxisAngleFromRotationMatrix() [1/2]

void ProSHADE_internal_maths::getAxisAngleFromRotationMatrix ( proshade_double *  rotMat,
proshade_double *  x,
proshade_double *  y,
proshade_double *  z,
proshade_double *  ang 
)

This function converts rotation matrix to the axis-angle representation.

This function takes a rotation matrix as an array of 9 numbers and converts it to the Angle-Axis representation, which is the main rotation representation used in ProSHADE. This function deals with both the North and South pole singularity of the rotation matrices.

Parameters
[in]rotMatRotation matrix as an array of 9 values.
[in]xPointer to which the x-axis value of the axis vector will be saved.
[in]yPointer to which the y-axis value of the axis vector will be saved.
[in]zPointer to which the z-axis value of the axis vector will be saved.
[in]angPointer to which the angle value will be saved.

Definition at line 1039 of file ProSHADE_maths.cpp.

1040 {
1041  //================================================ Initialise
1042  proshade_double singAtPiCheck = 0.01;
1043  proshade_double singAtIdentity = 0.05;
1044 
1045  //================================================ Check input for singularities
1046  if ( ( std::abs ( rotMat[1] - rotMat[3] ) < singAtPiCheck ) &&
1047  ( std::abs ( rotMat[2] - rotMat[6] ) < singAtPiCheck ) &&
1048  ( std::abs ( rotMat[5] - rotMat[7] ) < singAtPiCheck ) )
1049  {
1050  //============================================ Singularity in input! Check for identity matrix
1051  if ( ( std::abs ( rotMat[1] + rotMat[3] ) < singAtIdentity ) &&
1052  ( std::abs ( rotMat[2] + rotMat[6] ) < singAtIdentity ) &&
1053  ( std::abs ( rotMat[5] + rotMat[7] ) < singAtIdentity ) &&
1054  ( std::abs ( rotMat[0] + rotMat[4] + rotMat[8] - 3.0 ) < singAtIdentity ) )
1055  {
1056  //======================================== Identity matrix. Return 0 angle.
1057  *x = 1.0;
1058  *y = 0.0;
1059  *z = 0.0;
1060  *ang = 0.0;
1061 
1062  //======================================== Done
1063  return ;
1064  }
1065 
1066  //============================================ If we got here, this is the 180deg (pi rad) singularity. Find which axis should the rotation be done along
1067  *ang = M_PI;
1068 
1069  proshade_double xx = ( rotMat[0] + 1.0 ) / 2.0;
1070  proshade_double yy = ( rotMat[4] + 1.0 ) / 2.0;
1071  proshade_double zz = ( rotMat[8] + 1.0 ) / 2.0;
1072  proshade_double xy = ( rotMat[1] + rotMat[3] ) / 4.0;
1073  proshade_double xz = ( rotMat[2] + rotMat[6] ) / 4.0;
1074  proshade_double yz = ( rotMat[5] + rotMat[7] ) / 4.0;
1075 
1076  if ( ( xx > yy ) && ( xx > zz ) ) // XX is the largest diagonal
1077  {
1078  if ( xx < singAtPiCheck ) // and is still 0
1079  {
1080  *x = 0.0;
1081  *y = 1.0 / sqrt(2);
1082  *z = 1.0 / sqrt(2);
1083  }
1084  else
1085  {
1086  *x = sqrt ( xx );
1087  *y = xy / sqrt ( xx );
1088  *z = xz / sqrt ( xx );
1089  }
1090  }
1091 
1092  else if ( yy > zz ) // YY is the largest diagonal
1093  {
1094  if ( yy < singAtPiCheck ) // and is still 0
1095  {
1096  *x = 1.0 / sqrt(2);
1097  *y = 0.0;
1098  *z = 1.0 / sqrt(2);
1099  }
1100  else
1101  {
1102  *y = sqrt ( yy );
1103  *x = xy / sqrt ( yy );
1104  *z = yz / sqrt ( yy );
1105  }
1106  }
1107 
1108  else // ZZ is the largest diagonal
1109  {
1110  if ( zz < singAtPiCheck ) // and is still 0
1111  {
1112  *x = 1.0 / sqrt(2);
1113  *y = 1.0 / sqrt(2);
1114  *z = 0.0;
1115  }
1116  else
1117  {
1118  *z = sqrt ( zz );
1119  *x = xz / sqrt ( zz );
1120  *y = yz / sqrt ( zz );
1121  }
1122  }
1123 
1124  //============================================ Make sure largest axis is positive and so is the angle
1125  if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1126  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1127  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1128  {
1129  *x *= -1.0;
1130  *y *= -1.0;
1131  *z *= -1.0;
1132  *ang *= -1.0;
1133  }
1134  if ( *ang < 0.0 ) { *ang = ( 2.0 * M_PI ) + *ang; }
1135 
1136  //============================================ Done
1137  return ;
1138  }
1139 
1140  //================================================ No singularities! Now get angle
1141  *ang = std::acos ( ( std::max ( -1.0, std::min ( 3.0, rotMat[0] + rotMat[4] + rotMat[8] ) ) - 1.0 ) / 2.0 );
1142 
1143  //================================================ Init return values
1144  *x = 1.0;
1145  *y = 0.0;
1146  *z = 0.0;
1147 
1148  //================================================ Is angle 0? This should not happen, but will
1149  if ( std::abs ( *ang ) < singAtPiCheck )
1150  {
1151  *ang = 0.0;
1152  return ;
1153  }
1154 
1155  //================================================ Axis
1156  *x = rotMat[7] - rotMat[5];
1157  *y = rotMat[2] - rotMat[6];
1158  *z = rotMat[3] - rotMat[1];
1159 
1160  proshade_double normFactor = std::sqrt ( pow ( *x, 2.0 ) + pow ( *y, 2.0 ) + pow ( *z, 2.0 ) );
1161  *x /= normFactor;
1162  *y /= normFactor;
1163  *z /= normFactor;
1164 
1165  //================================================ Make sure largest axis is positive and so is the angle
1166  if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1167  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1168  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1169  {
1170  *x *= -1.0;
1171  *y *= -1.0;
1172  *z *= -1.0;
1173  *ang *= -1.0;
1174  }
1175  if ( *ang < 0.0 ) { *ang = ( 2.0 * M_PI ) + *ang; }
1176 
1177  //================================================ Done
1178  return ;
1179 
1180 }

◆ getAxisAngleFromRotationMatrix() [2/2]

void ProSHADE_internal_maths::getAxisAngleFromRotationMatrix ( std::vector< proshade_double > *  rotMat,
proshade_double *  x,
proshade_double *  y,
proshade_double *  z,
proshade_double *  ang 
)

This function converts rotation matrix to the axis-angle representation.

This function takes a rotation matrix as a pointer to a vector of doubles and converts it to the Angle-Axis representation, which is the main rotation representation used in ProSHADE. This function deals with both the North and South pole singularity of the rotation matrices.

Parameters
[in]rotMatRotation matrix as a pointer to a vector of doubles.
[in]xPointer to which the x-axis value of the axis vector will be saved.
[in]yPointer to which the y-axis value of the axis vector will be saved.
[in]zPointer to which the z-axis value of the axis vector will be saved.
[in]angPointer to which the angle value will be saved.

Definition at line 1194 of file ProSHADE_maths.cpp.

1195 {
1196  //================================================ Initialise
1197  proshade_double singAtPiCheck = 0.01;
1198  proshade_double singAtIdentity = 0.05;
1199 
1200  //================================================ Check input for singularities
1201  if ( ( std::abs ( rotMat->at(1) - rotMat->at(3) ) < singAtPiCheck ) &&
1202  ( std::abs ( rotMat->at(2) - rotMat->at(6) ) < singAtPiCheck ) &&
1203  ( std::abs ( rotMat->at(5) - rotMat->at(7) ) < singAtPiCheck ) )
1204  {
1205  //============================================ Singularity in input! Check for identity matrix
1206  if ( ( std::abs ( rotMat->at(1) + rotMat->at(3) ) < singAtIdentity ) &&
1207  ( std::abs ( rotMat->at(2) + rotMat->at(6) ) < singAtIdentity ) &&
1208  ( std::abs ( rotMat->at(5) + rotMat->at(7) ) < singAtIdentity ) &&
1209  ( std::abs ( rotMat->at(0) + rotMat->at(4) + rotMat->at(8) - 3.0 ) < singAtIdentity ) )
1210  {
1211  //======================================== Identity matrix. Return 0 angle.
1212  *x = 1.0;
1213  *y = 0.0;
1214  *z = 0.0;
1215  *ang = 0.0;
1216 
1217  //======================================== Done
1218  return ;
1219  }
1220 
1221  //============================================ If we got here, this is the 180deg (pi rad) singularity. Find which axis should the rotation be done along
1222  *ang = M_PI;
1223 
1224  proshade_double xx = ( rotMat->at(0) + 1.0 ) / 2.0;
1225  proshade_double yy = ( rotMat->at(4) + 1.0 ) / 2.0;
1226  proshade_double zz = ( rotMat->at(8) + 1.0 ) / 2.0;
1227  proshade_double xy = ( rotMat->at(1) + rotMat->at(3) ) / 4.0;
1228  proshade_double xz = ( rotMat->at(2) + rotMat->at(6) ) / 4.0;
1229  proshade_double yz = ( rotMat->at(5) + rotMat->at(7) ) / 4.0;
1230 
1231  if ( ( xx > yy ) && ( xx > zz ) ) // XX is the largest diagonal
1232  {
1233  if ( xx < singAtPiCheck ) // and is still 0
1234  {
1235  *x = 0.0;
1236  *y = 1.0 / sqrt(2);
1237  *z = 1.0 / sqrt(2);
1238  }
1239  else
1240  {
1241  *x = sqrt ( xx );
1242  *y = xy / sqrt ( xx );
1243  *z = xz / sqrt ( xx );
1244  }
1245  }
1246 
1247  else if ( yy > zz ) // YY is the largest diagonal
1248  {
1249  if ( yy < singAtPiCheck ) // and is still 0
1250  {
1251  *x = 1.0 / sqrt(2);
1252  *y = 0.0;
1253  *z = 1.0 / sqrt(2);
1254  }
1255  else
1256  {
1257  *y = sqrt ( yy );
1258  *x = xy / sqrt ( yy );
1259  *z = yz / sqrt ( yy );
1260  }
1261  }
1262 
1263  else // ZZ is the largest diagonal
1264  {
1265  if ( zz < singAtPiCheck ) // and is still 0
1266  {
1267  *x = 1.0 / sqrt(2);
1268  *y = 1.0 / sqrt(2);
1269  *z = 0.0;
1270  }
1271  else
1272  {
1273  *z = sqrt ( zz );
1274  *x = xz / sqrt ( zz );
1275  *y = yz / sqrt ( zz );
1276  }
1277  }
1278 
1279  //============================================ Make sure largest axis is positive and so is the angle
1280  if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1281  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1282  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1283  {
1284  *x *= -1.0;
1285  *y *= -1.0;
1286  *z *= -1.0;
1287  *ang *= -1.0;
1288  }
1289 
1290  //============================================ Done
1291  return ;
1292  }
1293 
1294  //================================================ No singularities! Now get angle
1295  *ang = std::acos ( ( std::max ( -1.0, std::min ( 3.0, rotMat->at(0) + rotMat->at(4) + rotMat->at(8) ) ) - 1.0 ) / 2.0 );
1296 
1297  //================================================ Init return values
1298  *x = 1.0;
1299  *y = 0.0;
1300  *z = 0.0;
1301 
1302  //================================================ Is angle 0? This should not happen, but will
1303  if ( std::abs ( *ang ) < singAtPiCheck )
1304  {
1305  *ang = 0.0;
1306  return ;
1307  }
1308 
1309  //================================================ Axis
1310  *x = rotMat->at(7) - rotMat->at(5);
1311  *y = rotMat->at(2) - rotMat->at(6);
1312  *z = rotMat->at(3) - rotMat->at(1);
1313 
1314  proshade_double normFactor = std::sqrt ( pow ( *x, 2.0 ) + pow ( *y, 2.0 ) + pow ( *z, 2.0 ) );
1315  *x /= normFactor;
1316  *y /= normFactor;
1317  *z /= normFactor;
1318 
1319  //================================================ Make sure largest axis is positive and so is the angle
1320  if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1321  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1322  ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1323  {
1324  *x *= -1.0;
1325  *y *= -1.0;
1326  *z *= -1.0;
1327  *ang *= -1.0;
1328  }
1329 
1330  //================================================ Done
1331  return ;
1332 
1333 }

◆ getEulerZXZFromAngleAxis()

void ProSHADE_internal_maths::getEulerZXZFromAngleAxis ( proshade_double  axX,
proshade_double  axY,
proshade_double  axZ,
proshade_double  axAng,
proshade_double *  eA,
proshade_double *  eB,
proshade_double *  eG,
proshade_unsign  angDim 
)

This function converts angle-axis representation to the Euler ZXZ angles representation.

This function does the angle-axis to Euler ZXZ conversion and if a problem around the Z axis arises, it deal with it.

Parameters
[in]axXAngle-axis representation axis x element.
[in]axYAngle-axis representation axis y element.
[in]axZAngle-axis representation axis z element.
[in]axAngAngle-axis representation angle.
[in]eAPointer to which the Euler angle alpha value will be saved.
[in]eBPointer to which the Euler angle beta value will be saved.
[in]eGPointer to which the Euler angle gamma value will be saved.

Definition at line 1437 of file ProSHADE_maths.cpp.

1438 {
1439  //================================================ If angle is 0 or infinity (anything divided by 0), return no rotation
1440  if ( ( axAng == 0.0 ) || ( std::isinf ( axAng ) ) )
1441  {
1442  //============================================ Return 0 ; 0 ; 0 for no angle
1443  *eA = 0.0;
1444  *eB = 0.0;
1445  *eG = 0.0;
1446 
1447  //============================================ Done
1448  return ;
1449  }
1450 
1451  //================================================ Compute required rotation matrix elements
1452  proshade_double cAng = std::cos ( axAng );
1453  proshade_double sAng = std::sin ( axAng );
1454  proshade_double tAng = 1.0 - cAng;
1455 
1456  proshade_double element22 = cAng + axZ * axZ * tAng;
1457 
1458  proshade_double tmp1 = axX * axZ * tAng;
1459  proshade_double tmp2 = axY * sAng;
1460  proshade_double element20 = tmp1 - tmp2;
1461  proshade_double element02 = tmp1 + tmp2;
1462 
1463  tmp1 = axY * axZ * tAng;
1464  tmp2 = axX * sAng;
1465  proshade_double element21 = tmp1 + tmp2;
1466  proshade_double element12 = tmp1 - tmp2;
1467 
1468  //================================================ Convert to Eulers
1469  if ( std::abs( element22 ) <= 0.99999 )
1470  {
1471  //============================================ This case occurs when there is no singularity in the rotation matrix (i.e. it does not have 0 or 180 degrees angle)
1472  *eA = std::atan2 ( element21, element20 );
1473  *eB = std::acos ( element22 );
1474  *eG = std::atan2 ( element12, -element02 );
1475  }
1476  else
1477  {
1478  //============================================ Compute some extra rotation matrix elements
1479  tmp1 = axX * axY * tAng;
1480  tmp2 = axZ * sAng;
1481  proshade_double element10 = tmp1 + tmp2;
1482  proshade_double element00 = cAng + axX * axX * tAng;
1483 
1484  //============================================ This case occurs when there is either 0 or 180 degrees rotation angle in the rotation matrix and therefore when beta is zero.
1485  if ( element22 >= 0.99999 )
1486  {
1487  //======================================== In this case, beta = 0 and alpha and gamma are only defined in terms of their sum. So we arbitrarily set gamma to 0 and solve alpha.
1488  *eA = std::atan2 ( element10, element00 );
1489  *eB = 0.0;
1490  *eG = 0.0;
1491  }
1492  if ( element22 <= -0.99999 )
1493  {
1494  //======================================== In this case, beta = 0 and alpha and gamma are only defined in terms of their difference. So we arbitrarily set gamma to 0 and solve alpha.
1495  *eA = std::atan2 ( element10, element00 );
1496  *eB = M_PI / 2.0;
1497  *eG = 0.0;
1498  }
1499  }
1500 
1501  //================================================ Get the angles to proper range
1502  if ( *eA < 0.0 ) { *eA = 2.0 * M_PI + *eA; }
1503  if ( *eB < 0.0 ) { *eB = M_PI + *eB; }
1504  if ( *eG < 0.0 ) { *eG = 2.0 * M_PI + *eG; }
1505 
1506  //================================================ Done
1507  return ;
1508 
1509 }

◆ getEulerZXZFromAngleAxisFullSearch()

void ProSHADE_internal_maths::getEulerZXZFromAngleAxisFullSearch ( proshade_double  axX,
proshade_double  axY,
proshade_double  axZ,
proshade_double  axAng,
proshade_double *  eA,
proshade_double *  eB,
proshade_double *  eG,
proshade_signed  angDim 
)

This function converts angle-axis representation to the Euler ZXZ angles representation using full search.

This function is meant for solving the issue of angle-axis conversion to Euler ZXZ convention for axis 0,0,1, where all the rotation matrix elements used for Euler alpha and gamma angles are 0.0. The function overcomes this by simply searching all the rotation function indices for having angle-axis value similar to the required one - a rather slow approach. Therefore, the getEulerZXZFromAngleAxis() function should be used instead and only if it fails (has all angles 0.0), then this function should be used instead.

Parameters
[in]axXAngle-axis representation axis x element.
[in]axYAngle-axis representation axis y element.
[in]axZAngle-axis representation axis z element.
[in]axAngAngle-axis representation angle.
[in]eAPointer to which the Euler angle alpha value will be saved.
[in]eBPointer to which the Euler angle beta value will be saved.
[in]eGPointer to which the Euler angle gamma value will be saved.

Definition at line 1526 of file ProSHADE_maths.cpp.

1527 {
1528  //================================================ Initialise variables
1529  proshade_double bestDist = 999.9;
1530  proshade_double eAHlp, eBHlp, eGHlp, axXHlp, axYHlp, axZHlp, axAngHlp, axDist;
1531 
1532  //================================================ Allocate memory
1533  proshade_double* rMat = new proshade_double[9];
1534  ProSHADE_internal_misc::checkMemoryAllocation ( rMat, __FILE__, __LINE__, __func__ );
1535 
1536  //================================================ For each rotation function index (i.e. existing Euler angles ZXZ combination)
1537  for ( proshade_signed xIt = 0; xIt < angDim; xIt++ )
1538  {
1539  for ( proshade_signed yIt = 0; yIt < angDim; yIt++ )
1540  {
1541  for ( proshade_signed zIt = 0; zIt < angDim; zIt++ )
1542  {
1543  //==================================== Speed up
1544  if ( bestDist < 0.001 ) { break; }
1545 
1546  //==================================== Get Euler ZXZ from the indices
1547  getEulerZXZFromSOFTPosition ( angDim/2, xIt, yIt, zIt, &eAHlp, &eBHlp, &eGHlp );
1548  getRotationMatrixFromEulerZXZAngles ( eAHlp, eBHlp, eGHlp, rMat );
1549  getAxisAngleFromRotationMatrix ( rMat, &axXHlp, &axYHlp, &axZHlp, &axAngHlp );
1550 
1551  //==================================== If angle is larger than 180 degrees
1552  if ( axAng > M_PI )
1553  {
1554  axAng = ( 2.0 * M_PI ) - axAng;
1555  axAng *= -1.0;
1556  }
1557 
1558  //==================================== Make sure vector direction is the same
1559  if ( ( ( std::max( std::abs( axXHlp ), std::max( std::abs( axYHlp ), std::abs( axZHlp ) ) ) == std::abs( axXHlp ) ) && ( axXHlp < 0.0 ) ) ||
1560  ( ( std::max( std::abs( axXHlp ), std::max( std::abs( axYHlp ), std::abs( axZHlp ) ) ) == std::abs( axYHlp ) ) && ( axYHlp < 0.0 ) ) ||
1561  ( ( std::max( std::abs( axXHlp ), std::max( std::abs( axYHlp ), std::abs( axZHlp ) ) ) == std::abs( axZHlp ) ) && ( axZHlp < 0.0 ) ) )
1562  {
1563  axXHlp *= -1.0;
1564  axYHlp *= -1.0;
1565  axZHlp *= -1.0;
1566  axAngHlp *= -1.0;
1567  }
1568 
1569  if ( ( ( std::max( std::abs( axX ), std::max( std::abs( axY ), std::abs( axZ ) ) ) == std::abs( axX ) ) && ( axX < 0.0 ) ) ||
1570  ( ( std::max( std::abs( axX ), std::max( std::abs( axY ), std::abs( axZ ) ) ) == std::abs( axY ) ) && ( axY < 0.0 ) ) ||
1571  ( ( std::max( std::abs( axX ), std::max( std::abs( axY ), std::abs( axZ ) ) ) == std::abs( axZ ) ) && ( axZ < 0.0 ) ) )
1572  {
1573  axX *= -1.0;
1574  axY *= -1.0;
1575  axZ *= -1.0;
1576  axAng *= -1.0;
1577  }
1578 
1579  //==================================== Compute distance to the requested angle-axis values
1580  axDist = std::abs( axAng - axAngHlp ) + ( 1.0 - std::abs ( ( ( axX * axXHlp ) + ( axY * axYHlp ) + ( axZ * axZHlp ) ) /
1581  ( sqrt( pow( axX, 2.0 ) + pow( axY, 2.0 ) + pow( axZ, 2.0 ) ) * sqrt( pow( axXHlp, 2.0 ) + pow( axYHlp, 2.0 ) + pow( axZHlp, 2.0 ) ) ) ) );
1582 
1583  //==================================== Is this point an improvement
1584  if ( std::abs ( axDist ) < bestDist )
1585  {
1586  //================================ If so, note it
1587  bestDist = std::abs ( axDist );
1588  *eA = eAHlp;
1589  *eB = eBHlp;
1590  *eG = eGHlp;
1591  }
1592  }
1593  }
1594  }
1595 
1596  //================================================ Release memory
1597  delete[] rMat;
1598 
1599  //================================================ Done
1600  return ;
1601 
1602 }

◆ getEulerZXZFromRotMatrix()

void ProSHADE_internal_maths::getEulerZXZFromRotMatrix ( proshade_double *  rotMat,
proshade_double *  eA,
proshade_double *  eB,
proshade_double *  eG 
)

This function converts rotation matrix to the Euler ZXZ angles representation.

Parameters
[in]rotMatRotation matrix as an array of 9 values.
[in]eAPointer to which the Euler angle alpha value will be saved.
[in]eBPointer to which the Euler angle beta value will be saved.
[in]eGPointer to which the Euler angle gamma value will be saved.

Definition at line 1394 of file ProSHADE_maths.cpp.

1395 {
1396  //================================================ Get ZXZ Euler from matrix
1397  *eA = atan2 ( rotMat[7], rotMat[6] );
1398  *eB = acos ( rotMat[8] );
1399  *eG = atan2 ( rotMat[5], -rotMat[2] );
1400 
1401  //================================================ Solve undefined 0,0 inputs (i.e. identity matrix)
1402  proshade_double errLimit = 0.001;
1403  if ( ( ( rotMat[7] < errLimit ) && ( rotMat[7] > -errLimit ) ) && ( ( rotMat[6] < errLimit ) && ( rotMat[6] > -errLimit ) ) )
1404  {
1405  //============================================ atan2 (0,0) is undefined, we want 0.0 here
1406  *eA = 0.0;
1407  }
1408 
1409  if ( ( ( rotMat[5] < errLimit ) && ( rotMat[5] > -errLimit ) ) && ( ( rotMat[2] < errLimit ) && ( rotMat[2] > -errLimit ) ) )
1410  {
1411  //============================================ atan2 (0,0) is undefined, we want 0.0 here
1412  *eG = 0.0;
1413  }
1414 
1415  //================================================ Get the angles to proper range
1416  if ( *eA < 0.0 ) { *eA = 2.0 * M_PI + *eA; }
1417  if ( *eB < 0.0 ) { *eB = M_PI + *eB; }
1418  if ( *eG < 0.0 ) { *eG = 2.0 * M_PI + *eG; }
1419 
1420  //================================================ Done
1421  return ;
1422 
1423 }

◆ getEulerZXZFromSOFTPosition()

void ProSHADE_internal_maths::getEulerZXZFromSOFTPosition ( proshade_signed  band,
proshade_signed  x,
proshade_signed  y,
proshade_signed  z,
proshade_double *  eulerAlpha,
proshade_double *  eulerBeta,
proshade_double *  eulerGamma 
)

Function to find Euler angles (ZXZ convention) from index position in the inverse SOFT map.

This function proceeds to convert the inverse SOFT map x, y and z position to Euler ZXZ convention angles, saving these into the supplied pointers.

Parameters
[in]bandThe maximum bandwidth of the computation.
[in]xThe x-axis position in the inverse SOFT map.
[in]yThe y-axis position in the inverse SOFT map.
[in]zThe z-axis position in the inverse SOFT map.
[in]eulerAlphaPointer to where the Euler alpha angle will be saved.
[in]eulerBetaPointer to where the Euler beta angle will be saved.
[in]eulerGammaPointer to where the Euler gamma angle will be saved.

Definition at line 961 of file ProSHADE_maths.cpp.

962 {
963  //================================================ Convert index to Euler angles
964  *eulerGamma = ( M_PI * y / static_cast<proshade_double> ( 1.0 * band ) );
965  *eulerBeta = ( M_PI * x / static_cast<proshade_double> ( 2.0 * band ) );
966  *eulerAlpha = ( M_PI * z / static_cast<proshade_double> ( 1.0 * band ) );
967 
968  //================================================ Done
969  return ;
970 
971 }

◆ getGLFirstEvenRoot()

void ProSHADE_internal_maths::getGLFirstEvenRoot ( proshade_double  polyAtZero,
proshade_unsign  order,
proshade_double *  abscAtZero,
proshade_double *  weighAtZero,
proshade_unsign  taylorSeriesCap 
)

This function finds the first root for Legendre polynomials of odd order.

The Legendre polynomials with odd order have zero as the first root, but the even oder polenomials have different value and this function serves the purpose of finding this value (i.e. the first root of the polynomial if the order is even).

Parameters
[in]polyAtZeroThe value of the polynomial at zero.
[in]orderThe positive integer value of the polynomial order.
[in]abscAtZeroPointer to variable storing the abscissa value at zero.
[in]weightAtZeroPointer to variable storing the weight value at zero.
[in]taylorSeriesCapThe limit on the Taylor series.

Definition at line 384 of file ProSHADE_maths.cpp.

385 {
386  //================================================ Sanity check
387  if ( taylorSeriesCap < 2 )
388  {
389  throw ProSHADE_exception ( "The Taylor series cap is too low.", "EI00020", __FILE__, __LINE__, __func__, "The Taylor series expansion limit is less than 2. This\n : seems very low; if you have a very small structure or very\n : low resolution, please manually increase the integration\n : order. Otherwise, please report this as a bug." );
390  }
391 
392  //================================================ Initialise variables
393  *abscAtZero = advanceGLPolyValue ( 0.0, -M_PI / 2.0, 0.0, order, taylorSeriesCap );
394  proshade_double hlp = 0.0;
395  proshade_double hlpVal = static_cast<proshade_double> ( order );
396  proshade_double *abscSteps;
397  proshade_double *weightSteps;
398 
399  //================================================ Allocate memory
400  abscSteps = new proshade_double [taylorSeriesCap+2];
401  weightSteps = new proshade_double [taylorSeriesCap+1];
402 
403  //================================================ Pre-set values
404  abscSteps[0] = 0.0;
405  abscSteps[1] = polyAtZero;
406  weightSteps[0] = 0.0;
407 
408  //================================================ Fill in abscissa and weight steps
409  for ( proshade_unsign iter = 0; iter <= taylorSeriesCap - 2; iter = iter + 2 )
410  {
411  hlp = static_cast<proshade_double> ( iter );
412 
413  abscSteps[iter+2] = 0.0;
414  abscSteps[iter+3] = ( hlp * ( hlp + 1.0 ) - hlpVal * ( hlpVal + 1.0 ) ) * abscSteps[iter+1] / (hlp + 1.0) / (hlp + 2.0 );
415 
416  weightSteps[iter+1] = 0.0;
417  weightSteps[iter+2] = ( hlp + 2.0 ) * abscSteps[iter+3];
418  }
419 
420  //================================================ Find abscissa and weights
421  for ( proshade_double iter = 0; iter < 5; iter++ )
422  {
423  *abscAtZero = *abscAtZero - evaluateGLSeries ( abscSteps, *abscAtZero, taylorSeriesCap ) / evaluateGLSeries ( weightSteps, *abscAtZero, taylorSeriesCap-1 );
424  }
425  *weighAtZero = evaluateGLSeries ( weightSteps, *abscAtZero, taylorSeriesCap-1 );
426 
427  //================================================ Free memory
428  delete abscSteps;
429  delete weightSteps;
430 
431  //================================================ Done
432  return ;
433 
434 }

◆ getGLPolyAtZero()

void ProSHADE_internal_maths::getGLPolyAtZero ( proshade_unsign  order,
proshade_double *  polyValue,
proshade_double *  deriValue 
)

This function obtains the Legendre polynomial values and its derivative at zero for any positive integer order polynomial.

This function takes the positive integer order of the Legendre polynomial and uses the recursive properties of the polynomials to work up to the order, computing the value at zero and its derivative for all lesser orders. It then returns the final values.

Parameters
[in]orderPositive integer order of the Legendre polynomial which value at zero we want.
[in]polyValuePointer to variable which will store the resulting polynomial value at zero.
[in]deriValuePointer to variable which will store the derivative of the zero value.

Definition at line 347 of file ProSHADE_maths.cpp.

348 {
349  //================================================ Initialise
350  proshade_double hlpVal = 0.0;
351  proshade_double prevPoly = 1.0;
352  proshade_double prevPrevPoly = 0.0;
353  proshade_double prevDeri = 0.0;
354  proshade_double prevPrevDeri = 0.0;
355 
356  for ( proshade_unsign ordIt = 0; ordIt < order; ordIt++ )
357  {
358  hlpVal = static_cast<proshade_double> ( ordIt );
359  *polyValue = -hlpVal * prevPrevPoly / ( hlpVal + 1.0 );
360  *deriValue = ( ( 2.0 * hlpVal + 1.0 ) * prevPoly - hlpVal * prevPrevDeri ) / ( hlpVal + 1.0 );
361  prevPrevPoly = prevPoly;
362  prevPoly = *polyValue;
363  prevPrevDeri = prevDeri;
364  prevDeri = *deriValue;
365  }
366 
367  //================================================ Done
368  return ;
369 
370 }

◆ getLegendreAbscAndWeights()

void ProSHADE_internal_maths::getLegendreAbscAndWeights ( proshade_unsign  order,
proshade_double *  abscissas,
proshade_double *  weights,
proshade_unsign  taylorSeriesCap 
)

Function to prepare abscissas and weights for Gauss-Legendre integration.

This function fills in the Gauss-Legendre interpolation points positions (abscissas) and their weights vectors, which will then be used for computing the Gauss-Legendre interpolation.

Parameters
[in]orderThe order to which the abscissas and weights should be prepared.
[in]abscissasThe array holding the abscissa values.
[in]weightsThe array holding the weight values.
[in]taylorSeriesCapThe limit on the Taylor series.

Definition at line 287 of file ProSHADE_maths.cpp.

288 {
289  //================================================ Sanity check
290  if ( order < 2 )
291  {
292  throw ProSHADE_exception ( "The integration order is too low.", "EI00019", __FILE__, __LINE__, __func__, "The Gauss-Legendre integration order is less than 2. This\n : seems very low; if you have a very small structure or very\n : low resolution, please manually increase the integration\n : order. Otherwise, please report this as a bug." );
293  }
294 
295  //================================================ Initialise
296  proshade_double polyValue = 0.0;
297  proshade_double deriValue = 0.0;
298  proshade_double weightSum = 0.0;
299 
300  //================================================ Find the polynomial and derivative values at 0
301  getGLPolyAtZero ( order,
302  &polyValue,
303  &deriValue );
304 
305  //================================================ If the order is odd, then 0 is a root ...
306  if ( order % 2 == 1 )
307  {
308  abscissas[((order-1)/2)] = polyValue;
309  weights[((order-1)/2)] = deriValue;
310  }
311  else
312  {
313  // ... and if order is even, find the first root
314  getGLFirstEvenRoot ( polyValue, order, &abscissas[(order/2)], &weights[(order/2)], taylorSeriesCap );
315  }
316 
317  //================================================ Now, having computed the first roots, complete the series
318  completeLegendreSeries ( order, abscissas, weights, taylorSeriesCap );
319 
320  //================================================ Correct weights by anscissa values
321  for ( proshade_unsign iter = 0; iter < order; iter++ )
322  {
323  weights[iter] = 2.0 / ( 1.0 - abscissas[iter] ) / ( 1.0 + abscissas[iter] ) / weights[iter] / weights[iter];
324  weightSum = weightSum + weights[iter];
325  }
326 
327  //================================================ Normalise weights
328  for ( proshade_unsign iter = 0; iter < order; iter++ )
329  {
330  weights[iter] = 2.0 * weights[iter] / weightSum;
331  }
332 
333  //================================================ Done
334  return ;
335 }

◆ getRotationMatrixFromAngleAxis()

void ProSHADE_internal_maths::getRotationMatrixFromAngleAxis ( proshade_double *  rotMat,
proshade_double  x,
proshade_double  y,
proshade_double  z,
proshade_double  ang 
)

This function converts the axis-angle representation to the rotation matrix representation.

Parameters
[in]rotMatRotation matrix as an array of 9 values will be saved to this pointer, must already be allocated.
[in]xThe x-axis value of the axis vector.
[in]yThe y-axis value of the axis vector.
[in]zThe z-axis value of the axis vector.
[in]angTheangle value.

Definition at line 1343 of file ProSHADE_maths.cpp.

1344 {
1345  //================================================ If angle is 0 or infinity (anything divided by 0), return identity matrix
1346  if ( ( ang == 0.0 ) || ( std::isinf ( ang ) ) )
1347  {
1348  //============================================ Create identity
1349  for ( proshade_unsign i = 0; i < 9; i++ ) { rotMat[i] = 0.0; }
1350  rotMat[0] = 1.0;
1351  rotMat[4] = 1.0;
1352  rotMat[8] = 1.0;
1353 
1354  //============================================ Done
1355  return ;
1356  }
1357 
1358  //================================================ Compute the matrix
1359  proshade_double cAng = cos ( ang );
1360  proshade_double sAng = sin ( ang );
1361  proshade_double tAng = 1.0 - cAng;
1362 
1363  rotMat[0] = cAng + x * x * tAng;
1364  rotMat[4] = cAng + y * y * tAng;
1365  rotMat[8] = cAng + z * z * tAng;
1366 
1367  proshade_double tmp1 = x * y * tAng;
1368  proshade_double tmp2 = z * sAng;
1369  rotMat[3] = tmp1 + tmp2;
1370  rotMat[1] = tmp1 - tmp2;
1371 
1372  tmp1 = x * z * tAng;
1373  tmp2 = y * sAng;
1374  rotMat[6] = tmp1 - tmp2;
1375  rotMat[2] = tmp1 + tmp2;
1376 
1377  tmp1 = y * z * tAng;
1378  tmp2 = x * sAng;
1379  rotMat[7] = tmp1 + tmp2;
1380  rotMat[5] = tmp1 - tmp2;
1381 
1382  //================================================ Done
1383  return ;
1384 
1385 }

◆ getRotationMatrixFromEulerZXZAngles()

void ProSHADE_internal_maths::getRotationMatrixFromEulerZXZAngles ( proshade_double  eulerAlpha,
proshade_double  eulerBeta,
proshade_double  eulerGamma,
proshade_double *  matrix 
)

Function to find the rotation matrix from Euler angles (ZXZ convention).

Parameters
[in]eulerAlphaThe Euler alpha angle value.
[in]eulerBetaThe Euler beta angle value.
[in]eulerGammaThe Euler gamma angle value.
[in]matrixA pointer to array of 9 values to which the results of the function will be saved.

Definition at line 1005 of file ProSHADE_maths.cpp.

1006 {
1007  //================================================ First row
1008  matrix[0] = cos ( eulerAlpha ) * cos ( eulerBeta ) * cos ( eulerGamma ) - sin ( eulerAlpha ) * sin ( eulerGamma );
1009  matrix[1] = sin ( eulerAlpha ) * cos ( eulerBeta ) * cos ( eulerGamma ) + cos ( eulerAlpha ) * sin ( eulerGamma );
1010  matrix[2] = -sin ( eulerBeta ) * cos ( eulerGamma );
1011 
1012  //================================================ Second row
1013  matrix[3] = -cos ( eulerAlpha ) * cos ( eulerBeta ) * sin ( eulerGamma ) - sin ( eulerAlpha ) * cos ( eulerGamma );
1014  matrix[4] = -sin ( eulerAlpha ) * cos ( eulerBeta ) * sin ( eulerGamma ) + cos ( eulerAlpha ) * cos ( eulerGamma );
1015  matrix[5] = sin ( eulerBeta ) * sin ( eulerGamma );
1016 
1017  //================================================ Third row
1018  matrix[6] = cos ( eulerAlpha ) * sin ( eulerBeta );
1019  matrix[7] = sin ( eulerAlpha ) * sin ( eulerBeta );
1020  matrix[8] = cos ( eulerBeta );
1021 
1022  //================================================ Done
1023  return ;
1024 
1025 }

◆ getSOFTPositionFromEulerZXZ()

void ProSHADE_internal_maths::getSOFTPositionFromEulerZXZ ( proshade_signed  band,
proshade_double  eulerAlpha,
proshade_double  eulerBeta,
proshade_double  eulerGamma,
proshade_double *  x,
proshade_double *  y,
proshade_double *  z 
)

Function to find the index position in the inverse SOFT map from given Euler angles (ZXZ convention).

This function does the conversion from Euler angles ZXZ convention to the SOFT map x, y and z position. It is not limitted to the SOFT map indices and instead if given Euler agnles between two indices will return a decimal point for the indices.

Parameters
[in]bandThe maximum bandwidth of the computation.
[in]eulerAlphaThe Euler alpha angle value.
[in]eulerBetaThe Euler beta angle value.
[in]eulerGammaThe Euler gamma angle value.
[in]xPointer to where the closest x-axis position in the inverse SOFT map will be saved to (position may be decimal!).
[in]yPointer to where the closest y-axis position in the inverse SOFT map will be saved to (position may be decimal!).
[in]zPointer to where the closest z-axis position in the inverse SOFT map will be saved to (position may be decimal!).

Definition at line 986 of file ProSHADE_maths.cpp.

987 {
988  //================================================ Convert Euler angles to indices
989  *x = static_cast<proshade_double> ( ( eulerBeta * static_cast<proshade_double> ( 2.0 * band ) ) / M_PI );
990  *y = static_cast<proshade_double> ( ( eulerGamma * static_cast<proshade_double> ( band ) ) / M_PI );
991  *z = static_cast<proshade_double> ( ( eulerAlpha * static_cast<proshade_double> ( band ) ) / M_PI );
992 
993  //================================================ Done
994  return ;
995 
996 }

◆ isAxisUnique() [1/2]

bool ProSHADE_internal_maths::isAxisUnique ( std::vector< proshade_double * > *  CSymList,
proshade_double *  axis,
proshade_double  tolerance = 0.1,
bool  improve = false 
)

This function checks if new axis is unique, or already detected.

This function compares the supplied axis against all members of the axes vector. If the axis has the same fold and very similar axis vector (i.e. all three elements are within tolerance), then the function returns false. If no such match is found, true is returned.

Parameters
[in]CSymListA vector containing the already detected Cyclic symmetries.
[in]axisThe axis to be checked against CSymList to see if it not already present.
[in]toleranceThe allowed error on each dimension of the axis.
[in]improveIf a similar axis is found and if this already existing axis has lower peak height, should the CSymList be updated with the higher peak height axis?
[out]retBoolean specifying whether a similar axis was found or not.

Definition at line 2613 of file ProSHADE_maths.cpp.

2614 {
2615  //================================================ Initialise variables
2616  bool ret = true;
2617  proshade_unsign whichImprove;
2618 
2619  //================================================ For each already detected member
2620  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( CSymList->size() ); grIt++ )
2621  {
2622  //============================================ Is fold the same?
2623  if ( CSymList->at(grIt)[0] == axis[0] )
2624  {
2625  if ( ProSHADE_internal_maths::vectorOrientationSimilarity ( CSymList->at(grIt)[1], CSymList->at(grIt)[2], CSymList->at(grIt)[3], axis[1], axis[2], axis[3], tolerance ) )
2626  {
2627  ret = false;
2628  whichImprove = grIt;
2629  break;
2630  }
2631  }
2632  }
2633 
2634  //================================================ Improve, if required
2635  if ( improve && !ret )
2636  {
2637  CSymList->at(whichImprove)[1] = axis[1];
2638  CSymList->at(whichImprove)[2] = axis[2];
2639  CSymList->at(whichImprove)[3] = axis[3];
2640  CSymList->at(whichImprove)[4] = axis[4];
2641  CSymList->at(whichImprove)[5] = axis[5];
2642  }
2643 
2644  //================================================ Done
2645  return ( ret );
2646 
2647 }

◆ isAxisUnique() [2/2]

bool ProSHADE_internal_maths::isAxisUnique ( std::vector< proshade_double * > *  CSymList,
proshade_double  X,
proshade_double  Y,
proshade_double  Z,
proshade_double  fold,
proshade_double  tolerance 
)

This function checks if new axis is unique, or already detected.

This function compares the supplied axis against all members of the axes vector. If the axis has the same fold and very similar axis vector (i.e. all three elements are within tolerance), then the function returns false. If no such match is found, true is returned.

Parameters
[in]CSymListA vector containing the already detected Cyclic symmetries.
[in]XThe axis x-element to be checked against CSymList to see if it not already present.
[in]YThe axis x-element to be checked against CSymList to see if it not already present.
[in]ZThe axis x-element to be checked against CSymList to see if it not already present.
[in]toleranceThe allowed error on each dimension of the axis.
[out]retBoolean specifying whether a similar axis was found or not.

Definition at line 2661 of file ProSHADE_maths.cpp.

2662 {
2663  //================================================ Initialise variables
2664  bool ret = true;
2665 
2666  //================================================ For each already detected member
2667  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( CSymList->size() ); grIt++ )
2668  {
2669  if ( fold == CSymList->at(grIt)[0] )
2670  {
2671  if ( ProSHADE_internal_maths::vectorOrientationSimilarity ( CSymList->at(grIt)[1], CSymList->at(grIt)[2], CSymList->at(grIt)[3], X, Y, Z, tolerance ) )
2672  {
2673  ret = false;
2674  break;
2675  }
2676  }
2677  }
2678  std::cout << std::endl;
2679 
2680  //================================================ Done
2681  return ( ret );
2682 
2683 }

◆ multiplyGroupElementMatrices()

std::vector< proshade_double > ProSHADE_internal_maths::multiplyGroupElementMatrices ( std::vector< proshade_double > *  el1,
std::vector< proshade_double > *  el2 
)

This function computes matrix multiplication using the ProSHADE group element matrix format as input and output.

Parameters
[in]el1Group element as rotation matrix in the group element matrix format.
[in]el2Group element as rotation matrix in the group element matrix format.
[out]retMatrix in the group element format resulting from the input matrices multiplication.

Definition at line 2141 of file ProSHADE_maths.cpp.

2142 {
2143  //================================================ Initialise variables
2144  std::vector< proshade_double > ret;
2145 
2146  //================================================ Compute
2147  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(0) * el2->at(0) ) +
2148  ( el1->at(1) * el2->at(3) ) +
2149  ( el1->at(2) * el2->at(6) ) );
2150  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(0) * el2->at(1) ) +
2151  ( el1->at(1) * el2->at(4) ) +
2152  ( el1->at(2) * el2->at(7) ) );
2153  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(0) * el2->at(2) ) +
2154  ( el1->at(1) * el2->at(5) ) +
2155  ( el1->at(2) * el2->at(8) ) );
2156 
2157  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(3) * el2->at(0) ) +
2158  ( el1->at(4) * el2->at(3) ) +
2159  ( el1->at(5) * el2->at(6) ) );
2160  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(3) * el2->at(1) ) +
2161  ( el1->at(4) * el2->at(4) ) +
2162  ( el1->at(5) * el2->at(7) ) );
2163  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(3) * el2->at(2) ) +
2164  ( el1->at(4) * el2->at(5) ) +
2165  ( el1->at(5) * el2->at(8) ) );
2166 
2167  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(6) * el2->at(0) ) +
2168  ( el1->at(7) * el2->at(3) ) +
2169  ( el1->at(8) * el2->at(6) ) );
2170  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(6) * el2->at(1) ) +
2171  ( el1->at(7) * el2->at(4) ) +
2172  ( el1->at(8) * el2->at(7) ) );
2173  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(6) * el2->at(2) ) +
2174  ( el1->at(7) * el2->at(5) ) +
2175  ( el1->at(8) * el2->at(8) ) );
2176 
2177 
2178  //================================================ Done
2179  return ( ret );
2180 
2181 }

◆ multiplyTwoSquareMatrices()

void ProSHADE_internal_maths::multiplyTwoSquareMatrices ( proshade_double *  A,
proshade_double *  B,
proshade_double *  res,
proshade_unsign  dim = 3 
)

Function to compute matrix multiplication.

Parameters
[in]AThe left matrix of the matrix multiplication to be solved.
[in]BThe right matrix of the matrix multiplication to be solved. (Assuming it already has been transposed).
[in]resMatrix containing the results.
[in]dimThe dimension of all the matrices (i.e. assuming square dim*dim matrices).
Warning
This function assumes the second matrix has been transposed already!

Definition at line 1613 of file ProSHADE_maths.cpp.

1614 {
1615  //================================================ Set res to 0.0s
1616  for ( proshade_unsign iter = 0; iter < 9; iter++ ) { res[iter] = 0.0; }
1617 
1618  //================================================ Compute the matrix multiplication
1619  for ( proshade_unsign row = 0; row < dim; row++ )
1620  {
1621  for ( proshade_unsign col = 0; col < dim; col++ )
1622  {
1623  for ( proshade_unsign inner = 0; inner < dim; inner++ )
1624  {
1625  res[(row*dim)+col] += A[(inner*dim)+row] * B[(col*dim)+inner];
1626  }
1627  }
1628  }
1629 
1630  //================================================ Done
1631  return ;
1632 
1633 }

◆ normalDistributionValue()

proshade_double ProSHADE_internal_maths::normalDistributionValue ( proshade_double  mean,
proshade_double  standardDev,
proshade_double  value 
)

Function to the heiht of normal distribution given by mean and standard deviation for a given value.

Parameters
[in]meanThe mean of the normal distribution.
[in]standardDevThe standard deviation of the normal distribution.
[in]valueThe value on the axis for which the height of the normal distribution is to be obtained.
[out]XThe height of the normal distribution at point given by the value.

Definition at line 1688 of file ProSHADE_maths.cpp.

1689 {
1690  //================================================ Compute and return
1691  return ( ( 1.0 / sqrt ( 2.0 * M_PI * pow(standardDev,2.0) ) ) * std::exp ( - pow( value - mean, 2.0 ) / 2.0 * pow(standardDev,2.0) ) );
1692 
1693 }

◆ optimiseAxisBiCubicInterpolation()

void ProSHADE_internal_maths::optimiseAxisBiCubicInterpolation ( proshade_double *  bestLattitude,
proshade_double *  bestLongitude,
proshade_double *  bestSum,
std::vector< proshade_unsign > *  sphereList,
std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *  sphereMappedRotFun,
proshade_double  step = 0.05 
)

This function provides axis optimisation given starting lattitude and longitude indices.

This function takes the initial lattitude and longitude indices as well as the current best sum over all appropriate spheres and the list of the spheres and proceeds to use bi-cubic interpolation and a sort of gradient ascend algorithm to search the space around the given indices for interpolated values, which would have higher sum of the rotation function values than the initial position. If any improvement is found, it will over-write the input variables.

Parameters
[in]bestLattitudeProshade double pointer to variable containing the best lattitude index value and to which the optimised result will be saved into.
[in]bestLongitudeProshade double pointer to variable containing the best longitude index value and to which the optimised result will be saved into.
[in]bestSumProshade double pointer to variable containing the best position rotation function values sum and to which the optimised result will be saved into.
[in]sphereListA vector containing the list of spheres which form the set for this symmetry.
[in]stepThe size of the step.

Definition at line 2296 of file ProSHADE_maths.cpp.

2297 {
2298  //================================================ Initialise variables
2299  proshade_double lonM, lonP, latM, latP, movSum;
2300  std::vector<proshade_double> latVals ( 3 );
2301  std::vector<proshade_double> lonVals ( 3 );
2302  proshade_double learningRate = 0.1;
2303  proshade_double prevVal = *bestSum;
2304  proshade_double valChange = 999.9;
2305  proshade_double origBestLat = std::round ( *bestLattitude );
2306  proshade_double origBestLon = std::round ( *bestLongitude );
2307  proshade_double tmpVal;
2308 
2309  //================================================ Initialise interpolators in all directions around the point of interest
2310  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsMinusMinus;
2311  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsMinusPlus;
2312  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsPlusMinus;
2313  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsPlusPlus;
2314  prepareBiCubicInterpolatorsMinusMinus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsMinusMinus, sphereMappedRotFun );
2315  prepareBiCubicInterpolatorsMinusPlus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsMinusPlus, sphereMappedRotFun );
2316  prepareBiCubicInterpolatorsPlusMinus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsPlusMinus, sphereMappedRotFun );
2317  prepareBiCubicInterpolatorsPlusPlus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsPlusPlus, sphereMappedRotFun );
2318 
2319  //================================================ Start the pseudo gradient ascent (while there is some change)
2320  while ( valChange > 0.0001 )
2321  {
2322  //============================================ Find the surrounding points to the currently best position
2323  lonM = *bestLongitude - step;
2324  lonP = *bestLongitude + step;
2325  latM = *bestLattitude - step;
2326  latP = *bestLattitude + step;
2327 
2328  //============================================ Deal with optimising outside of prepared range - recursion
2329  if ( latM < ( origBestLat - 1.0 ) ) { tmpVal = *bestLattitude; *bestLattitude = origBestLat - 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( *bestLattitude == origBestLat - 1.0 ) { *bestLattitude = tmpVal; } break; }
2330  if ( latP > ( origBestLat + 1.0 ) ) { tmpVal = *bestLattitude; *bestLattitude = origBestLat + 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( *bestLattitude == origBestLat + 1.0 ) { *bestLattitude = tmpVal; } break; }
2331  if ( lonM < ( origBestLon - 1.0 ) ) { tmpVal = *bestLongitude; *bestLongitude = origBestLon - 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( *bestLongitude == origBestLon - 1.0 ) { *bestLongitude = tmpVal; } break; }
2332  if ( lonP > ( origBestLon + 1.0 ) ) { tmpVal = *bestLongitude; *bestLongitude = origBestLon + 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( *bestLongitude == origBestLon + 1.0 ) { *bestLongitude = tmpVal; } break; }
2333 
2334  //============================================ Prepare vectors of tested positions
2335  latVals.at(0) = latM; latVals.at(1) = *bestLattitude; latVals.at(2) = latP;
2336  lonVals.at(0) = lonM; lonVals.at(1) = *bestLongitude; lonVals.at(2) = lonP;
2337 
2338  //============================================ Find the best change
2339  for ( proshade_unsign laIt = 0; laIt < static_cast<proshade_unsign> ( latVals.size() ); laIt++ )
2340  {
2341  for ( proshade_unsign loIt = 0; loIt < static_cast<proshade_unsign> ( lonVals.size() ); loIt++ )
2342  {
2343  //==================================== For this combination of lat and lon, find sum over spheres
2344  movSum = 1.0;
2345  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sphereList->size() ); iter++ )
2346  {
2347  //================================ Interpolate using correct interpolators
2348  if ( ( latVals.at(laIt) <= origBestLat ) && ( lonVals.at(loIt) <= origBestLon ) ) { movSum += interpolsMinusMinus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2349  if ( ( latVals.at(laIt) <= origBestLat ) && ( lonVals.at(loIt) > origBestLon ) ) { movSum += interpolsMinusPlus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2350  if ( ( latVals.at(laIt) > origBestLat ) && ( lonVals.at(loIt) <= origBestLon ) ) { movSum += interpolsPlusMinus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2351  if ( ( latVals.at(laIt) > origBestLat ) && ( lonVals.at(loIt) > origBestLon ) ) { movSum += interpolsPlusPlus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2352  }
2353 
2354  //==================================== If position has improved, save it
2355  if ( *bestSum < movSum )
2356  {
2357  *bestSum = movSum;
2358  *bestLongitude = lonVals.at(loIt);
2359  *bestLattitude = latVals.at(laIt);
2360  }
2361  }
2362  }
2363 
2364  //============================================ Prepare for next iteration
2365  valChange = std::floor ( 100000.0 * ( *bestSum - prevVal ) ) / 100000.0;
2366  prevVal = std::floor ( 100000.0 * ( *bestSum ) ) / 100000.0;
2367  step = std::max ( ( valChange / step ) * learningRate, 0.01 );
2368  if ( learningRate >= 0.02 ) { learningRate -= 0.01; }
2369 
2370  }
2371 
2372  //================================================ Release interpolators memory
2373  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsMinusMinus.size() ); intIt++ ) { delete interpolsMinusMinus.at(intIt); }
2374  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsMinusPlus.size() ); intIt++ ) { delete interpolsMinusPlus.at(intIt); }
2375  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsPlusMinus.size() ); intIt++ ) { delete interpolsPlusMinus.at(intIt); }
2376  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsPlusPlus.size() ); intIt++ ) { delete interpolsPlusPlus.at(intIt); }
2377 
2378  //================================================ Done
2379  return ;
2380 }

◆ pearsonCorrCoeff()

proshade_double ProSHADE_internal_maths::pearsonCorrCoeff ( proshade_double *  valSet1,
proshade_double *  valSet2,
proshade_unsign  length 
)

Function for computing the Pearson's correlation coefficient.

This function takes two numerical arrays of same length and proceeds to compute the Pearson's correlation coefficient, which it then returns.

Parameters
[in]valSet1This is the set of x-values.
[in]valSet2This is the set of y-values.
[in]lengthThe length of both arrays (both arrays have to have the same length).
[out]XThe Pearson's correlation coefficient value.

Definition at line 244 of file ProSHADE_maths.cpp.

245 {
246  //================================================ Find vector means
247  proshade_double xMean = 0.0;
248  proshade_double yMean = 0.0;
249  proshade_double zeroCount = 0.0;
250  for ( proshade_unsign iter = 0; iter < length; iter++ )
251  {
252  xMean += valSet1[iter];
253  yMean += valSet2[iter];
254  }
255  xMean /= static_cast<proshade_double> ( length - zeroCount );
256  yMean /= static_cast<proshade_double> ( length - zeroCount );
257 
258  //================================================ Get Pearson's correlation coefficient
259  proshade_double xmmymm = 0.0;
260  proshade_double xmmsq = 0.0;
261  proshade_double ymmsq = 0.0;
262  for ( proshade_unsign iter = 0; iter < length; iter++ )
263  {
264  xmmymm += ( valSet1[iter] - xMean ) * ( valSet2[iter] - yMean );
265  xmmsq += pow( valSet1[iter] - xMean, 2.0 );
266  ymmsq += pow( valSet2[iter] - yMean, 2.0 );
267  }
268 
269  proshade_double ret = xmmymm / ( sqrt(xmmsq) * sqrt(ymmsq) );
270 
271  //================================================ Done
272  if ( std::isnan ( ret ) ) { return ( 0.0 ); }
273  return ( ret );
274 
275 }

◆ prepareBiCubicInterpolatorsMinusMinus()

void ProSHADE_internal_maths::prepareBiCubicInterpolatorsMinusMinus ( proshade_double  bestLattitude,
proshade_double  bestLongitude,
std::vector< proshade_unsign > *  sphereList,
std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *  interpols,
std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *  sphereMappedRotFun 
)

This function prepares the interpolation objects for the bi-cubic interpolation.

This function takes the position around which the interpolation is to be done and proceeds to create the interpolator objects for bi-cubic interpolation in the – direction (i.e. when both interpolated values will be lower than the best lattitude and longitude) using the correct spheres in the correct ranges.

Parameters
[in]bestLattitudeThe lattitude index value around which interpolation is to be prepared.
[in]bestLongitudeThe longitude index value around which interpolation is to be prepared.
[in]sphereListA vector containing the list of spheres which form the set for this symmetry.
[in]interpolsA pointer to a vector of ProSHADE interpolator objects to which the interpolators will be saved into.

Definition at line 2393 of file ProSHADE_maths.cpp.

2394 {
2395  //================================================ Initialise local variables
2396  proshade_signed latHlp, lonHlp;
2397  proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2398 
2399  //================================================ Prepare the interpolator objects for interpolation around the position
2400  for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2401  {
2402  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along first dimension)
2403  proshade_double** interpGrid = new proshade_double*[4];
2404  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid, __FILE__, __LINE__, __func__ );
2405 
2406  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along second dimension)
2407  for ( proshade_unsign iter = 0; iter < 4; iter++ )
2408  {
2409  interpGrid[iter] = new proshade_double[4];
2410  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid[iter], __FILE__, __LINE__, __func__ );
2411  }
2412 
2413  //============================================ Fill in the value grid on which the interpolation is to be done
2414  for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2415  {
2416  for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2417  {
2418  latHlp = bestLattitude - 2 + latIt; if ( latHlp < 0.0 ) { latHlp += angDim; } if ( latHlp >= angDim ) { latHlp -= angDim; }
2419  lonHlp = bestLongitude - 2 + lonIt; if ( lonHlp < 0.0 ) { lonHlp += angDim; } if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2420  interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2421  }
2422  }
2423 
2424  //============================================ Create the interpolators
2425  ProSHADE_internal_maths::BicubicInterpolator* biCubInterp = new ProSHADE_internal_maths::BicubicInterpolator ( interpGrid, bestLattitude - 1.0, bestLongitude - 1.0 );
2426  interpols->emplace_back ( biCubInterp );
2427 
2428  //============================================ Release memory
2429  for ( proshade_unsign iter = 0; iter < 4; iter++ ) { delete[] interpGrid[iter]; }
2430  delete[] interpGrid;
2431  }
2432 
2433  //================================================ Done
2434  return ;
2435 }

◆ prepareBiCubicInterpolatorsMinusPlus()

void ProSHADE_internal_maths::prepareBiCubicInterpolatorsMinusPlus ( proshade_double  bestLattitude,
proshade_double  bestLongitude,
std::vector< proshade_unsign > *  sphereList,
std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *  interpols,
std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *  sphereMappedRotFun 
)

This function prepares the interpolation objects for the bi-cubic interpolation.

This function takes the position around which the interpolation is to be done and proceeds to create the interpolator objects for bi-cubic interpolation in the -+ direction (i.e. when interpolated lattitude is lower than best lattitude, but interpolated longitude is higher than best longitude) using the correct spheres in the correct ranges.

Parameters
[in]bestLattitudeThe lattitude index value around which interpolation is to be prepared.
[in]bestLongitudeThe longitude index value around which interpolation is to be prepared.
[in]sphereListA vector containing the list of spheres which form the set for this symmetry.
[in]interpolsA pointer to a vector of ProSHADE interpolator objects to which the interpolators will be saved into.

Definition at line 2448 of file ProSHADE_maths.cpp.

2449 {
2450  //================================================ Initialise local variables
2451  proshade_signed latHlp, lonHlp;
2452  proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2453 
2454  //================================================ Prepare the interpolator objects for interpolation around the position
2455  for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2456  {
2457  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along first dimension)
2458  proshade_double** interpGrid = new proshade_double*[4];
2459  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid, __FILE__, __LINE__, __func__ );
2460 
2461  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along second dimension)
2462  for ( proshade_unsign iter = 0; iter < 4; iter++ )
2463  {
2464  interpGrid[iter] = new proshade_double[4];
2465  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid[iter], __FILE__, __LINE__, __func__ );
2466  }
2467 
2468  //============================================ Fill in the value grid on which the interpolation is to be done
2469  for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2470  {
2471  for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2472  {
2473  latHlp = bestLattitude - 2 + latIt; if ( latHlp < 0.0 ) { latHlp += angDim; } if ( latHlp >= angDim ) { latHlp -= angDim; }
2474  lonHlp = bestLongitude - 1 + lonIt; if ( lonHlp < 0.0 ) { lonHlp += angDim; } if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2475  interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2476  }
2477  }
2478 
2479  //============================================ Create the interpolators
2480  ProSHADE_internal_maths::BicubicInterpolator* biCubInterp = new ProSHADE_internal_maths::BicubicInterpolator ( interpGrid, bestLattitude - 1.0, bestLongitude );
2481  interpols->emplace_back ( biCubInterp );
2482 
2483  //============================================ Release memory
2484  for ( proshade_unsign iter = 0; iter < 4; iter++ ) { delete[] interpGrid[iter]; }
2485  delete[] interpGrid;
2486  }
2487 
2488  //================================================ Done
2489  return ;
2490 }

◆ prepareBiCubicInterpolatorsPlusMinus()

void ProSHADE_internal_maths::prepareBiCubicInterpolatorsPlusMinus ( proshade_double  bestLattitude,
proshade_double  bestLongitude,
std::vector< proshade_unsign > *  sphereList,
std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *  interpols,
std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *  sphereMappedRotFun 
)

This function prepares the interpolation objects for the bi-cubic interpolation.

This function takes the position around which the interpolation is to be done and proceeds to create the interpolator objects for bi-cubic interpolation in the +- direction (i.e. when interpolated lattitude is higher than best lattitude, but interpolated longitude is lower than best longitude) using the correct spheres in the correct ranges.

Parameters
[in]bestLattitudeThe lattitude index value around which interpolation is to be prepared.
[in]bestLongitudeThe longitude index value around which interpolation is to be prepared.
[in]sphereListA vector containing the list of spheres which form the set for this symmetry.
[in]interpolsA pointer to a vector of ProSHADE interpolator objects to which the interpolators will be saved into.

Definition at line 2503 of file ProSHADE_maths.cpp.

2504 {
2505  //================================================ Initialise local variables
2506  proshade_signed latHlp, lonHlp;
2507  proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2508 
2509  //================================================ Prepare the interpolator objects for interpolation around the position
2510  for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2511  {
2512  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along first dimension)
2513  proshade_double** interpGrid = new proshade_double*[4];
2514  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid, __FILE__, __LINE__, __func__ );
2515 
2516  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along second dimension)
2517  for ( proshade_unsign iter = 0; iter < 4; iter++ )
2518  {
2519  interpGrid[iter] = new proshade_double[4];
2520  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid[iter], __FILE__, __LINE__, __func__ );
2521  }
2522 
2523  //============================================ Fill in the value grid on which the interpolation is to be done
2524  for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2525  {
2526  for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2527  {
2528  latHlp = bestLattitude - 1 + latIt; if ( latHlp < 0.0 ) { latHlp += angDim; } if ( latHlp >= angDim ) { latHlp -= angDim; }
2529  lonHlp = bestLongitude - 2 + lonIt; if ( lonHlp < 0.0 ) { lonHlp += angDim; } if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2530  interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2531  }
2532  }
2533 
2534  //============================================ Create the interpolators
2535  ProSHADE_internal_maths::BicubicInterpolator* biCubInterp = new ProSHADE_internal_maths::BicubicInterpolator ( interpGrid, bestLattitude, bestLongitude - 1.0 );
2536  interpols->emplace_back ( biCubInterp );
2537 
2538  //============================================ Release memory
2539  for ( proshade_unsign iter = 0; iter < 4; iter++ ) { delete[] interpGrid[iter]; }
2540  delete[] interpGrid;
2541  }
2542 
2543  //================================================ Done
2544  return ;
2545 }

◆ prepareBiCubicInterpolatorsPlusPlus()

void ProSHADE_internal_maths::prepareBiCubicInterpolatorsPlusPlus ( proshade_double  bestLattitude,
proshade_double  bestLongitude,
std::vector< proshade_unsign > *  sphereList,
std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *  interpols,
std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *  sphereMappedRotFun 
)

This function prepares the interpolation objects for the bi-cubic interpolation.

This function takes the position around which the interpolation is to be done and proceeds to create the interpolator objects for bi-cubic interpolation in the ++ direction (i.e. when both interpolated values will be larger than the best lattitude and longitude) using the correct spheres in the correct ranges.

Parameters
[in]bestLattitudeThe lattitude index value around which interpolation is to be prepared.
[in]bestLongitudeThe longitude index value around which interpolation is to be prepared.
[in]sphereListA vector containing the list of spheres which form the set for this symmetry.
[in]interpolsA pointer to a vector of ProSHADE interpolator objects to which the interpolators will be saved into.

Definition at line 2558 of file ProSHADE_maths.cpp.

2559 {
2560  //================================================ Initialise local variables
2561  proshade_signed latHlp, lonHlp;
2562  proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2563 
2564  //================================================ Prepare the interpolator objects for interpolation around the position
2565  for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2566  {
2567  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along first dimension)
2568  proshade_double** interpGrid = new proshade_double*[4];
2569  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid, __FILE__, __LINE__, __func__ );
2570 
2571  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along second dimension)
2572  for ( proshade_unsign iter = 0; iter < 4; iter++ )
2573  {
2574  interpGrid[iter] = new proshade_double[4];
2575  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid[iter], __FILE__, __LINE__, __func__ );
2576  }
2577 
2578  //============================================ Fill in the value grid on which the interpolation is to be done
2579  for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2580  {
2581  for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2582  {
2583  latHlp = bestLattitude - 1 + latIt; if ( latHlp < 0.0 ) { latHlp += angDim; } if ( latHlp >= angDim ) { latHlp -= angDim; }
2584  lonHlp = bestLongitude - 1 + lonIt; if ( lonHlp < 0.0 ) { lonHlp += angDim; } if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2585  interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2586  }
2587  }
2588 
2589  //============================================ Create the interpolators
2590  ProSHADE_internal_maths::BicubicInterpolator* biCubInterp = new ProSHADE_internal_maths::BicubicInterpolator ( interpGrid, bestLattitude, bestLongitude );
2591  interpols->emplace_back ( biCubInterp );
2592 
2593  //============================================ Release memory
2594  for ( proshade_unsign iter = 0; iter < 4; iter++ ) { delete[] interpGrid[iter]; }
2595  delete[] interpGrid;
2596  }
2597 
2598  //================================================ Done
2599  return ;
2600 }

◆ primeFactorsDecomp()

std::vector< proshade_signed > ProSHADE_internal_maths::primeFactorsDecomp ( proshade_signed  number)

Function to find prime factors of an integer.

Parameters
[in]numberA single integer number to be decomposed into its prime factors.

Definition at line 1639 of file ProSHADE_maths.cpp.

1640 {
1641  //================================================ Initialise variables
1642  std::vector < proshade_signed > ret;
1643 
1644  //================================================ Deal with negative numbers
1645  bool changeSign = false;
1646  if ( number < 0 ) { changeSign = true; number = -number; }
1647 
1648  //================================================ Deal with zero and one
1649  if ( number == 0 ) { ProSHADE_internal_misc::addToSignedVector ( &ret, 0 ); return ( ret ); }
1650  if ( number == 1 ) { ProSHADE_internal_misc::addToSignedVector ( &ret, 1 ); return ( ret ); }
1651 
1652  //================================================ Divide by 2 as long as you can
1653  while ( number % 2 == 0 )
1654  {
1656  number = number / 2;
1657  }
1658 
1659  //================================================ Check all odd numbers up to the square root
1660  for ( proshade_unsign posDiv = 3; posDiv <= sqrt ( number ); posDiv += 2)
1661  {
1662  // If posDiv is a divisor of the number, save the result
1663  while ( number % posDiv == 0 )
1664  {
1666  number = number / posDiv;
1667  }
1668  }
1669 
1670  //================================================ If the number was a large prime number, save it as it is
1671  if ( number > 2 ) { ProSHADE_internal_misc::addToSignedVector ( &ret, number ); }
1672 
1673  //================================================ Finish dealing with negative numbers
1674  if ( changeSign ) { ret.at(0) = -ret.at(0); }
1675 
1676  //================================================ Done
1677  return ( ret );
1678 
1679 }

◆ rotationMatrixSimilarity()

bool ProSHADE_internal_maths::rotationMatrixSimilarity ( std::vector< proshade_double > *  mat1,
std::vector< proshade_double > *  mat2,
proshade_double  tolerance = 0.1 
)

This function compares the distance between two rotation matrices and decides if they are similar using tolerance.

This function computes the distance between two rotation matrices, specifically by computing the trace of (R1 * R2^T). This measure will be 3.0 if the two matrices are identical and will decrease the more the rotation matrices difference diverges from identity. Therefore, from this trace 3.0 is subtracted and the absolute value of the result is compared to the tolerance. If the difference is less than the tolerance, true is returned, while false is returned otherwise.

Parameters
[in]mat1Vector of 9 numbers representing first rotation matrix.
[in]mat1Vector of 9 numbers representing second rotation matrix.
[in]toleranceDouble number representing the maximum allowed error on the distance.
[out]resBoolean decision if the two matrices are similar or not.

Definition at line 2195 of file ProSHADE_maths.cpp.

2196 {
2197  //================================================ Initialise variables
2198  bool ret = false;
2199 
2200  //================================================ Compute trace of mat1 * mat2^T
2201  proshade_double trace = ( mat1->at(0) * mat2->at(0) ) + ( mat1->at(1) * mat2->at(1) ) + ( mat1->at(2) * mat2->at(2) );
2202  trace += ( mat1->at(3) * mat2->at(3) ) + ( mat1->at(4) * mat2->at(4) ) + ( mat1->at(5) * mat2->at(5) );
2203  trace += ( mat1->at(6) * mat2->at(6) ) + ( mat1->at(7) * mat2->at(7) ) + ( mat1->at(8) * mat2->at(8) );
2204 
2205  //================================================ Subtract 3 (so that we would have 0 in case of idenity matrix)
2206  trace -= 3.0;
2207 
2208  //================================================ Compare to tolerance
2209  if ( tolerance > std::abs ( trace ) ) { ret = true; }
2210 
2211  //================================================ Done
2212  return ( ret );
2213 
2214 }

◆ transpose3x3MatrixInPlace()

void ProSHADE_internal_maths::transpose3x3MatrixInPlace ( proshade_double *  mat)

Transposes 3x3 matrix in place.

Parameters
[in]matThe matrix to be transposed.

Definition at line 1843 of file ProSHADE_maths.cpp.

1844 {
1845  //================================================ Initialise variables
1846  proshade_double tmp;
1847 
1848  //================================================ Transpose the non-diagonal values
1849  tmp = mat[1];
1850  mat[1] = mat[3];
1851  mat[3] = tmp;
1852 
1853  tmp = mat[2];
1854  mat[2] = mat[6];
1855  mat[6] = tmp;
1856 
1857  tmp = mat[5];
1858  mat[5] = mat[7];
1859  mat[7] = tmp;
1860 
1861  //================================================ Done
1862  return ;
1863 
1864 }

◆ vectorMeanAndSD()

void ProSHADE_internal_maths::vectorMeanAndSD ( std::vector< proshade_double > *  vec,
proshade_double *&  ret 
)

Function to get vector mean and standard deviation.

This function takes a pointer to a vector of proshade_double's and returns the mean and standard deviation of such vector.

Parameters
[in]vecPointer to a vector of proshade_double's for which mean and sd should be obtained.
[in]retPointer to array of 2 proshade_double's, which will be the return values - first mean and second sd.

Definition at line 121 of file ProSHADE_maths.cpp.

122 {
123  //================================================ Get mean
124  ret[0] = std::accumulate ( vec->begin(), vec->end(), 0.0 ) / static_cast<proshade_double> ( vec->size() );
125 
126  //================================================ Get standard deviation
127  proshade_double squaredSum = std::inner_product ( vec->begin(), vec->end(), vec->begin(), 0.0 );
128  ret[1] = std::sqrt ( ( squaredSum / static_cast<proshade_double> ( vec->size() ) ) - std::pow ( ret[0], 2.0 ) );
129 
130  //================================================ Check for NaN's
131  if ( ret[0] != ret[0] ) { ret[0] = 0.0; }
132  if ( ret[1] != ret[1] ) { ret[1] = 0.0; }
133 
134  //================================================ Return
135  return ;
136 
137 }

◆ vectorMedianAndIQR()

void ProSHADE_internal_maths::vectorMedianAndIQR ( std::vector< proshade_double > *  vec,
proshade_double *&  ret 
)

Function to get vector median and inter-quartile range.

This function takes a pointer to a vector of proshade_double's and returns the median and the inter-quartile range of such vector.

Parameters
[in]vecPointer to a vector of proshade_double's for which median and IQR should be obtained.
[in]retPointer to array of 2 proshade_double's, which will be the return values - first median and second IQR.

Definition at line 147 of file ProSHADE_maths.cpp.

148 {
149  //================================================ Sanity check
150  if ( vec->size() < 3 ) { ret[0] = 0.0; ret[1] = 0.0; return; }
151 
152  //================================================ Sort the vector
153  std::sort ( vec->begin(), vec->end() );
154 
155  //================================================ Get median
156  if ( static_cast<proshade_unsign> ( vec->size() ) % 2 == 0)
157  {
158  ret[0] = ( vec->at( ( static_cast<proshade_unsign> ( vec->size() ) / 2 ) - 1 ) +
159  vec->at( static_cast<proshade_unsign> ( vec->size() ) / 2 ) ) / 2.0;
160  }
161  else
162  {
163  ret[0] = vec->at( static_cast<proshade_unsign> ( vec->size() ) / 2 );
164  }
165 
166  //================================================ Get first and third quartile
167  proshade_double Q1, Q3;
168  if ( static_cast<proshade_unsign> ( vec->size() ) % 2 == 0)
169  {
170  Q1 = ( vec->at( ( static_cast<proshade_unsign> ( vec->size() ) / 4 ) - 1 ) +
171  vec->at( static_cast<proshade_unsign> ( vec->size() ) / 4 ) ) / 2.0;
172  Q3 = ( vec->at( ( ( static_cast<proshade_unsign> ( vec->size() ) / 4 ) * 3 ) - 1 ) +
173  vec->at( ( static_cast<proshade_unsign> ( vec->size() ) / 4 ) * 3 ) ) / 2.0;
174  }
175  else
176  {
177  Q1 = vec->at( static_cast<proshade_unsign> ( vec->size() ) / 4 );
178  Q3 = vec->at( ( static_cast<proshade_unsign> ( vec->size() ) / 4 ) * 3 );
179  }
180 
181  //================================================ And now save the IQR
182  ret[1] = Q3 - Q1;
183 
184  //================================================ Return
185  return ;
186 
187 }

◆ vectorOrientationSimilarity()

bool ProSHADE_internal_maths::vectorOrientationSimilarity ( proshade_double  a1,
proshade_double  a2,
proshade_double  a3,
proshade_double  b1,
proshade_double  b2,
proshade_double  b3,
proshade_double  tolerance = 0.1 
)

This function compares two vectors using cosine distance and decides if they are similar using tolerance.

This function computes the distance between two vectors, specifically by computing the cosine distance ( ( dot( A, B ) ) / ( mag(A) x mag(B) ) ). This measure will be 1.0 if the two vectors are identically oriented, 0.0 if they are perpendicular and -1.0 if they have opposite direction. Given that opposite direction must be regarded as same direction with opposite angles for symmetry axes detection purposes, this function uses the absolute value of this measure and checks if the two supplied vectors (supplied element by element) have cosine distance within 1.0 - tolerance, returning true if they do and false otherwise.

Parameters
[in]a1The first element of the first vector.
[in]a2The second element of the first vector.
[in]a3The third element of the first vector.
[in]b1The first element of the second vector.
[in]b2The second element of the second vector.
[in]b3The third element of the second vector.
[in]toleranceThe allowed difference of the distance measure from the 1.0 for the vectors to still be considered similar.
[out]resBoolean decision if the two vectors are similar or not.

Definition at line 2232 of file ProSHADE_maths.cpp.

2233 {
2234  //================================================ Initialise variables
2235  bool ret = false;
2236 
2237  //================================================ Cosine distance
2238  proshade_double cosDist = ( ( a1 * b1 ) + ( a2 * b2 ) + ( a3 * b3 ) ) /
2239  ( sqrt( pow( a1, 2.0 ) + pow( a2, 2.0 ) + pow( a3, 2.0 ) ) *
2240  sqrt( pow( b1, 2.0 ) + pow( b2, 2.0 ) + pow( b3, 2.0 ) ) );
2241 
2242  //================================================ Compare the absolute value of distance to 1.0 - tolerance
2243  if ( std::abs( cosDist ) > ( 1.0 - tolerance ) ) { ret = true; }
2244 
2245  //================================================ Done
2246  return ( ret );
2247 
2248 }

◆ vectorOrientationSimilaritySameDirection()

bool ProSHADE_internal_maths::vectorOrientationSimilaritySameDirection ( proshade_double  a1,
proshade_double  a2,
proshade_double  a3,
proshade_double  b1,
proshade_double  b2,
proshade_double  b3,
proshade_double  tolerance = 0.1 
)

This function compares two vectors using cosine distance and decides if they are similar using tolerance.

This function computes the distance between two vectors, specifically by computing the cosine distance ( ( dot( A, B ) ) / ( mag(A) x mag(B) ) ). This measure will be 1.0 if the two vectors are identically oriented, 0.0 if they are perpendicular and -1.0 if they have opposite direction. Given that opposite direction must be regarded as different vector for peak detection purposes (spheres with different angles are covered separately), this function does not use the absolute value of this measure and checks if the two supplied vectors (supplied element by element) have cosine distance within 1.0 - tolerance, returning true if they do and false otherwise.

Parameters
[in]a1The first element of the first vector.
[in]a2The second element of the first vector.
[in]a3The third element of the first vector.
[in]b1The first element of the second vector.
[in]b2The second element of the second vector.
[in]b3The third element of the second vector.
[in]toleranceThe allowed difference of the distance measure from the 1.0 for the vectors to still be considered similar.
[out]resBoolean decision if the two vectors are similar or not.

Definition at line 2266 of file ProSHADE_maths.cpp.

2267 {
2268  //================================================ Initialise variables
2269  bool ret = false;
2270 
2271  //================================================ Cosine distance
2272  proshade_double cosDist = ( ( a1 * b1 ) + ( a2 * b2 ) + ( a3 * b3 ) ) /
2273  ( sqrt( pow( a1, 2.0 ) + pow( a2, 2.0 ) + pow( a3, 2.0 ) ) *
2274  sqrt( pow( b1, 2.0 ) + pow( b2, 2.0 ) + pow( b3, 2.0 ) ) );
2275 
2276  //================================================ Compare the absolute value of distance to 1.0 - tolerance
2277  if ( cosDist > ( 1.0 - tolerance ) ) { ret = true; }
2278 
2279  //================================================ Done
2280  return ( ret );
2281 
2282 }
ProSHADE_internal_maths::prepareBiCubicInterpolatorsMinusPlus
void prepareBiCubicInterpolatorsMinusPlus(proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
This function prepares the interpolation objects for the bi-cubic interpolation.
Definition: ProSHADE_maths.cpp:2448
ProSHADE_internal_maths::computeCrossProduct
proshade_double * computeCrossProduct(proshade_double *x1, proshade_double *y1, proshade_double *z1, proshade_double *x2, proshade_double *y2, proshade_double *z2)
Simple 3D vector cross product computation.
Definition: ProSHADE_maths.cpp:1737
ProSHADE_exception
This class is the representation of ProSHADE exception.
Definition: ProSHADE_exceptions.hpp:37
ProSHADE_internal_maths::prepareBiCubicInterpolatorsPlusMinus
void prepareBiCubicInterpolatorsPlusMinus(proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
This function prepares the interpolation objects for the bi-cubic interpolation.
Definition: ProSHADE_maths.cpp:2503
ProSHADE_internal_maths::evaluateGLSeries
proshade_double evaluateGLSeries(proshade_double *series, proshade_double target, proshade_unsign terms)
This function evaluates the Taylor expansion.
Definition: ProSHADE_maths.cpp:447
ProSHADE_internal_maths::getRotationMatrixFromEulerZXZAngles
void getRotationMatrixFromEulerZXZAngles(proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma, proshade_double *matrix)
Function to find the rotation matrix from Euler angles (ZXZ convention).
Definition: ProSHADE_maths.cpp:1005
ProSHADE_internal_maths::getGLFirstEvenRoot
void getGLFirstEvenRoot(proshade_double polyAtZero, proshade_unsign order, proshade_double *abscAtZero, proshade_double *weighAtZero, proshade_unsign taylorSeriesCap)
This function finds the first root for Legendre polynomials of odd order.
Definition: ProSHADE_maths.cpp:384
ProSHADE_internal_maths::optimiseAxisBiCubicInterpolation
void optimiseAxisBiCubicInterpolation(proshade_double *bestLattitude, proshade_double *bestLongitude, proshade_double *bestSum, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun, proshade_double step=0.05)
This function provides axis optimisation given starting lattitude and longitude indices.
Definition: ProSHADE_maths.cpp:2296
ProSHADE_internal_maths::vectorOrientationSimilarity
bool vectorOrientationSimilarity(proshade_double a1, proshade_double a2, proshade_double a3, proshade_double b1, proshade_double b2, proshade_double b3, proshade_double tolerance=0.1)
This function compares two vectors using cosine distance and decides if they are similar using tolera...
Definition: ProSHADE_maths.cpp:2232
ProSHADE_internal_maths::prepareBiCubicInterpolatorsMinusMinus
void prepareBiCubicInterpolatorsMinusMinus(proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
This function prepares the interpolation objects for the bi-cubic interpolation.
Definition: ProSHADE_maths.cpp:2393
ProSHADE_internal_misc::addToDoubleVector
void addToDoubleVector(std::vector< proshade_double > *vecToAddTo, proshade_double elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:77
ProSHADE_internal_maths::computeDotProduct
proshade_double computeDotProduct(proshade_double *x1, proshade_double *y1, proshade_double *z1, proshade_double *x2, proshade_double *y2, proshade_double *z2)
Simple 3D vector dot product computation.
Definition: ProSHADE_maths.cpp:1705
ProSHADE_internal_misc::addToSignedVector
void addToSignedVector(std::vector< proshade_signed > *vecToAddTo, proshade_signed elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:121
ProSHADE_internal_maths::BicubicInterpolator
Definition: ProSHADE_maths.hpp:102
ProSHADE_internal_maths::getGLPolyAtZero
void getGLPolyAtZero(proshade_unsign order, proshade_double *polyValue, proshade_double *deriValue)
This function obtains the Legendre polynomial values and its derivative at zero for any positive inte...
Definition: ProSHADE_maths.cpp:347
ProSHADE_internal_maths::getAxisAngleFromRotationMatrix
void getAxisAngleFromRotationMatrix(proshade_double *rotMat, proshade_double *x, proshade_double *y, proshade_double *z, proshade_double *ang)
This function converts rotation matrix to the axis-angle representation.
Definition: ProSHADE_maths.cpp:1039
ProSHADE_internal_maths::advanceGLPolyValue
proshade_double advanceGLPolyValue(proshade_double from, proshade_double to, proshade_double valAtFrom, proshade_unsign noSteps, proshade_unsign taylorSeriesCap)
This function finds the next value of the polynomial.
Definition: ProSHADE_maths.cpp:477
ProSHADE_internal_misc::checkMemoryAllocation
void checkMemoryAllocation(chVar checkVar, std::string fileP, unsigned int lineP, std::string funcP, std::string infoP="This error may occurs when ProSHADE requests memory to be\n : allocated to it and this operation fails. This could\n : happen when not enough memory is available, either due to\n : other processes using a lot of memory, or when the machine\n : does not have sufficient memory available. Re-run to see\n : if this problem persists.")
Checks if memory was allocated properly.
Definition: ProSHADE_misc.hpp:65
ProSHADE_internal_maths::getEulerZXZFromSOFTPosition
void getEulerZXZFromSOFTPosition(proshade_signed band, proshade_signed x, proshade_signed y, proshade_signed z, proshade_double *eulerAlpha, proshade_double *eulerBeta, proshade_double *eulerGamma)
Function to find Euler angles (ZXZ convention) from index position in the inverse SOFT map.
Definition: ProSHADE_maths.cpp:961
ProSHADE_internal_misc::addToUnsignVector
void addToUnsignVector(std::vector< proshade_unsign > *vecToAddTo, proshade_unsign elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:99
ProSHADE_internal_maths::completeLegendreSeries
void completeLegendreSeries(proshade_unsign order, proshade_double *abscissa, proshade_double *weights, proshade_unsign taylorSeriesCap)
This function completes the Legendre polynomial series assuming you have obtained the first values.
Definition: ProSHADE_maths.cpp:521
ProSHADE_internal_maths::prepareBiCubicInterpolatorsPlusPlus
void prepareBiCubicInterpolatorsPlusPlus(proshade_double bestLattitude, proshade_double bestLongitude, std::vector< proshade_unsign > *sphereList, std::vector< ProSHADE_internal_maths::BicubicInterpolator * > *interpols, std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_sphere * > *sphereMappedRotFun)
This function prepares the interpolation objects for the bi-cubic interpolation.
Definition: ProSHADE_maths.cpp:2558
ProSHADE_internal_maths::compute3x3MatrixInverse
proshade_double * compute3x3MatrixInverse(proshade_double *mat)
Function for computing a 3x3 matrix inverse.
Definition: ProSHADE_maths.cpp:1810
ProSHADE_internal_maths::compute3x3MatrixMultiplication
proshade_double * compute3x3MatrixMultiplication(proshade_double *mat1, proshade_double *mat2)
Function for computing a 3x3 matrix multiplication.
Definition: ProSHADE_maths.cpp:1759