ProSHADE  0.7.6.1 (AUG 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 getRotationMatrixFromEulerZXZAngles (proshade_single eulerAlpha, proshade_single eulerBeta, proshade_single eulerGamma, proshade_single *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 getRotationMatrixFromAngleAxis (proshade_single *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)
 This function converts angle-axis representation to the Euler ZXZ angles representation. 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_single * compute3x3MatrixVectorMultiplication (proshade_single *mat, proshade_single x, proshade_single y, proshade_single 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...
 
proshade_double computeGaussian (proshade_double val, proshade_double sigma)
 This function computes a Gaussian (normal) distribution value given distance from mean and sigma. More...
 
std::vector< proshade_double > smoothen1D (proshade_double step, proshade_signed windowSize, proshade_double sigma, std::vector< proshade_double > data)
 This function takes a 1D vector and computes smoothened version based on the parameters. More...
 
proshade_single getResolutionOfReflection (proshade_single h, proshade_single k, proshade_single l, proshade_single xDim, proshade_single yDim, proshade_single zDim)
 This function computes the resolution of a particular reflection. More...
 
void binReciprocalSpaceReflections (proshade_unsign xInds, proshade_unsign yInds, proshade_unsign zInds, proshade_signed *noBin, proshade_signed *&binIndexing)
 This function does binning of the reciprocal space reflections. More...
 
proshade_double computeFSC (fftw_complex *fCoeffs1, fftw_complex *fCoeffs2, proshade_unsign xInds, proshade_unsign yInds, proshade_unsign zInds, proshade_signed noBins, proshade_signed *binIndexing, proshade_double **&binData, proshade_signed *&binCounts)
 This function computes the FSC. 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 479 of file ProSHADE_maths.cpp.

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

◆ 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 200 of file ProSHADE_maths.cpp.

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

◆ binReciprocalSpaceReflections()

void ProSHADE_internal_maths::binReciprocalSpaceReflections ( proshade_unsign  xInds,
proshade_unsign  yInds,
proshade_unsign  zInds,
proshade_signed *  noBin,
proshade_signed *&  binIndexing 
)

This function does binning of the reciprocal space reflections.

This funcion uses the knowledge of the cell dimensions to firstly decide which dimension is the limitting one in terms of FSC binning and then it proceeds to compute the bins, their positions in terms of distance from F000. Finally, it uses this information to compute a "mask" map, where each position has the value corresponding to the appropriate bin index, thus allowing fast binning of any map with the same dimensions as supplied to this function.

Parameters
[in]xIndsThe number of indices along the x-axis.
[in]yIndsThe number of indices along the y-axis.
[in]zIndsThe number of indices along the z-axis.
[in]noBinVariable to which the number of binds found will be saved into.
[in]binIndexingA pointer to which the map of bin belonging for each reflection will be saved into.

Definition at line 2939 of file ProSHADE_maths.cpp.

2940 {
2941  //================================================ Allocate output bin indexing memory and set to -100
2942  binIndexing = new proshade_signed [xInds * yInds * zInds];
2943  ProSHADE_internal_misc::checkMemoryAllocation ( binIndexing, __FILE__, __LINE__, __func__ );
2944  for ( size_t iter = 0; iter < static_cast< size_t > ( xInds * yInds * zInds ); iter++ ) { binIndexing[iter] = -100; }
2945  proshade_single xIndsF = static_cast< proshade_single > ( xInds ), yIndsF = static_cast< proshade_single > ( yInds ), zIndsF = static_cast< proshade_single > ( zInds );
2946 
2947  //================================================ Allocate local memory
2948  proshade_single *mins = new proshade_single[3];
2949  proshade_single *maxs = new proshade_single[3];
2950  proshade_single *resMins = new proshade_single[3];
2951  proshade_signed *resMinLoc = new proshade_signed[3];
2952  proshade_single *steps = new proshade_single[3];
2953 
2954  //================================================ Check local memory
2955  ProSHADE_internal_misc::checkMemoryAllocation ( mins, __FILE__, __LINE__, __func__ );
2956  ProSHADE_internal_misc::checkMemoryAllocation ( maxs, __FILE__, __LINE__, __func__ );
2957  ProSHADE_internal_misc::checkMemoryAllocation ( resMins, __FILE__, __LINE__, __func__ );
2958  ProSHADE_internal_misc::checkMemoryAllocation ( resMinLoc, __FILE__, __LINE__, __func__ );
2959  ProSHADE_internal_misc::checkMemoryAllocation ( steps, __FILE__, __LINE__, __func__ );
2960 
2961  //================================================ Initialise local variables
2962  proshade_single resol = 0.0f;
2963  proshade_signed reciX, reciY, reciZ, arrPos = 0, minLoc = -1;
2964  *noBin = 0;
2965 
2966  //================================================ Determine reciprocal space indexing
2967  mins[0] = std::floor ( xIndsF / -2.0f );
2968  mins[1] = std::floor ( yIndsF / -2.0f );
2969  mins[2] = std::floor ( zIndsF / -2.0f );
2970 
2971  maxs[0] = -mins[0];
2972  maxs[1] = -mins[1];
2973  maxs[2] = -mins[2];
2974 
2975  if ( xInds % 2 == 0 ) { mins[0] += 1.0f; }
2976  if ( yInds % 2 == 0 ) { mins[1] += 1.0f; }
2977  if ( zInds % 2 == 0 ) { mins[2] += 1.0f; }
2978 
2979  //================================================ Get minimum resolution based on dims for each dimension
2980  resMins[0] = ProSHADE_internal_maths::getResolutionOfReflection ( maxs[0], 0.0f, 0.0f, xIndsF, yIndsF, zIndsF );
2981  resMins[1] = ProSHADE_internal_maths::getResolutionOfReflection ( 0.0f, maxs[1], 0.0f, xIndsF, yIndsF, zIndsF );
2982  resMins[2] = ProSHADE_internal_maths::getResolutionOfReflection ( 0.0f, 0.0f, maxs[2], xIndsF, yIndsF, zIndsF );
2983 
2984  //================================================ Decide which dimension to work with (the one with the lowest resolution)
2985  resMinLoc[0] = 0; resMinLoc[1] = 0; resMinLoc[2] = 0;
2986  const FloatingPoint< proshade_single > lhs1 ( resMins[0] ), lhs2 ( resMins[1] ), lhs3 ( resMins[2] ), rhs1 ( std::min( resMins[0], std::min( resMins[1], resMins[2] ) ) );
2987  if ( lhs1.AlmostEquals ( rhs1 ) ) { resMinLoc[0] = 1; minLoc = 0; }
2988  if ( lhs2.AlmostEquals ( rhs1 ) ) { resMinLoc[1] = 1; minLoc = 1; }
2989  if ( lhs3.AlmostEquals ( rhs1 ) ) { resMinLoc[2] = 1; minLoc = 2; }
2990 
2991  //================================================ Find the bins and corresponding cut-offs
2992  std::vector< proshade_single > resArray ( static_cast< size_t > ( maxs[minLoc] - 1 ), 0.0f );
2993  std::vector< proshade_single > binArray ( static_cast< size_t > ( maxs[minLoc] - 1 ), 0.0f );
2994  for ( proshade_signed dimIt = 0; dimIt < static_cast< proshade_signed > ( maxs[minLoc] - 1 ); dimIt++ )
2995  {
2996  //============================================ Prepare steps
2997  steps[0] = ( static_cast< proshade_single > ( dimIt ) + 2.5f ) * static_cast< proshade_single > ( resMinLoc[0] );
2998  steps[1] = ( static_cast< proshade_single > ( dimIt ) + 2.5f ) * static_cast< proshade_single > ( resMinLoc[1] );
2999  steps[2] = ( static_cast< proshade_single > ( dimIt ) + 2.5f ) * static_cast< proshade_single > ( resMinLoc[2] );
3000 
3001  //============================================ Find resolution
3002  resol = ProSHADE_internal_maths::getResolutionOfReflection ( steps[0], steps[1], steps[2], xIndsF, yIndsF, zIndsF );
3003 
3004  //============================================ Assign to arrays
3005  resArray.at( static_cast< size_t > ( dimIt ) ) = resol;
3006  binArray.at( static_cast< size_t > ( dimIt ) ) = static_cast< proshade_single > ( dimIt ) + 2.5f;
3007  *noBin = dimIt + 1;
3008  }
3009 
3010  //================================================ Assign reflections to bins
3011  for ( proshade_signed xIt = 0; xIt < static_cast< proshade_signed > ( xInds ); xIt++ )
3012  {
3013  for ( proshade_signed yIt = 0; yIt < static_cast< proshade_signed > ( yInds ); yIt++ )
3014  {
3015  for ( proshade_signed zIt = 0; zIt < static_cast< proshade_signed > ( zInds / 2 ) + 1; zIt++ )
3016  {
3017  //==================================== Deal with reciprocal indices ordering
3018  reciX = xIt; if ( reciX > static_cast< proshade_signed > ( maxs[0] ) ) { reciX -= static_cast< proshade_signed > ( xInds ); }
3019  reciY = yIt; if ( reciY > static_cast< proshade_signed > ( maxs[1] ) ) { reciY -= static_cast< proshade_signed > ( yInds ); }
3020  reciZ = zIt; if ( reciZ > static_cast< proshade_signed > ( maxs[2] ) ) { reciZ -= static_cast< proshade_signed > ( zInds ); }
3021 
3022  //==================================== For each bin, check if this reflection belongs to it
3023  for ( proshade_signed binIt = 0; binIt < (*noBin); binIt++ )
3024  {
3025  //================================ Check by comparing distances
3026  if ( std::sqrt ( std::pow ( static_cast< proshade_single > ( reciX ), 2.0f ) +
3027  std::pow ( static_cast< proshade_single > ( reciY ), 2.0f ) +
3028  std::pow ( static_cast< proshade_single > ( reciZ ), 2.0f ) ) <= binArray.at( static_cast< size_t > ( binIt ) ) )
3029  {
3030  //============================ This is the bin for this reflection. Assign it.
3031  arrPos = zIt + static_cast< proshade_signed > ( zInds ) * ( yIt + static_cast< proshade_signed > ( yInds ) * xIt );
3032  binIndexing[ static_cast< size_t > ( arrPos ) ] = binIt;
3033 
3034  //============================ If one of the uneven ends, do not use Friedel's Law
3035  if ( reciX == static_cast< proshade_signed > ( mins[0] ) || -reciX == static_cast< proshade_signed > ( mins[0] ) ) { break; }
3036  if ( reciY == static_cast< proshade_signed > ( mins[1] ) || -reciY == static_cast< proshade_signed > ( mins[1] ) ) { break; }
3037  if ( reciZ == static_cast< proshade_signed > ( mins[2] ) || -reciZ == static_cast< proshade_signed > ( mins[2] ) ) { break; }
3038 
3039  //============================ Use Friedel's Law to find the second index (this is why we can use zDim / 2)
3040  reciX *= -1; if ( reciX < 0 ) { reciX += static_cast< proshade_signed > ( xInds ); }
3041  reciY *= -1; if ( reciY < 0 ) { reciY += static_cast< proshade_signed > ( yInds ); }
3042  reciZ *= -1; if ( reciZ < 0 ) { reciZ += static_cast< proshade_signed > ( zInds ); }
3043 
3044  //============================ Apply Friedel's Law
3045  arrPos = reciZ + static_cast< proshade_signed > ( zInds ) * ( reciY + static_cast< proshade_signed > ( yInds ) * reciX );
3046  binIndexing[ static_cast< size_t > ( arrPos ) ] = binIt;
3047 
3048  //============================ Done, exit bins loop
3049  break;
3050  }
3051  }
3052  }
3053  }
3054  }
3055 
3056  //================================================ Release memory
3057  delete[] mins;
3058  delete[] maxs;
3059  delete[] resMins;
3060  delete[] resMinLoc;
3061  delete[] steps;
3062 
3063  //================================================ Done
3064  return ;
3065 
3066 }

◆ 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 523 of file ProSHADE_maths.cpp.

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

◆ 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 804 of file ProSHADE_maths.cpp.

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

◆ 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 871 of file ProSHADE_maths.cpp.

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

◆ 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 1902 of file ProSHADE_maths.cpp.

1903 {
1904  //================================================ Allocate memory
1905  proshade_double* inverse = new proshade_double[9];
1906  ProSHADE_internal_misc::checkMemoryAllocation ( inverse, __FILE__, __LINE__, __func__ );
1907 
1908  //================================================ Compute determinant
1909  proshade_double matDet = ( mat[0] * mat[4] * mat[8] ) +
1910  ( mat[1] * mat[5] * mat[6] ) +
1911  ( mat[2] * mat[3] * mat[7] ) -
1912  ( mat[0] * mat[5] * mat[7] ) -
1913  ( mat[1] * mat[3] * mat[8] ) -
1914  ( mat[2] * mat[4] * mat[6] );
1915 
1916  //================================================ Compute inverse matrix
1917  inverse[0] = ( mat[4] * mat[8] - mat[5] * mat[7] ) / matDet;
1918  inverse[1] = ( mat[2] * mat[7] - mat[1] * mat[8] ) / matDet;
1919  inverse[2] = ( mat[1] * mat[5] - mat[2] * mat[4] ) / matDet;
1920  inverse[3] = ( mat[5] * mat[6] - mat[3] * mat[8] ) / matDet;
1921  inverse[4] = ( mat[0] * mat[8] - mat[2] * mat[6] ) / matDet;
1922  inverse[5] = ( mat[2] * mat[3] - mat[0] * mat[5] ) / matDet;
1923  inverse[6] = ( mat[3] * mat[7] - mat[4] * mat[6] ) / matDet;
1924  inverse[7] = ( mat[1] * mat[6] - mat[0] * mat[7] ) / matDet;
1925  inverse[8] = ( mat[0] * mat[4] - mat[1] * mat[3] ) / matDet;
1926 
1927  //================================================ Done
1928  return ( inverse );
1929 }

◆ 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 1827 of file ProSHADE_maths.cpp.

1828 {
1829  //================================================ Allocate memory
1830  proshade_double* ret = new proshade_double[9];
1831  ProSHADE_internal_misc::checkMemoryAllocation ( ret, __FILE__, __LINE__, __func__ );
1832 
1833  //================================================ Multiply
1834  ret[0] = ( mat1[0] * mat2[0] ) + ( mat1[1] * mat2[3] ) + ( mat1[2] * mat2[6] );
1835  ret[1] = ( mat1[0] * mat2[1] ) + ( mat1[1] * mat2[4] ) + ( mat1[2] * mat2[7] );
1836  ret[2] = ( mat1[0] * mat2[2] ) + ( mat1[1] * mat2[5] ) + ( mat1[2] * mat2[8] );
1837  ret[3] = ( mat1[3] * mat2[0] ) + ( mat1[4] * mat2[3] ) + ( mat1[5] * mat2[6] );
1838  ret[4] = ( mat1[3] * mat2[1] ) + ( mat1[4] * mat2[4] ) + ( mat1[5] * mat2[7] );
1839  ret[5] = ( mat1[3] * mat2[2] ) + ( mat1[4] * mat2[5] ) + ( mat1[5] * mat2[8] );
1840  ret[6] = ( mat1[6] * mat2[0] ) + ( mat1[7] * mat2[3] ) + ( mat1[8] * mat2[6] );
1841  ret[7] = ( mat1[6] * mat2[1] ) + ( mat1[7] * mat2[4] ) + ( mat1[8] * mat2[7] );
1842  ret[8] = ( mat1[6] * mat2[2] ) + ( mat1[7] * mat2[5] ) + ( mat1[8] * mat2[8] );
1843 
1844  //================================================ Done
1845  return ( ret );
1846 
1847 }

◆ compute3x3MatrixVectorMultiplication() [1/2]

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 1857 of file ProSHADE_maths.cpp.

1858 {
1859  //================================================ Allocate memory
1860  proshade_double* ret = new proshade_double[3];
1861  ProSHADE_internal_misc::checkMemoryAllocation ( ret, __FILE__, __LINE__, __func__ );
1862 
1863  //================================================ Compute the multiplication
1864  ret[0] = ( x * mat[0] ) + ( y * mat[1] ) + ( z * mat[2] );
1865  ret[1] = ( x * mat[3] ) + ( y * mat[4] ) + ( z * mat[5] );
1866  ret[2] = ( x * mat[6] ) + ( y * mat[7] ) + ( z * mat[8] );
1867 
1868  //================================================ Done
1869  return ( ret );
1870 
1871 }

◆ compute3x3MatrixVectorMultiplication() [2/2]

proshade_single * ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( proshade_single *  mat,
proshade_single  x,
proshade_single  y,
proshade_single  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 1881 of file ProSHADE_maths.cpp.

1882 {
1883  //================================================ Allocate memory
1884  proshade_single* ret = new proshade_single[3];
1885  ProSHADE_internal_misc::checkMemoryAllocation ( ret, __FILE__, __LINE__, __func__ );
1886 
1887  //================================================ Compute the multiplication
1888  ret[0] = ( x * mat[0] ) + ( y * mat[1] ) + ( z * mat[2] );
1889  ret[1] = ( x * mat[3] ) + ( y * mat[4] ) + ( z * mat[5] );
1890  ret[2] = ( x * mat[6] ) + ( y * mat[7] ) + ( z * mat[8] );
1891 
1892  //================================================ Done
1893  return ( ret );
1894 
1895 }

◆ 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 1805 of file ProSHADE_maths.cpp.

1806 {
1807  //================================================ Allocate memory
1808  proshade_double* crossProd = new proshade_double[3];
1809  ProSHADE_internal_misc::checkMemoryAllocation ( crossProd, __FILE__, __LINE__, __func__ );
1810 
1811  //================================================ Compute
1812  crossProd[0] = ( (*y1) * (*z2) ) - ( (*z1) * (*y2) );
1813  crossProd[1] = ( (*z1) * (*x2) ) - ( (*x1) * (*z2) );
1814  crossProd[2] = ( (*x1) * (*y2) ) - ( (*y1) * (*x2) );
1815 
1816  //================================================ Done
1817  return ( crossProd );
1818 
1819 }

◆ 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 1773 of file ProSHADE_maths.cpp.

1774 {
1775  //================================================ Compute and return
1776  return ( (*x1 * *x2) + (*y1 * *y2) + (*z1 * *z2) );
1777 }

◆ 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 1789 of file ProSHADE_maths.cpp.

1790 {
1791  //================================================ Compute and return
1792  return ( (x1 * x2) + (y1 * y2) + (z1 * z2) );
1793 }

◆ computeFSC()

proshade_double ProSHADE_internal_maths::computeFSC ( fftw_complex *  fCoeffs1,
fftw_complex *  fCoeffs2,
proshade_unsign  xInds,
proshade_unsign  yInds,
proshade_unsign  zInds,
proshade_signed  noBins,
proshade_signed *  binIndexing,
proshade_double **&  binData,
proshade_signed *&  binCounts 
)

This function computes the FSC.

This funcion computes the Fourier Shell Correlation weighted average from two identically sized arrays of Fourier coefficients. It requires these to have been pre-computed as well as the number of bins and bin mapping to be pre-computed (using the binReciprocalSpaceReflections function). Given all these inputs, this function simply computes all the required sums for each bin, processes them and outputs the weighted average FSC over all bins.

Parameters
[in]fCoeffs1The Fourier coefficients of the first map.
[in]fCoeffs2The Fourier coefficients of the second map.
[in]xIndsThe number of indices along the x-axis.
[in]yIndsThe number of indices along the y-axis.
[in]zIndsThe number of indices along the z-axis.
[in]noBinNumber of bins.
[in]binIndexingThe map of bin belonging for each reflection.
[in]binDataArray of arrays for holding temporary results of the FSC computation. It needs to have been already allocated and have dimensions of noBins x 12. This array is modified by the function in case the caller would like access to these.
[in]binCountsArray of counts for each bin. It needs to be pre-allocated and have dimension of noBins. This array is modified by the function in case the caller would like access to these.
[out]fscThe Fourier Shell Correlation between the two supplied Fourier coefficient maps.

Definition at line 3086 of file ProSHADE_maths.cpp.

3087 {
3088  //================================================ Initialise local variables
3089  proshade_double realOrig, realRot, imagOrig, imagRot, fsc = 0.0;;
3090  proshade_signed indx, arrPos;
3091  std::vector< proshade_double > covarByBin ( static_cast< size_t > ( noBins ), 0.0 );
3092  std::vector< proshade_double > fscByBin ( static_cast< size_t > ( noBins ), 0.0 );
3093 
3094  //================================================ Compute bin sums
3095  for ( proshade_signed xIt = 0; xIt < static_cast< proshade_signed > ( xInds ); xIt++ )
3096  {
3097  for ( proshade_signed yIt = 0; yIt < static_cast< proshade_signed > ( yInds ); yIt++ )
3098  {
3099  for ( proshade_signed zIt = 0; zIt < static_cast< proshade_signed > ( zInds ); zIt++ )
3100  {
3101  //==================================== Find array position
3102  arrPos = zIt + static_cast< proshade_signed > ( zInds ) * ( yIt + static_cast< proshade_signed > ( yInds ) * xIt );
3103 
3104  //==================================== If no bin is associated, skip this reflection
3105  indx = binIndexing[ static_cast< size_t > ( arrPos ) ];
3106  if ( ( indx < 0 ) || ( indx > noBins ) ) { continue; }
3107 
3108  //==================================== Calculate the sums
3109  realOrig = fCoeffs1[arrPos][0];
3110  imagOrig = fCoeffs1[arrPos][1];
3111  realRot = fCoeffs2[arrPos][0];
3112  imagRot = fCoeffs2[arrPos][1];
3113 
3114  binData[indx][0] += realOrig;
3115  binData[indx][1] += imagOrig;
3116  binData[indx][2] += realRot;
3117  binData[indx][3] += imagRot;
3118  binData[indx][4] += realOrig * realRot;
3119  binData[indx][5] += imagOrig * imagRot;
3120  binData[indx][6] += std::pow ( realOrig, 2.0 );
3121  binData[indx][7] += std::pow ( imagOrig, 2.0 );
3122  binData[indx][8] += std::pow ( realRot, 2.0 );
3123  binData[indx][9] += std::pow ( imagRot, 2.0 );
3124 
3125  //==================================== Update bin counts
3126  binCounts[indx] += 1;
3127  }
3128  }
3129  }
3130 
3131  //================================================ Compute covariance by bin
3132  for ( size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ )
3133  {
3134  covarByBin.at(binIt) = ( ( binData[binIt][4] + binData[binIt][5] ) / static_cast< proshade_double > ( binCounts[binIt] ) -
3135  ( ( binData[binIt][0] / static_cast< proshade_double > ( binCounts[binIt] ) *
3136  binData[binIt][2] / static_cast< proshade_double > ( binCounts[binIt] ) ) +
3137  ( binData[binIt][1] / static_cast< proshade_double > ( binCounts[binIt] ) *
3138  binData[binIt][3] / static_cast< proshade_double > ( binCounts[binIt] ) ) ) );
3139  }
3140 
3141  //================================================ Get FSC by bin
3142  for ( size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ )
3143  {
3144  binData[binIt][10] = ( binData[binIt][6] + binData[binIt][7] ) / static_cast< proshade_double > ( binCounts[binIt] ) -
3145  ( std::pow ( binData[binIt][0] / static_cast< proshade_double > ( binCounts[binIt] ), 2.0 ) +
3146  std::pow ( binData[binIt][1] / static_cast< proshade_double > ( binCounts[binIt] ), 2.0 ) );
3147  binData[binIt][11] = ( binData[binIt][8] + binData[binIt][9] ) / static_cast< proshade_double > ( binCounts[binIt] ) -
3148  ( std::pow ( binData[binIt][2] / static_cast< proshade_double > ( binCounts[binIt] ), 2.0 ) +
3149  std::pow ( binData[binIt][3] / static_cast< proshade_double > ( binCounts[binIt] ), 2.0 ) );
3150  fscByBin.at(binIt) = covarByBin.at(binIt) / ( std::sqrt ( binData[binIt][10] ) * std::sqrt ( binData[binIt][11] ) );
3151  }
3152 
3153  //================================================ Get average FSC over all bins
3154  proshade_double binSizeSum = 0.0;
3155  for ( size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ )
3156  {
3157  fsc += fscByBin.at(binIt) * static_cast< proshade_double > ( binCounts[binIt] );
3158  binSizeSum += static_cast< proshade_double > ( binCounts[binIt] );
3159  }
3160  fsc /= static_cast< proshade_double > ( binSizeSum );
3161 
3162  //================================================ Done
3163  return ( fsc );
3164 
3165 }

◆ computeGaussian()

proshade_double ProSHADE_internal_maths::computeGaussian ( proshade_double  val,
proshade_double  sigma 
)

This function computes a Gaussian (normal) distribution value given distance from mean and sigma.

This function simply returns the height of a normal distribution with a given sigma for a value specific distance from the mean.

Parameters
[in]valThe distance from the mean for which the Gaussian height should be computed.
[in]sigmaThe standard deviation of the Gaussian for which the computation is done.
[out]heightThe height of the Gaussian distribution as desctibed by the sigma.

Definition at line 2843 of file ProSHADE_maths.cpp.

2844 {
2845  //================================================ Compute cumulative probability from Z-score
2846  proshade_double zScore = ( val / sigma );
2847  proshade_double cumulativeProbability = 0.5 * std::erfc ( zScore * M_SQRT1_2 );
2848 
2849  //================================================ Symmetrise
2850  if ( cumulativeProbability > 0.5 ) { cumulativeProbability = 1.0 - cumulativeProbability; }
2851 
2852  //================================================ Done
2853  return ( cumulativeProbability );
2854 
2855 }

◆ 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 449 of file ProSHADE_maths.cpp.

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

◆ 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 2798 of file ProSHADE_maths.cpp.

2799 {
2800  //================================================ Initialise variables
2801  std::vector< proshade_unsign > ret;
2802  std::vector< std::pair< proshade_unsign, bool > > sieveOfEratosthenesArray;
2803 
2804  //================================================ Sanity check
2805  if ( upTo < 2 ) { return ( ret ); }
2806 
2807  //================================================ Initialise the sieve array up to the required number
2808  for ( proshade_unsign iter = 2; iter <= upTo; iter++ ) { sieveOfEratosthenesArray.emplace_back ( std::pair< proshade_unsign, bool > ( iter, true ) ); }
2809 
2810  //================================================ For each entry in the array
2811  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); iter++ )
2812  {
2813  //============================================ If this entry is still true
2814  if ( sieveOfEratosthenesArray.at(iter).second )
2815  {
2816  //======================================== Set all entries with the position x * [this entry value] to false
2817  for ( proshade_unsign it = iter + sieveOfEratosthenesArray.at(iter).first; it < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); it += sieveOfEratosthenesArray.at(iter).first )
2818  {
2819  sieveOfEratosthenesArray.at(it).second = false;
2820  }
2821  }
2822  }
2823 
2824  //================================================ Copy passing results to return vector
2825  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); iter++ )
2826  {
2827  if ( sieveOfEratosthenesArray.at(iter).second ) { ProSHADE_internal_misc::addToUnsignVector ( &ret, sieveOfEratosthenesArray.at(iter).first ); }
2828  }
2829 
2830  //================================================ Done
2831  return ( ret );
2832 
2833 }

◆ 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 1975 of file ProSHADE_maths.cpp.

1976 {
1977  //================================================ Allocate required memory
1978  proshade_double* inPlaneRotation = new proshade_double[9];
1979  proshade_double* basisChangeMat = new proshade_double[9];
1980  ProSHADE_internal_misc::checkMemoryAllocation ( inPlaneRotation, __FILE__, __LINE__, __func__ );
1981  ProSHADE_internal_misc::checkMemoryAllocation ( basisChangeMat, __FILE__, __LINE__, __func__ );
1982 
1983  //================================================ Normalise inputs
1984  proshade_double normF = std::sqrt( std::pow( x1, 2.0 ) + std::pow ( y1, 2.0 ) + std::pow ( z1, 2.0 ) );
1985  x1 /= normF; y1 /= normF; z1 /= normF;
1986 
1987  normF = std::sqrt( std::pow( x2, 2.0 ) + std::pow ( y2, 2.0 ) + std::pow ( z2, 2.0 ) );
1988  x2 /= normF; y2 /= normF; z2 /= normF;
1989 
1990  //================================================ Compute cross product's magnitude
1991  proshade_double* crossProd = ProSHADE_internal_maths::computeCrossProduct( &x1, &y1, &z1, &x2, &y2, &z2 );
1992  proshade_double crossProdMag = std::sqrt( std::pow( crossProd[0], 2.0 ) + std::pow ( crossProd[1], 2.0 ) + std::pow ( crossProd[2], 2.0 ) );
1993  delete[] crossProd;
1994 
1995  //================================================ Compute dot product
1996  proshade_double dotProd = ProSHADE_internal_maths::computeDotProduct ( &x1, &y1, &z1, &x2, &y2, &z2 );
1997 
1998  //================================================ Construct the in-plane rotation matrix
1999  inPlaneRotation[0] = dotProd; inPlaneRotation[1] = -crossProdMag; inPlaneRotation[2] = 0.0;
2000  inPlaneRotation[3] = crossProdMag; inPlaneRotation[4] = dotProd; inPlaneRotation[5] = 0.0;
2001  inPlaneRotation[6] = 0.0; inPlaneRotation[7] = 0.0; inPlaneRotation[8] = 1.0;
2002 
2003  //================================================ Construct change of basis matrix
2004  crossProd = ProSHADE_internal_maths::computeCrossProduct( &x2, &y2, &z2, &x1, &y1, &z1 );
2005  normF = std::sqrt ( std::pow ( x2 - ( dotProd * x1 ), 2.0 ) + std::pow ( y2 - ( dotProd * y1 ), 2.0 ) + std::pow ( z2 - ( dotProd * z1 ), 2.0 ) );
2006 
2007  basisChangeMat[0] = x1; basisChangeMat[1] = ( x2 - ( dotProd * x1 ) ) / normF; basisChangeMat[2] = crossProd[0];
2008  basisChangeMat[3] = y1; basisChangeMat[4] = ( y2 - ( dotProd * y1 ) ) / normF; basisChangeMat[5] = crossProd[1];
2009  basisChangeMat[6] = z1; basisChangeMat[7] = ( z2 - ( dotProd * z1 ) ) / normF; basisChangeMat[8] = crossProd[2];
2010 
2011  //================================================ Invert the change of basis matrix
2012  proshade_double* basisChangeMatInverse = ProSHADE_internal_maths::compute3x3MatrixInverse ( basisChangeMat );
2013 
2014  //================================================ Multiply inverse of change of basis matrix with the in plane rotation matrix, then multiply the result with the inverse
2015  proshade_double* tmpMat = ProSHADE_internal_maths::compute3x3MatrixMultiplication ( basisChangeMat, inPlaneRotation );
2016  proshade_double* rotMat = ProSHADE_internal_maths::compute3x3MatrixMultiplication ( tmpMat, basisChangeMatInverse );
2017 
2018  //================================================ Release memory
2019  delete[] crossProd;
2020  delete[] inPlaneRotation;
2021  delete[] basisChangeMat;
2022  delete[] basisChangeMatInverse;
2023  delete[] tmpMat;
2024 
2025  //================================================ Done
2026  return ( rotMat );
2027 
2028 }

◆ 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 2197 of file ProSHADE_maths.cpp.

2198 {
2199  //================================================ Initialise variables
2200  std::vector < proshade_double > ret;
2201 
2202  //================================================ Solution
2203  proshade_double solX = - ( y1 * dot2 * z3 - y1 * dot3 * z2 - z1 * dot2 * y3 + z1 * dot3 * y2 + dot1 * y3 * z2 - dot1 * z3 * y2 ) /
2204  ( -x1 * y3 * z2 + x1 * z3 * y2 + y1 * x3 * z2 - y1 * z3 * x2 - z1 * x3 * y2 + z1 * y3 * x2 );
2205  proshade_double solY = - ( x1 * dot2 * z3 - x1 * dot3 * z2 - z1 * dot2 * x3 + z1 * dot3 * x2 + dot1 * x3 * z2 - dot1 * z3 * x2 ) /
2206  ( x1 * y3 * z2 - x1 * z3 * y2 - y1 * x3 * z2 + y1 * z3 * x2 + z1 * x3 * y2 - z1 * y3 * x2 );
2207  proshade_double solZ = - ( x1 * dot2 * y3 - x1 * dot3 * y2 - y1 * dot2 * x3 + y1 * dot3 * x2 + dot1 * x3 * y2 - dot1 * y3 * x2 ) /
2208  ( -x1 * y3 * z2 + x1 * z3 * y2 + y1 * x3 * z2 - y1 * z3 * x2 - z1 * x3 * y2 + z1 * y3 * x2 );
2209 
2210  //================================================ Normalise the axis to magnitude 1
2211  proshade_double normFactor = sqrt ( pow ( solX, 2.0 ) + pow ( solY, 2.0 ) + pow ( solZ, 2.0 ) );
2212  solX /= normFactor;
2213  solY /= normFactor;
2214  solZ /= normFactor;
2215 
2216  //================================================ Set largest axis element to positive (ProSHADE standard)
2217  const FloatingPoint< proshade_double > lhs1 ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) );
2218  const FloatingPoint< proshade_double > rhs1 ( std::abs ( solX ) );
2219  const FloatingPoint< proshade_double > rhs2 ( std::abs ( solY ) );
2220  const FloatingPoint< proshade_double > rhs3 ( std::abs ( solZ ) );
2221  if ( ( ( lhs1.AlmostEquals ( rhs1 ) ) && ( solX < 0.0 ) ) ||
2222  ( ( lhs1.AlmostEquals ( rhs2 ) ) && ( solY < 0.0 ) ) ||
2223  ( ( lhs1.AlmostEquals ( rhs3 ) ) && ( solZ < 0.0 ) ) ) { solX *= -1.0; solY *= -1.0; solZ *= -1.0; }
2224 
2225  //================================================ Save solutions
2229 
2230  //================================================ Done
2231  return ( ret );
2232 
2233 }

◆ 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 2051 of file ProSHADE_maths.cpp.

2052 {
2053  //================================================ Initialise variables
2054  std::vector < proshade_double > ret;
2055 
2056  //================================================ Solution
2057  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 ) -
2058  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 ) ) *
2059  ( 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 ) ) ) -
2060  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 ) /
2061  ( 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 ) ) );
2062  proshade_double solY = ( ( dot2 * pow ( x2, 2.0 ) * pow ( z1, 3.0 ) ) /
2063  ( 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 ) -
2064  ( dot1 * pow ( x2, 2.0 ) * z2 * pow ( z1, 2.0 ) ) /
2065  ( 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 ) -
2066  ( 2.0 * x1 * dot2 * x2 * z2 * pow ( z1, 2.0 ) ) /
2067  ( 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 -
2068  ( 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 ) -
2069  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 ) *
2070  ( 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 ) /
2071  ( 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 ) ) +
2072  ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) * z1 ) /
2073  ( 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 ) +
2074  ( x1 * dot1 * x2 * pow ( y2, 2.0 ) * z1 ) /
2075  ( 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 ) +
2076  ( pow ( x1, 2.0 ) * dot2 * pow ( z2, 2.0 ) * z1 ) /
2077  ( 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 ) +
2078  ( 2.0 * x1 * dot1 * x2 * pow ( z2, 2.0 ) * z1 ) /
2079  ( 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 ) -
2080  ( y1 * dot1 * pow ( x2, 2.0 ) * y2 * z1 ) /
2081  ( 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 ) -
2082  ( x1 * y1 * dot2 * x2 * y2 * z1 ) /
2083  ( 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 +
2084  ( 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 ) -
2085  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 ) *
2086  ( 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 ) ) ) /
2087  ( 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 ) ) -
2088  ( pow ( x1, 2.0 ) * dot1 * pow ( z2, 3.0 ) ) /
2089  ( 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 ) -
2090  ( pow ( x1, 2.0 ) * dot1 * pow ( y2, 2.0 ) * z2 ) /
2091  ( 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 ) -
2092  ( x1 * pow ( y1, 2.0 ) * dot2 * x2 * z2 ) /
2093  ( 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 ) +
2094  ( pow ( x1, 2.0 ) * y1 * dot2 * y2 * z2 ) /
2095  ( 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 ) +
2096  ( x1 * y1 * dot1 * x2 * y2 * z2 ) /
2097  ( 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 );
2098  proshade_double solZ = ( - ( dot2 * pow ( x2, 2.0 ) * y2 * pow ( z1, 3.0 ) ) /
2099  ( ( 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 ) ) -
2100  ( dot2 * pow ( x2, 2.0 ) * pow ( z1, 2.0 ) ) /
2101  ( 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 ) +
2102  ( dot1 * pow ( x2, 2.0 ) * y2 * z2 * pow ( z1, 2.0 ) ) /
2103  ( ( 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 ) ) +
2104  ( 2.0 * x1 * dot2 * x2 * y2 * z2 * pow ( z1, 2.0 ) ) /
2105  ( ( 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 ) ) +
2106  ( 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 ) -
2107  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 ) *
2108  ( 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 ) /
2109  ( 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 ) ) +
2110  ( dot2 * y2 * z1 ) / ( y1 * z2 - z1 * y2 ) +
2111  ( dot1 * pow ( x2, 2.0 ) * z2 * z1 ) /
2112  ( 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 ) +
2113  ( x1 * dot2 * x2 * z2 * z1 ) /
2114  ( 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 ) -
2115  ( x1 * dot1 * x2 * pow ( y2, 3.0 ) * z1 ) /
2116  ( ( 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 ) ) +
2117  ( y1 * dot1 * pow ( x2, 2.0 ) * pow ( y2, 2.0 ) * z1 ) /
2118  ( ( 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 ) ) +
2119  ( x1 * y1 * dot2 * x2 * pow ( y2, 2.0 ) * z1 ) /
2120  ( ( 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 ) ) -
2121  ( pow ( x1, 2.0 ) * dot2 * y2 * pow ( z2, 2.0 ) * z1 ) /
2122  ( ( 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 ) ) -
2123  ( 2.0 * x1 * dot1 * x2 * y2 * pow ( z2, 2.0 ) * z1 ) /
2124  ( ( 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 ) ) -
2125  ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) * y2 * z1 ) /
2126  ( ( 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 +
2127  ( 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 ) -
2128  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 ) *
2129  ( 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 ) ) ) /
2130  ( 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 ) ) -
2131  ( 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 ) -
2132  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 ) *
2133  ( 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 ) ) ) /
2134  ( 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 ) ) -
2135  ( dot1 * y2 * z2 ) / ( y1 * z2 - z1 * y2 ) -
2136  ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) ) /
2137  ( 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 ) -
2138  ( x1 * dot1 * x2 * pow ( y2, 2.0 ) ) /
2139  ( 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 ) -
2140  ( x1 * dot1 * x2 * pow ( z2, 2.0 ) ) /
2141  ( 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 ) +
2142  ( y1 * dot1 * pow ( x2, 2.0 ) * y2 ) /
2143  ( 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 ) +
2144  ( x1 * y1 * dot2 * x2 * y2 ) /
2145  ( 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 ) +
2146  ( pow ( x1, 2.0 ) * dot1 * y2 * pow ( z2, 3.0 ) ) /
2147  ( ( 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 ) ) +
2148  ( pow ( x1, 2.0 ) * dot1 * pow ( y2, 3.0 ) * z2 ) /
2149  ( ( 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 ) ) -
2150  ( pow ( x1, 2.0 ) * y1 * dot2 * pow ( y2, 2.0 ) * z2 ) /
2151  ( ( 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 ) ) -
2152  ( x1 * y1 * dot1 * x2 * pow ( y2, 2.0 ) * z2 ) /
2153  ( ( 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 ) ) +
2154  ( x1 * pow ( y1, 2.0 ) * dot2 * x2 * y2 * z2 ) /
2155  ( ( 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;
2156 
2157  //================================================ Set largest axis element to positive (ProSHADE standard)
2158  const FloatingPoint< proshade_double > lhs1 ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) );
2159  const FloatingPoint< proshade_double > rhs1 ( std::abs ( solX ) );
2160  const FloatingPoint< proshade_double > rhs2 ( std::abs ( solY ) );
2161  const FloatingPoint< proshade_double > rhs3 ( std::abs ( solZ ) );
2162  if ( ( ( lhs1.AlmostEquals ( rhs1 ) ) && ( solX < 0.0 ) ) ||
2163  ( ( lhs1.AlmostEquals ( rhs2 ) ) && ( solY < 0.0 ) ) ||
2164  ( ( lhs1.AlmostEquals ( rhs3 ) ) && ( solZ < 0.0 ) ) ) { solX *= -1.0; solY *= -1.0; solZ *= -1.0; }
2165 
2166  //================================================ Save solutions
2170 
2171  //================================================ Done
2172  return ( ret );
2173 
2174 }

◆ 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 711 of file ProSHADE_maths.cpp.

712 {
713  //================================================ Initialise local variables
714  proshade_triplet* intData = new proshade_triplet [order];
715  ProSHADE_internal_misc::checkMemoryAllocation ( intData, __FILE__, __LINE__, __func__ );
716  proshade_triplet posVals;
717  proshade_unsign lesserPos = 0;
718  proshade_unsign upperPos = 0;
719  proshade_double lesserWeight = 0.0;
720  proshade_double upperWeight = 0.0;
721 
722  //================================================ Rescale to <order> points
723  for ( proshade_unsign absIter = 0; absIter < order; absIter++ )
724  {
725  //============================================ Init loop
726  posVals[0] = 0.0;
727  posVals[1] = 0.0;
728  posVals[2] = 0.0;
729 
730  //============================================ Find real position of abscissas
731  posVals[0] = ( ( abscissas[absIter] + 1.0 ) / 2.0 ) * integralOverRange;
732 
733 
734  //============================================ Find lesser and upper bounds
735  for ( proshade_unsign valIt = 0; valIt < valsSize; valIt++ )
736  {
737  if ( ( ( static_cast< proshade_double > ( valIt ) * maxSphereDists ) <= posVals[0] ) && ( ( ( static_cast< proshade_double > ( valIt ) + 1.0 ) * maxSphereDists ) > posVals[0] ) )
738  {
739  lesserPos = static_cast<proshade_unsign> ( valIt );
740  upperPos = static_cast<proshade_unsign> ( valIt + 1 );
741  break;
742  }
743  }
744 
745  //============================================ Linear Interpolation
746  lesserWeight = 0.0;
747  upperWeight = 0.0;
748  if ( lesserPos != 0 )
749  {
750  //======================================== 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 ...
751  lesserWeight = static_cast< proshade_double > ( upperPos ) - ( posVals[0] / maxSphereDists );
752  upperWeight = 1.0 - lesserWeight;
753 
754  posVals[1] = ( lesserWeight * vals[lesserPos-1][0] ) + ( upperWeight * vals[upperPos-1][0] );
755  posVals[2] = ( lesserWeight * vals[lesserPos-1][1] ) + ( upperWeight * vals[upperPos-1][1] );
756  }
757  else
758  {
759  //======================================== ... 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:
760  upperWeight = 1.0 - ( static_cast< proshade_double > ( upperPos ) - ( posVals[0] / maxSphereDists ) );
761 
762  posVals[1] = ( upperWeight * vals[upperPos-1][0] );
763  posVals[2] = ( upperWeight * vals[upperPos-1][1] );
764  }
765 
766  intData[absIter][0] = posVals[0];
767  intData[absIter][1] = posVals[1];
768  intData[absIter][2] = posVals[2];
769  }
770 
771  //================================================ Integrate
772  *retReal = 0.0;
773  *retImag = 0.0;
774  for ( proshade_unsign absPoint = 0; absPoint < order; absPoint++ )
775  {
776  *retReal += ( weights[absPoint] * intData[absPoint][1] );
777  *retImag += ( weights[absPoint] * intData[absPoint][2] );
778  }
779 
780  //================================================ Normalise
781  *retReal *= ( integralOverRange / 2.0 );
782  *retImag *= ( integralOverRange / 2.0 );
783 
784  //================================================ Release memory
785  delete[] intData;
786 
787  //================================================ Done
788  return ;
789 
790 }

◆ 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 621 of file ProSHADE_maths.cpp.

622 {
623  //================================================ Initialise local variables
624  proshade_double ret = 0.0;
625  proshade_complex* intData = new proshade_complex[order];
626  ProSHADE_internal_misc::checkMemoryAllocation ( intData, __FILE__, __LINE__, __func__ );
627  proshade_complex posVals;
628  proshade_unsign lesserPos = 0;
629  proshade_unsign upperPos = 0;
630  proshade_double lesserWeight = 0.0;
631  proshade_double upperWeight = 0.0;
632 
633  //================================================ Rescale to <order> points
634  for ( proshade_unsign absIter = 0; absIter < order; absIter++ )
635  {
636  //============================================ Init loop
637  posVals[0] = 0.0;
638  posVals[1] = 0.0;
639 
640  //============================================ Find real position of abscissas
641  posVals[0] = ( ( abscissas[absIter] + 1.0 ) / 2.0 ) * integralOverRange;
642 
643 
644  //============================================ Find lesser and upper bounds
645  for ( proshade_unsign valIt = 0; valIt < valsSize; valIt++ )
646  {
647  if ( ( ( static_cast< proshade_double > ( valIt ) * maxSphereDists ) <= posVals[0] ) && ( ( ( static_cast< proshade_double > ( valIt ) + 1.0 ) * maxSphereDists ) > posVals[0] ) )
648  {
649  lesserPos = static_cast<proshade_unsign> ( valIt );
650  upperPos = static_cast<proshade_unsign> ( valIt + 1 );
651  break;
652  }
653  }
654 
655  //============================================ Linear Interpolation
656  lesserWeight = 0.0;
657  upperWeight = 0.0;
658  if ( lesserPos != 0 )
659  {
660  //======================================== 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 ...
661  lesserWeight = static_cast< proshade_double > ( upperPos ) - ( posVals[0] / maxSphereDists );
662  upperWeight = 1.0 - lesserWeight;
663 
664  posVals[1] = ( lesserWeight * vals[lesserPos-1] ) + ( upperWeight * vals[upperPos-1] );
665  }
666  else
667  {
668  //======================================== ... 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:
669  upperWeight = 1.0 - ( static_cast< proshade_double > ( upperPos ) - ( posVals[0] / maxSphereDists ) );
670 
671  posVals[1] = ( upperWeight * vals[upperPos-1] );
672  }
673 
674  intData[absIter][0] = posVals[0];
675  intData[absIter][1] = posVals[1];
676  }
677 
678  //================================================ Integrate
679  for ( proshade_unsign absPoint = 0; absPoint < order; absPoint++ )
680  {
681  ret += ( weights[absPoint] * intData[absPoint][1] );
682  }
683 
684  //================================================ Normalise
685  ret *= ( integralOverRange / 2.0 );
686 
687  //================================================ Release memory
688  delete[] intData;
689 
690  //================================================ Done
691  return ( ret );
692 
693 }

◆ 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 1124 of file ProSHADE_maths.cpp.

1125 {
1126  //================================================ Initialise
1127  proshade_double singAtPiCheck = 0.01;
1128  proshade_double singAtIdentity = 0.05;
1129 
1130  //================================================ Check input for singularities
1131  if ( ( std::abs ( rotMat[1] - rotMat[3] ) < singAtPiCheck ) &&
1132  ( std::abs ( rotMat[2] - rotMat[6] ) < singAtPiCheck ) &&
1133  ( std::abs ( rotMat[5] - rotMat[7] ) < singAtPiCheck ) )
1134  {
1135  //============================================ Singularity in input! Check for identity matrix
1136  if ( ( std::abs ( rotMat[1] + rotMat[3] ) < singAtIdentity ) &&
1137  ( std::abs ( rotMat[2] + rotMat[6] ) < singAtIdentity ) &&
1138  ( std::abs ( rotMat[5] + rotMat[7] ) < singAtIdentity ) &&
1139  ( std::abs ( rotMat[0] + rotMat[4] + rotMat[8] - 3.0 ) < singAtIdentity ) )
1140  {
1141  //======================================== Identity matrix. Return 0 angle.
1142  *x = 1.0;
1143  *y = 0.0;
1144  *z = 0.0;
1145  *ang = 0.0;
1146 
1147  //======================================== Done
1148  return ;
1149  }
1150 
1151  //============================================ If we got here, this is the 180deg (pi rad) singularity. Find which axis should the rotation be done along
1152  *ang = M_PI;
1153 
1154  proshade_double xx = ( rotMat[0] + 1.0 ) / 2.0;
1155  proshade_double yy = ( rotMat[4] + 1.0 ) / 2.0;
1156  proshade_double zz = ( rotMat[8] + 1.0 ) / 2.0;
1157  proshade_double xy = ( rotMat[1] + rotMat[3] ) / 4.0;
1158  proshade_double xz = ( rotMat[2] + rotMat[6] ) / 4.0;
1159  proshade_double yz = ( rotMat[5] + rotMat[7] ) / 4.0;
1160 
1161  if ( ( xx > yy ) && ( xx > zz ) ) // XX is the largest diagonal
1162  {
1163  if ( xx < singAtPiCheck ) // and is still 0
1164  {
1165  *x = 0.0;
1166  *y = 1.0 / sqrt(2);
1167  *z = 1.0 / sqrt(2);
1168  }
1169  else
1170  {
1171  *x = sqrt ( xx );
1172  *y = xy / sqrt ( xx );
1173  *z = xz / sqrt ( xx );
1174  }
1175  }
1176 
1177  else if ( yy > zz ) // YY is the largest diagonal
1178  {
1179  if ( yy < singAtPiCheck ) // and is still 0
1180  {
1181  *x = 1.0 / sqrt(2);
1182  *y = 0.0;
1183  *z = 1.0 / sqrt(2);
1184  }
1185  else
1186  {
1187  *y = sqrt ( yy );
1188  *x = xy / sqrt ( yy );
1189  *z = yz / sqrt ( yy );
1190  }
1191  }
1192 
1193  else // ZZ is the largest diagonal
1194  {
1195  if ( zz < singAtPiCheck ) // and is still 0
1196  {
1197  *x = 1.0 / sqrt(2);
1198  *y = 1.0 / sqrt(2);
1199  *z = 0.0;
1200  }
1201  else
1202  {
1203  *z = sqrt ( zz );
1204  *x = xz / sqrt ( zz );
1205  *y = yz / sqrt ( zz );
1206  }
1207  }
1208 
1209  //============================================ Make sure largest axis is positive and so is the angle
1210  const FloatingPoint< proshade_double > lhs1 ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) );
1211  const FloatingPoint< proshade_double > rhs1 ( std::abs ( *x ) );
1212  const FloatingPoint< proshade_double > rhs2 ( std::abs ( *y ) );
1213  const FloatingPoint< proshade_double > rhs3 ( std::abs ( *z ) );
1214  if ( ( ( lhs1.AlmostEquals ( rhs1 ) ) && ( *x < 0.0 ) ) ||
1215  ( ( lhs1.AlmostEquals ( rhs2 ) ) && ( *y < 0.0 ) ) ||
1216  ( ( lhs1.AlmostEquals ( rhs3 ) ) && ( *z < 0.0 ) ) )
1217  {
1218  *x *= -1.0;
1219  *y *= -1.0;
1220  *z *= -1.0;
1221  *ang *= -1.0;
1222  }
1223  if ( *ang < 0.0 ) { *ang = ( 2.0 * M_PI ) + *ang; }
1224 
1225  //============================================ Done
1226  return ;
1227  }
1228 
1229  //================================================ No singularities! Now get angle
1230  *ang = std::acos ( ( std::max ( -1.0, std::min ( 3.0, rotMat[0] + rotMat[4] + rotMat[8] ) ) - 1.0 ) / 2.0 );
1231 
1232  //================================================ Init return values
1233  *x = 1.0;
1234  *y = 0.0;
1235  *z = 0.0;
1236 
1237  //================================================ Is angle 0? This should not happen, but will
1238  if ( std::abs ( *ang ) < singAtPiCheck )
1239  {
1240  *ang = 0.0;
1241  return ;
1242  }
1243 
1244  //================================================ Axis
1245  *x = rotMat[7] - rotMat[5];
1246  *y = rotMat[2] - rotMat[6];
1247  *z = rotMat[3] - rotMat[1];
1248 
1249  proshade_double normFactor = std::sqrt ( pow ( *x, 2.0 ) + pow ( *y, 2.0 ) + pow ( *z, 2.0 ) );
1250  *x /= normFactor;
1251  *y /= normFactor;
1252  *z /= normFactor;
1253 
1254  //================================================ Make sure largest axis is positive and so is the angle
1255  const FloatingPoint< proshade_double > lhs1 ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) );
1256  const FloatingPoint< proshade_double > rhs1 ( std::abs ( *x ) );
1257  const FloatingPoint< proshade_double > rhs2 ( std::abs ( *y ) );
1258  const FloatingPoint< proshade_double > rhs3 ( std::abs ( *z ) );
1259  if ( ( ( lhs1.AlmostEquals ( rhs1 ) ) && ( *x < 0.0 ) ) ||
1260  ( ( lhs1.AlmostEquals ( rhs2 ) ) && ( *y < 0.0 ) ) ||
1261  ( ( lhs1.AlmostEquals ( rhs3 ) ) && ( *z < 0.0 ) ) )
1262  {
1263  *x *= -1.0;
1264  *y *= -1.0;
1265  *z *= -1.0;
1266  *ang *= -1.0;
1267  }
1268  if ( *ang < 0.0 ) { *ang = ( 2.0 * M_PI ) + *ang; }
1269 
1270  //================================================ Done
1271  return ;
1272 
1273 }

◆ 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 1287 of file ProSHADE_maths.cpp.

1288 {
1289  //================================================ Initialise
1290  proshade_double singAtPiCheck = 0.01;
1291  proshade_double singAtIdentity = 0.05;
1292 
1293  //================================================ Check input for singularities
1294  if ( ( std::abs ( rotMat->at(1) - rotMat->at(3) ) < singAtPiCheck ) &&
1295  ( std::abs ( rotMat->at(2) - rotMat->at(6) ) < singAtPiCheck ) &&
1296  ( std::abs ( rotMat->at(5) - rotMat->at(7) ) < singAtPiCheck ) )
1297  {
1298  //============================================ Singularity in input! Check for identity matrix
1299  if ( ( std::abs ( rotMat->at(1) + rotMat->at(3) ) < singAtIdentity ) &&
1300  ( std::abs ( rotMat->at(2) + rotMat->at(6) ) < singAtIdentity ) &&
1301  ( std::abs ( rotMat->at(5) + rotMat->at(7) ) < singAtIdentity ) &&
1302  ( std::abs ( rotMat->at(0) + rotMat->at(4) + rotMat->at(8) - 3.0 ) < singAtIdentity ) )
1303  {
1304  //======================================== Identity matrix. Return 0 angle.
1305  *x = 1.0;
1306  *y = 0.0;
1307  *z = 0.0;
1308  *ang = 0.0;
1309 
1310  //======================================== Done
1311  return ;
1312  }
1313 
1314  //============================================ If we got here, this is the 180deg (pi rad) singularity. Find which axis should the rotation be done along
1315  *ang = M_PI;
1316 
1317  proshade_double xx = ( rotMat->at(0) + 1.0 ) / 2.0;
1318  proshade_double yy = ( rotMat->at(4) + 1.0 ) / 2.0;
1319  proshade_double zz = ( rotMat->at(8) + 1.0 ) / 2.0;
1320  proshade_double xy = ( rotMat->at(1) + rotMat->at(3) ) / 4.0;
1321  proshade_double xz = ( rotMat->at(2) + rotMat->at(6) ) / 4.0;
1322  proshade_double yz = ( rotMat->at(5) + rotMat->at(7) ) / 4.0;
1323 
1324  if ( ( xx > yy ) && ( xx > zz ) ) // XX is the largest diagonal
1325  {
1326  if ( xx < singAtPiCheck ) // and is still 0
1327  {
1328  *x = 0.0;
1329  *y = 1.0 / sqrt(2);
1330  *z = 1.0 / sqrt(2);
1331  }
1332  else
1333  {
1334  *x = sqrt ( xx );
1335  *y = xy / sqrt ( xx );
1336  *z = xz / sqrt ( xx );
1337  }
1338  }
1339 
1340  else if ( yy > zz ) // YY is the largest diagonal
1341  {
1342  if ( yy < singAtPiCheck ) // and is still 0
1343  {
1344  *x = 1.0 / sqrt(2);
1345  *y = 0.0;
1346  *z = 1.0 / sqrt(2);
1347  }
1348  else
1349  {
1350  *y = sqrt ( yy );
1351  *x = xy / sqrt ( yy );
1352  *z = yz / sqrt ( yy );
1353  }
1354  }
1355 
1356  else // ZZ is the largest diagonal
1357  {
1358  if ( zz < singAtPiCheck ) // and is still 0
1359  {
1360  *x = 1.0 / sqrt(2);
1361  *y = 1.0 / sqrt(2);
1362  *z = 0.0;
1363  }
1364  else
1365  {
1366  *z = sqrt ( zz );
1367  *x = xz / sqrt ( zz );
1368  *y = yz / sqrt ( zz );
1369  }
1370  }
1371 
1372  //============================================ Make sure largest axis is positive and so is the angle
1373  const FloatingPoint< proshade_double > lhs1 ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) );
1374  const FloatingPoint< proshade_double > rhs1 ( std::abs ( *x ) );
1375  const FloatingPoint< proshade_double > rhs2 ( std::abs ( *y ) );
1376  const FloatingPoint< proshade_double > rhs3 ( std::abs ( *z ) );
1377  if ( ( ( lhs1.AlmostEquals ( rhs1 ) ) && ( *x < 0.0 ) ) ||
1378  ( ( lhs1.AlmostEquals ( rhs2 ) ) && ( *y < 0.0 ) ) ||
1379  ( ( lhs1.AlmostEquals ( rhs3 ) ) && ( *z < 0.0 ) ) )
1380  {
1381  *x *= -1.0;
1382  *y *= -1.0;
1383  *z *= -1.0;
1384  *ang *= -1.0;
1385  }
1386 
1387  //============================================ Done
1388  return ;
1389  }
1390 
1391  //================================================ No singularities! Now get angle
1392  *ang = std::acos ( ( std::max ( -1.0, std::min ( 3.0, rotMat->at(0) + rotMat->at(4) + rotMat->at(8) ) ) - 1.0 ) / 2.0 );
1393 
1394  //================================================ Init return values
1395  *x = 1.0;
1396  *y = 0.0;
1397  *z = 0.0;
1398 
1399  //================================================ Is angle 0? This should not happen, but will
1400  if ( std::abs ( *ang ) < singAtPiCheck )
1401  {
1402  *ang = 0.0;
1403  return ;
1404  }
1405 
1406  //================================================ Axis
1407  *x = rotMat->at(7) - rotMat->at(5);
1408  *y = rotMat->at(2) - rotMat->at(6);
1409  *z = rotMat->at(3) - rotMat->at(1);
1410 
1411  proshade_double normFactor = std::sqrt ( pow ( *x, 2.0 ) + pow ( *y, 2.0 ) + pow ( *z, 2.0 ) );
1412  *x /= normFactor;
1413  *y /= normFactor;
1414  *z /= normFactor;
1415 
1416  //================================================ Make sure largest axis is positive and so is the angle
1417  const FloatingPoint< proshade_double > lhs1 ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) );
1418  const FloatingPoint< proshade_double > rhs1 ( std::abs ( *x ) );
1419  const FloatingPoint< proshade_double > rhs2 ( std::abs ( *y ) );
1420  const FloatingPoint< proshade_double > rhs3 ( std::abs ( *z ) );
1421  if ( ( ( lhs1.AlmostEquals ( rhs1 ) ) && ( *x < 0.0 ) ) ||
1422  ( ( lhs1.AlmostEquals ( rhs2 ) ) && ( *y < 0.0 ) ) ||
1423  ( ( lhs1.AlmostEquals ( rhs3 ) ) && ( *z < 0.0 ) ) )
1424  {
1425  *x *= -1.0;
1426  *y *= -1.0;
1427  *z *= -1.0;
1428  *ang *= -1.0;
1429  }
1430 
1431  //================================================ Done
1432  return ;
1433 
1434 }

◆ 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 
)

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 1598 of file ProSHADE_maths.cpp.

1599 {
1600  //================================================ If angle is 0 or infinity (anything divided by 0), return no rotation
1601  if ( ( axAng == 0.0 ) || ( std::isinf ( axAng ) ) )
1602  {
1603  //============================================ Return 0 ; 0 ; 0 for no angle
1604  *eA = 0.0;
1605  *eB = 0.0;
1606  *eG = 0.0;
1607 
1608  //============================================ Done
1609  return ;
1610  }
1611 
1612  //================================================ Compute required rotation matrix elements
1613  proshade_double cAng = std::cos ( axAng );
1614  proshade_double sAng = std::sin ( axAng );
1615  proshade_double tAng = 1.0 - cAng;
1616 
1617  proshade_double element22 = cAng + axZ * axZ * tAng;
1618 
1619  proshade_double tmp1 = axX * axZ * tAng;
1620  proshade_double tmp2 = axY * sAng;
1621  proshade_double element20 = tmp1 - tmp2;
1622  proshade_double element02 = tmp1 + tmp2;
1623 
1624  tmp1 = axY * axZ * tAng;
1625  tmp2 = axX * sAng;
1626  proshade_double element21 = tmp1 + tmp2;
1627  proshade_double element12 = tmp1 - tmp2;
1628 
1629  //================================================ Convert to Eulers
1630  if ( std::abs( element22 ) <= 0.99999 )
1631  {
1632  //============================================ This case occurs when there is no singularity in the rotation matrix (i.e. it does not have 0 or 180 degrees angle)
1633  *eA = std::atan2 ( element21, element20 );
1634  *eB = std::acos ( element22 );
1635  *eG = std::atan2 ( element12, -element02 );
1636  }
1637  else
1638  {
1639  //============================================ Compute some extra rotation matrix elements
1640  tmp1 = axX * axY * tAng;
1641  tmp2 = axZ * sAng;
1642  proshade_double element10 = tmp1 + tmp2;
1643  proshade_double element00 = cAng + axX * axX * tAng;
1644 
1645  //============================================ This case occurs when there is either 0 or 180 degrees rotation angle in the rotation matrix and therefore when beta is zero.
1646  if ( element22 >= 0.99999 )
1647  {
1648  //======================================== 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.
1649  *eA = std::atan2 ( element10, element00 );
1650  *eB = 0.0;
1651  *eG = 0.0;
1652  }
1653  if ( element22 <= -0.99999 )
1654  {
1655  //======================================== In this case, beta = PI and alpha and gamma are only defined in terms of their difference. So we arbitrarily set gamma to 0 and solve alpha.
1656  *eA = std::atan2 ( element10, element00 );
1657  *eB = M_PI;
1658  *eG = 0.0;
1659  }
1660  }
1661 
1662  //================================================ Get the angles to proper range
1663  if ( *eA < 0.0 ) { *eA = 2.0 * M_PI + *eA; }
1664  if ( *eB < 0.0 ) { *eB = M_PI + *eB; }
1665  if ( *eG < 0.0 ) { *eG = 2.0 * M_PI + *eG; }
1666 
1667  //================================================ Done
1668  return ;
1669 
1670 }

◆ 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 1547 of file ProSHADE_maths.cpp.

1548 {
1549  //================================================ Convert to Eulers
1550  if ( std::abs( rotMat[8] ) <= 0.99999 )
1551  {
1552  //============================================ This case occurs when there is no singularity in the rotation matrix (i.e. it does not have 0 or 180 degrees angle)
1553  *eA = std::atan2 ( rotMat[7], rotMat[6] );
1554  *eB = std::acos ( rotMat[8] );
1555  *eG = std::atan2 ( rotMat[5], -rotMat[2] );
1556  }
1557  else
1558  {
1559  //============================================ This case occurs when there is either 0 or 180 degrees rotation angle in the rotation matrix and therefore when beta is zero.
1560  if ( rotMat[8] >= 0.99999 )
1561  {
1562  //======================================== 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.
1563  *eA = std::atan2 ( rotMat[3], rotMat[0] );
1564  *eB = 0.0;
1565  *eG = 0.0;
1566  }
1567  if ( rotMat[8] <= -0.99999 )
1568  {
1569  //======================================== In this case, beta = PI and alpha and gamma are only defined in terms of their difference. So we arbitrarily set gamma to 0 and solve alpha.
1570  *eA = std::atan2 ( rotMat[3], rotMat[0] );
1571  *eB = M_PI;
1572  *eG = 0.0;
1573  }
1574  }
1575 
1576  //================================================ Get the angles to proper range
1577  if ( *eA < 0.0 ) { *eA = 2.0 * M_PI + *eA; }
1578  if ( *eB < 0.0 ) { *eB = M_PI + *eB; }
1579  if ( *eG < 0.0 ) { *eG = 2.0 * M_PI + *eG; }
1580 
1581  //================================================ Done
1582  return ;
1583 
1584 }

◆ 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 963 of file ProSHADE_maths.cpp.

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

◆ 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 386 of file ProSHADE_maths.cpp.

387 {
388  //================================================ Sanity check
389  if ( taylorSeriesCap < 2 )
390  {
391  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." );
392  }
393 
394  //================================================ Initialise variables
395  *abscAtZero = advanceGLPolyValue ( 0.0, -M_PI / 2.0, 0.0, order, taylorSeriesCap );
396  proshade_double hlp = 0.0;
397  proshade_double hlpVal = static_cast<proshade_double> ( order );
398  proshade_double *abscSteps;
399  proshade_double *weightSteps;
400 
401  //================================================ Allocate memory
402  abscSteps = new proshade_double [taylorSeriesCap+2];
403  weightSteps = new proshade_double [taylorSeriesCap+1];
404 
405  //================================================ Pre-set values
406  abscSteps[0] = 0.0;
407  abscSteps[1] = polyAtZero;
408  weightSteps[0] = 0.0;
409 
410  //================================================ Fill in abscissa and weight steps
411  for ( proshade_unsign iter = 0; iter <= taylorSeriesCap - 2; iter = iter + 2 )
412  {
413  hlp = static_cast<proshade_double> ( iter );
414 
415  abscSteps[iter+2] = 0.0;
416  abscSteps[iter+3] = ( hlp * ( hlp + 1.0 ) - hlpVal * ( hlpVal + 1.0 ) ) * abscSteps[iter+1] / (hlp + 1.0) / (hlp + 2.0 );
417 
418  weightSteps[iter+1] = 0.0;
419  weightSteps[iter+2] = ( hlp + 2.0 ) * abscSteps[iter+3];
420  }
421 
422  //================================================ Find abscissa and weights
423  for ( proshade_double iter = 0; iter < 5; iter++ )
424  {
425  *abscAtZero = *abscAtZero - evaluateGLSeries ( abscSteps, *abscAtZero, taylorSeriesCap ) / evaluateGLSeries ( weightSteps, *abscAtZero, taylorSeriesCap-1 );
426  }
427  *weighAtZero = evaluateGLSeries ( weightSteps, *abscAtZero, taylorSeriesCap-1 );
428 
429  //================================================ Free memory
430  delete abscSteps;
431  delete weightSteps;
432 
433  //================================================ Done
434  return ;
435 
436 }

◆ 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 349 of file ProSHADE_maths.cpp.

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

◆ 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 289 of file ProSHADE_maths.cpp.

290 {
291  //================================================ Sanity check
292  if ( order < 2 )
293  {
294  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." );
295  }
296 
297  //================================================ Initialise
298  proshade_double polyValue = 0.0;
299  proshade_double deriValue = 0.0;
300  proshade_double weightSum = 0.0;
301 
302  //================================================ Find the polynomial and derivative values at 0
303  getGLPolyAtZero ( order,
304  &polyValue,
305  &deriValue );
306 
307  //================================================ If the order is odd, then 0 is a root ...
308  if ( order % 2 == 1 )
309  {
310  abscissas[((order-1)/2)] = polyValue;
311  weights[((order-1)/2)] = deriValue;
312  }
313  else
314  {
315  // ... and if order is even, find the first root
316  getGLFirstEvenRoot ( polyValue, order, &abscissas[(order/2)], &weights[(order/2)], taylorSeriesCap );
317  }
318 
319  //================================================ Now, having computed the first roots, complete the series
320  completeLegendreSeries ( order, abscissas, weights, taylorSeriesCap );
321 
322  //================================================ Correct weights by anscissa values
323  for ( proshade_unsign iter = 0; iter < order; iter++ )
324  {
325  weights[iter] = 2.0 / ( 1.0 - abscissas[iter] ) / ( 1.0 + abscissas[iter] ) / weights[iter] / weights[iter];
326  weightSum = weightSum + weights[iter];
327  }
328 
329  //================================================ Normalise weights
330  for ( proshade_unsign iter = 0; iter < order; iter++ )
331  {
332  weights[iter] = 2.0 * weights[iter] / weightSum;
333  }
334 
335  //================================================ Done
336  return ;
337 }

◆ getResolutionOfReflection()

proshade_single ProSHADE_internal_maths::getResolutionOfReflection ( proshade_single  h,
proshade_single  k,
proshade_single  l,
proshade_single  xDim,
proshade_single  yDim,
proshade_single  zDim 
)

This function computes the resolution of a particular reflection.

Parameters
[in]hThe index of the reflection in reciprocal space along the x-axis.
[in]kThe index of the reflection in reciprocal space along the y-axis.
[in]lThe index of the reflection in reciprocal space along the z-axis.
[in]xDimThe dimension of the cell along the x-axis in Angstroms.
[in]yDimThe dimension of the cell along the y-axis in Angstroms.
[in]zDimThe dimension of the cell along the z-axis in Angstroms.
[out]retThe resolution of the particular reflection.

Definition at line 2905 of file ProSHADE_maths.cpp.

2906 {
2907  //================================================ Compute volume and proportions
2908  proshade_single vol = ( xDim * yDim * zDim );
2909  proshade_single sa = ( yDim * zDim ) / vol;
2910  proshade_single sb = ( xDim * zDim ) / vol;
2911  proshade_single sc = ( xDim * yDim ) / vol;
2912 
2913  //================================================ Compute distance
2914  proshade_single s2 = ( std::pow ( h * sa, 2.0f ) +
2915  std::pow ( k * sb, 2.0f ) +
2916  std::pow ( l * sc, 2.0f ) ) / 4.0f;
2917 
2918  //================================================ Deal with F000
2919  if ( s2 == 0.0f ) { s2 = 0.0000000001f; }
2920 
2921  //================================================ Done
2922  return ( 1.0f / ( 2.0f * std::sqrt ( s2 ) ) );
2923 
2924 }

◆ getRotationMatrixFromAngleAxis() [1/2]

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 1444 of file ProSHADE_maths.cpp.

1445 {
1446  //================================================ If angle is 0 or infinity (anything divided by 0), return identity matrix
1447  if ( ( ang == 0.0 ) || ( std::isinf ( ang ) ) )
1448  {
1449  //============================================ Create identity
1450  for ( proshade_unsign i = 0; i < 9; i++ ) { rotMat[i] = 0.0; }
1451  rotMat[0] = 1.0;
1452  rotMat[4] = 1.0;
1453  rotMat[8] = 1.0;
1454 
1455  //============================================ Done
1456  return ;
1457  }
1458 
1459  //================================================ Compute the matrix
1460  proshade_double cAng = cos ( ang );
1461  proshade_double sAng = sin ( ang );
1462  proshade_double tAng = 1.0 - cAng;
1463 
1464  rotMat[0] = cAng + x * x * tAng;
1465  rotMat[4] = cAng + y * y * tAng;
1466  rotMat[8] = cAng + z * z * tAng;
1467 
1468  proshade_double tmp1 = x * y * tAng;
1469  proshade_double tmp2 = z * sAng;
1470  rotMat[3] = tmp1 + tmp2;
1471  rotMat[1] = tmp1 - tmp2;
1472 
1473  tmp1 = x * z * tAng;
1474  tmp2 = y * sAng;
1475  rotMat[6] = tmp1 - tmp2;
1476  rotMat[2] = tmp1 + tmp2;
1477 
1478  tmp1 = y * z * tAng;
1479  tmp2 = x * sAng;
1480  rotMat[7] = tmp1 + tmp2;
1481  rotMat[5] = tmp1 - tmp2;
1482 
1483  //================================================ Done
1484  return ;
1485 
1486 }

◆ getRotationMatrixFromAngleAxis() [2/2]

void ProSHADE_internal_maths::getRotationMatrixFromAngleAxis ( proshade_single *  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 1496 of file ProSHADE_maths.cpp.

1497 {
1498  //================================================ If angle is 0 or infinity (anything divided by 0), return identity matrix
1499  if ( ( ang == 0.0 ) || ( std::isinf ( ang ) ) )
1500  {
1501  //============================================ Create identity
1502  for ( size_t i = 0; i < 9; i++ ) { rotMat[i] = 0.0f; }
1503  rotMat[0] = 1.0f;
1504  rotMat[4] = 1.0f;
1505  rotMat[8] = 1.0f;
1506 
1507  //============================================ Done
1508  return ;
1509  }
1510 
1511  //================================================ Compute the matrix
1512  proshade_single cAng = cos ( static_cast< proshade_single > ( ang ) );
1513  proshade_single sAng = sin ( static_cast< proshade_single > ( ang ) );
1514  proshade_single tAng = 1.0f - cAng;
1515 
1516  rotMat[0] = cAng + static_cast< proshade_single > ( x ) * static_cast< proshade_single > ( x ) * tAng;
1517  rotMat[4] = cAng + static_cast< proshade_single > ( y ) * static_cast< proshade_single > ( y ) * tAng;
1518  rotMat[8] = cAng + static_cast< proshade_single > ( z ) * static_cast< proshade_single > ( z ) * tAng;
1519 
1520  proshade_single tmp1 = static_cast< proshade_single > ( x ) * static_cast< proshade_single > ( y ) * tAng;
1521  proshade_single tmp2 = static_cast< proshade_single > ( z ) * sAng;
1522  rotMat[3] = tmp1 + tmp2;
1523  rotMat[1] = tmp1 - tmp2;
1524 
1525  tmp1 = static_cast< proshade_single > ( x ) * static_cast< proshade_single > ( z ) * tAng;
1526  tmp2 = static_cast< proshade_single > ( y ) * sAng;
1527  rotMat[6] = tmp1 - tmp2;
1528  rotMat[2] = tmp1 + tmp2;
1529 
1530  tmp1 = static_cast< proshade_single > ( y ) * static_cast< proshade_single > ( z ) * tAng;
1531  tmp2 = static_cast< proshade_single > ( x ) * sAng;
1532  rotMat[7] = tmp1 + tmp2;
1533  rotMat[5] = tmp1 - tmp2;
1534 
1535  //================================================ Done
1536  return ;
1537 
1538 }

◆ getRotationMatrixFromEulerZXZAngles() [1/2]

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 1007 of file ProSHADE_maths.cpp.

1008 {
1009  //================================================ No singularity/glimbal lock present
1010  if ( std::abs ( std::cos ( eulerBeta ) ) <= 0.98 )
1011  {
1012  //============================================ First row
1013  matrix[0] = cos ( eulerAlpha ) * cos ( eulerBeta ) * cos ( eulerGamma ) - sin ( eulerAlpha ) * sin ( eulerGamma );
1014  matrix[1] = sin ( eulerAlpha ) * cos ( eulerBeta ) * cos ( eulerGamma ) + cos ( eulerAlpha ) * sin ( eulerGamma );
1015  matrix[2] = -sin ( eulerBeta ) * cos ( eulerGamma );
1016 
1017  //============================================ Second row
1018  matrix[3] = -cos ( eulerAlpha ) * cos ( eulerBeta ) * sin ( eulerGamma ) - sin ( eulerAlpha ) * cos ( eulerGamma );
1019  matrix[4] = -sin ( eulerAlpha ) * cos ( eulerBeta ) * sin ( eulerGamma ) + cos ( eulerAlpha ) * cos ( eulerGamma );
1020  matrix[5] = sin ( eulerBeta ) * sin ( eulerGamma );
1021 
1022  //============================================ Third row
1023  matrix[6] = cos ( eulerAlpha ) * sin ( eulerBeta );
1024  matrix[7] = sin ( eulerAlpha ) * sin ( eulerBeta );
1025  matrix[8] = cos ( eulerBeta );
1026  }
1027  else
1028  {
1029  //============================================ Beta is either 0 or pi, making the alpha and gamma dimensions collapse into one (either only alpha+gamma or alpha-gamma are defined). In this case, we use conversion through quatermions.
1030  proshade_double qi = std::cos ( ( eulerAlpha - eulerGamma ) / 2.0 ) * std::sin ( eulerBeta / 2.0 );
1031  proshade_double qj = std::sin ( ( eulerAlpha - eulerGamma ) / 2.0 ) * std::sin ( eulerBeta / 2.0 );
1032  proshade_double qk = std::sin ( ( eulerAlpha + eulerGamma ) / 2.0 ) * std::cos ( eulerBeta / 2.0 );
1033  proshade_double qr = std::cos ( ( eulerAlpha + eulerGamma ) / 2.0 ) * std::cos ( eulerBeta / 2.0 );
1034 
1035  //============================================ First row
1036  matrix[0] = -1.0 + 2.0 * std::pow ( qi, 2.0 ) + 2.0 * std::pow ( qr, 2.0 );
1037  matrix[1] = 2.0 * ( qi * qj - qk * qr );
1038  matrix[2] = 2.0 * ( qi * qk + qj * qr );
1039 
1040  //============================================ Second row
1041  matrix[3] = 2.0 * ( qi * qj + qk * qr );
1042  matrix[4] = -1.0 + 2.0 * std::pow ( qj, 2.0 ) + 2.0 * std::pow ( qr, 2.0 );
1043  matrix[5] = 2.0 * ( qj * qk - qi * qr );
1044 
1045  //============================================ Third row
1046  matrix[6] = 2.0 * ( qi * qk - qj * qr );
1047  matrix[7] = 2.0 * ( qj * qk + qi * qr );
1048  matrix[8] = -1.0 + 2.0 * std::pow ( qk, 2.0 ) + 2.0 * std::pow ( qr, 2.0 );
1049  }
1050 
1051  //================================================ Done
1052  return ;
1053 
1054 }

◆ getRotationMatrixFromEulerZXZAngles() [2/2]

void ProSHADE_internal_maths::getRotationMatrixFromEulerZXZAngles ( proshade_single  eulerAlpha,
proshade_single  eulerBeta,
proshade_single  eulerGamma,
proshade_single *  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 1063 of file ProSHADE_maths.cpp.

1064 {
1065  //================================================ No singularity/glimbal lock present
1066  if ( std::abs ( std::cos ( eulerBeta ) ) <= 0.9999999f )
1067  {
1068  //============================================ First row
1069  matrix[0] = cos ( eulerAlpha ) * cos ( eulerBeta ) * cos ( eulerGamma ) - sin ( eulerAlpha ) * sin ( eulerGamma );
1070  matrix[1] = sin ( eulerAlpha ) * cos ( eulerBeta ) * cos ( eulerGamma ) + cos ( eulerAlpha ) * sin ( eulerGamma );
1071  matrix[2] = -sin ( eulerBeta ) * cos ( eulerGamma );
1072 
1073  //============================================ Second row
1074  matrix[3] = -cos ( eulerAlpha ) * cos ( eulerBeta ) * sin ( eulerGamma ) - sin ( eulerAlpha ) * cos ( eulerGamma );
1075  matrix[4] = -sin ( eulerAlpha ) * cos ( eulerBeta ) * sin ( eulerGamma ) + cos ( eulerAlpha ) * cos ( eulerGamma );
1076  matrix[5] = sin ( eulerBeta ) * sin ( eulerGamma );
1077 
1078  //============================================ Third row
1079  matrix[6] = cos ( eulerAlpha ) * sin ( eulerBeta );
1080  matrix[7] = sin ( eulerAlpha ) * sin ( eulerBeta );
1081  matrix[8] = cos ( eulerBeta );
1082  }
1083  else
1084  {
1085  //============================================ Beta is either 0 or pi, making the alpha and gamma dimensions collapse into one (either only alpha+gamma or alpha-gamma are defined). In this case, we use conversion through quatermions.
1086  proshade_single qi = std::cos ( ( eulerAlpha - eulerGamma ) / 2.0f ) * std::sin ( eulerBeta / 2.0f );
1087  proshade_single qj = std::sin ( ( eulerAlpha - eulerGamma ) / 2.0f ) * std::sin ( eulerBeta / 2.0f );
1088  proshade_single qk = std::sin ( ( eulerAlpha + eulerGamma ) / 2.0f ) * std::cos ( eulerBeta / 2.0f );
1089  proshade_single qr = std::cos ( ( eulerAlpha + eulerGamma ) / 2.0f ) * std::cos ( eulerBeta / 2.0f );
1090 
1091  //============================================ First row
1092  matrix[0] = -1.0f + 2.0f * std::pow ( qi, 2.0f ) + 2.0f * std::pow ( qr, 2.0f );
1093  matrix[1] = 2.0f * ( qi * qj - qk * qr );
1094  matrix[2] = 2.0f * ( qi * qk + qj * qr );
1095 
1096  //============================================ Second row
1097  matrix[3] = 2.0f * ( qi * qj + qk * qr );
1098  matrix[4] = -1.0f + 2.0f * std::pow ( qj, 2.0f ) + 2.0f * std::pow ( qr, 2.0f );
1099  matrix[5] = 2.0f * ( qj * qk - qi * qr );
1100 
1101  //============================================ Third row
1102  matrix[6] = 2.0f * ( qi * qk - qj * qr );
1103  matrix[7] = 2.0f * ( qj * qk + qi * qr );
1104  matrix[8] = -1.0f + 2.0f * std::pow ( qk, 2.0f ) + 2.0f * std::pow ( qr, 2.0f );
1105  }
1106 
1107  //================================================ Done
1108  return ;
1109 
1110 }

◆ 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 988 of file ProSHADE_maths.cpp.

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

◆ 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 2717 of file ProSHADE_maths.cpp.

2718 {
2719  //================================================ Initialise variables
2720  bool ret = true;
2721  proshade_unsign whichImprove = 0;
2722 
2723  //================================================ For each already detected member
2724  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( CSymList->size() ); grIt++ )
2725  {
2726  //============================================ Is fold the same?
2727  const FloatingPoint< proshade_double > lhs ( CSymList->at(grIt)[0] ), rhs ( axis[0] );
2728  if ( lhs.AlmostEquals ( rhs ) )
2729  {
2730  if ( ProSHADE_internal_maths::vectorOrientationSimilarity ( CSymList->at(grIt)[1], CSymList->at(grIt)[2], CSymList->at(grIt)[3], axis[1], axis[2], axis[3], tolerance ) )
2731  {
2732  ret = false;
2733  whichImprove = grIt;
2734  break;
2735  }
2736  }
2737  }
2738 
2739  //================================================ Improve, if required
2740  if ( improve && !ret )
2741  {
2742  CSymList->at(whichImprove)[1] = axis[1];
2743  CSymList->at(whichImprove)[2] = axis[2];
2744  CSymList->at(whichImprove)[3] = axis[3];
2745  CSymList->at(whichImprove)[4] = axis[4];
2746  CSymList->at(whichImprove)[5] = axis[5];
2747  }
2748 
2749  //================================================ Done
2750  return ( ret );
2751 
2752 }

◆ 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 2766 of file ProSHADE_maths.cpp.

2767 {
2768  //================================================ Initialise variables
2769  bool ret = true;
2770 
2771  //================================================ For each already detected member
2772  for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( CSymList->size() ); grIt++ )
2773  {
2774  const FloatingPoint< proshade_double > lhs ( fold ), rhs ( CSymList->at(grIt)[0] );
2775  if ( lhs.AlmostEquals ( rhs ) )
2776  {
2777  if ( ProSHADE_internal_maths::vectorOrientationSimilarity ( CSymList->at(grIt)[1], CSymList->at(grIt)[2], CSymList->at(grIt)[3], X, Y, Z, tolerance ) )
2778  {
2779  ret = false;
2780  break;
2781  }
2782  }
2783  }
2784 
2785  //================================================ Done
2786  return ( ret );
2787 
2788 }

◆ 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 2241 of file ProSHADE_maths.cpp.

2242 {
2243  //================================================ Initialise variables
2244  std::vector< proshade_double > ret;
2245 
2246  //================================================ Compute
2247  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(0) * el2->at(0) ) +
2248  ( el1->at(1) * el2->at(3) ) +
2249  ( el1->at(2) * el2->at(6) ) );
2250  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(0) * el2->at(1) ) +
2251  ( el1->at(1) * el2->at(4) ) +
2252  ( el1->at(2) * el2->at(7) ) );
2253  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(0) * el2->at(2) ) +
2254  ( el1->at(1) * el2->at(5) ) +
2255  ( el1->at(2) * el2->at(8) ) );
2256 
2257  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(3) * el2->at(0) ) +
2258  ( el1->at(4) * el2->at(3) ) +
2259  ( el1->at(5) * el2->at(6) ) );
2260  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(3) * el2->at(1) ) +
2261  ( el1->at(4) * el2->at(4) ) +
2262  ( el1->at(5) * el2->at(7) ) );
2263  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(3) * el2->at(2) ) +
2264  ( el1->at(4) * el2->at(5) ) +
2265  ( el1->at(5) * el2->at(8) ) );
2266 
2267  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(6) * el2->at(0) ) +
2268  ( el1->at(7) * el2->at(3) ) +
2269  ( el1->at(8) * el2->at(6) ) );
2270  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(6) * el2->at(1) ) +
2271  ( el1->at(7) * el2->at(4) ) +
2272  ( el1->at(8) * el2->at(7) ) );
2273  ProSHADE_internal_misc::addToDoubleVector ( &ret, ( el1->at(6) * el2->at(2) ) +
2274  ( el1->at(7) * el2->at(5) ) +
2275  ( el1->at(8) * el2->at(8) ) );
2276 
2277 
2278  //================================================ Done
2279  return ( ret );
2280 
2281 }

◆ 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 1681 of file ProSHADE_maths.cpp.

1682 {
1683  //================================================ Set res to 0.0s
1684  for ( proshade_unsign iter = 0; iter < 9; iter++ ) { res[iter] = 0.0; }
1685 
1686  //================================================ Compute the matrix multiplication
1687  for ( proshade_unsign row = 0; row < dim; row++ )
1688  {
1689  for ( proshade_unsign col = 0; col < dim; col++ )
1690  {
1691  for ( proshade_unsign inner = 0; inner < dim; inner++ )
1692  {
1693  res[(row*dim)+col] += A[(inner*dim)+row] * B[(col*dim)+inner];
1694  }
1695  }
1696  }
1697 
1698  //================================================ Done
1699  return ;
1700 
1701 }

◆ 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 1756 of file ProSHADE_maths.cpp.

1757 {
1758  //================================================ Compute and return
1759  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) ) );
1760 
1761 }

◆ 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 2396 of file ProSHADE_maths.cpp.

2397 {
2398  //================================================ Initialise variables
2399  proshade_double lonM, lonP, latM, latP, movSum;
2400  std::vector<proshade_double> latVals ( 3 );
2401  std::vector<proshade_double> lonVals ( 3 );
2402  proshade_double learningRate = 0.1;
2403  proshade_double prevVal = *bestSum;
2404  proshade_double valChange = 999.9;
2405  proshade_double origBestLat = std::round ( *bestLattitude );
2406  proshade_double origBestLon = std::round ( *bestLongitude );
2407  proshade_double tmpVal;
2408 
2409  //================================================ Initialise interpolators in all directions around the point of interest
2410  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsMinusMinus;
2411  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsMinusPlus;
2412  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsPlusMinus;
2413  std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsPlusPlus;
2414  prepareBiCubicInterpolatorsMinusMinus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsMinusMinus, sphereMappedRotFun );
2415  prepareBiCubicInterpolatorsMinusPlus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsMinusPlus, sphereMappedRotFun );
2416  prepareBiCubicInterpolatorsPlusMinus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsPlusMinus, sphereMappedRotFun );
2417  prepareBiCubicInterpolatorsPlusPlus ( std::round ( *bestLattitude ), std::round ( *bestLongitude ), sphereList, &interpolsPlusPlus, sphereMappedRotFun );
2418 
2419  //================================================ Start the pseudo gradient ascent (while there is some change)
2420  while ( valChange > 0.0001 )
2421  {
2422  //============================================ Find the surrounding points to the currently best position
2423  lonM = *bestLongitude - step;
2424  lonP = *bestLongitude + step;
2425  latM = *bestLattitude - step;
2426  latP = *bestLattitude + step;
2427 
2428  //============================================ Deal with optimising outside of prepared range - recursion
2429  const FloatingPoint< proshade_double > lhs1 ( *bestLattitude ), rhs1 ( origBestLat - 1.0 );
2430  const FloatingPoint< proshade_double > lhs2 ( *bestLattitude ), rhs2 ( origBestLat + 1.0 );
2431  const FloatingPoint< proshade_double > lhs3 ( *bestLongitude ), rhs3 ( origBestLon - 1.0 );
2432  const FloatingPoint< proshade_double > lhs4 ( *bestLongitude ), rhs4 ( origBestLon + 1.0 );
2433  if ( latM < ( origBestLat - 1.0 ) ) { tmpVal = *bestLattitude; *bestLattitude = origBestLat - 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( lhs1.AlmostEquals ( rhs1 ) ) { *bestLattitude = tmpVal; } break; }
2434  if ( latP > ( origBestLat + 1.0 ) ) { tmpVal = *bestLattitude; *bestLattitude = origBestLat + 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( lhs2.AlmostEquals ( rhs2 ) ) { *bestLattitude = tmpVal; } break; }
2435  if ( lonM < ( origBestLon - 1.0 ) ) { tmpVal = *bestLongitude; *bestLongitude = origBestLon - 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( lhs3.AlmostEquals ( rhs3 ) ) { *bestLongitude = tmpVal; } break; }
2436  if ( lonP > ( origBestLon + 1.0 ) ) { tmpVal = *bestLongitude; *bestLongitude = origBestLon + 1.0; optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step ); if ( lhs4.AlmostEquals ( rhs4 ) ) { *bestLongitude = tmpVal; } break; }
2437 
2438  //============================================ Prepare vectors of tested positions
2439  latVals.at(0) = latM; latVals.at(1) = *bestLattitude; latVals.at(2) = latP;
2440  lonVals.at(0) = lonM; lonVals.at(1) = *bestLongitude; lonVals.at(2) = lonP;
2441 
2442  //============================================ Find the best change
2443  for ( proshade_unsign laIt = 0; laIt < static_cast<proshade_unsign> ( latVals.size() ); laIt++ )
2444  {
2445  for ( proshade_unsign loIt = 0; loIt < static_cast<proshade_unsign> ( lonVals.size() ); loIt++ )
2446  {
2447  //==================================== For this combination of lat and lon, find sum over spheres
2448  movSum = 1.0;
2449  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sphereList->size() ); iter++ )
2450  {
2451  //================================ Interpolate using correct interpolators
2452  if ( ( latVals.at(laIt) <= origBestLat ) && ( lonVals.at(loIt) <= origBestLon ) ) { movSum += interpolsMinusMinus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2453  if ( ( latVals.at(laIt) <= origBestLat ) && ( lonVals.at(loIt) > origBestLon ) ) { movSum += interpolsMinusPlus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2454  if ( ( latVals.at(laIt) > origBestLat ) && ( lonVals.at(loIt) <= origBestLon ) ) { movSum += interpolsPlusMinus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2455  if ( ( latVals.at(laIt) > origBestLat ) && ( lonVals.at(loIt) > origBestLon ) ) { movSum += interpolsPlusPlus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2456  }
2457 
2458  //==================================== If position has improved, save it
2459  if ( *bestSum < movSum )
2460  {
2461  *bestSum = movSum;
2462  *bestLongitude = lonVals.at(loIt);
2463  *bestLattitude = latVals.at(laIt);
2464  }
2465  }
2466  }
2467 
2468  //============================================ Prepare for next iteration
2469  valChange = std::floor ( 100000.0 * ( *bestSum - prevVal ) ) / 100000.0;
2470  prevVal = std::floor ( 100000.0 * ( *bestSum ) ) / 100000.0;
2471  step = std::max ( ( valChange / step ) * learningRate, 0.01 );
2472  if ( learningRate >= 0.02 ) { learningRate -= 0.01; }
2473 
2474  }
2475 
2476  //================================================ Release interpolators memory
2477  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsMinusMinus.size() ); intIt++ ) { delete interpolsMinusMinus.at(intIt); }
2478  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsMinusPlus.size() ); intIt++ ) { delete interpolsMinusPlus.at(intIt); }
2479  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsPlusMinus.size() ); intIt++ ) { delete interpolsPlusMinus.at(intIt); }
2480  for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsPlusPlus.size() ); intIt++ ) { delete interpolsPlusPlus.at(intIt); }
2481 
2482  //================================================ Done
2483  return ;
2484 }

◆ 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 246 of file ProSHADE_maths.cpp.

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

◆ 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 2497 of file ProSHADE_maths.cpp.

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

◆ 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 2552 of file ProSHADE_maths.cpp.

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

◆ 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 2607 of file ProSHADE_maths.cpp.

2608 {
2609  //================================================ Initialise local variables
2610  proshade_signed latHlp, lonHlp;
2611  proshade_signed angDim = static_cast< proshade_signed > ( sphereMappedRotFun->at(0)->getAngularDim() );
2612 
2613  //================================================ Prepare the interpolator objects for interpolation around the position
2614  for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2615  {
2616  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along first dimension)
2617  proshade_double** interpGrid = new proshade_double*[4];
2618  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid, __FILE__, __LINE__, __func__ );
2619 
2620  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along second dimension)
2621  for ( proshade_unsign iter = 0; iter < 4; iter++ )
2622  {
2623  interpGrid[iter] = new proshade_double[4];
2624  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid[iter], __FILE__, __LINE__, __func__ );
2625  }
2626 
2627  //============================================ Fill in the value grid on which the interpolation is to be done
2628  for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2629  {
2630  for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2631  {
2632  latHlp = static_cast< proshade_signed > ( bestLattitude - 1 + static_cast< proshade_double > ( latIt ) ); if ( latHlp < 0 ) { latHlp += angDim; } if ( latHlp >= angDim ) { latHlp -= angDim; }
2633  lonHlp = static_cast< proshade_signed > ( bestLongitude - 2 + static_cast< proshade_double > ( lonIt ) ); if ( lonHlp < 0 ) { lonHlp += angDim; } if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2634  interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( static_cast< proshade_unsign > ( latHlp ), static_cast< proshade_unsign > ( lonHlp ) );
2635  }
2636  }
2637 
2638  //============================================ Create the interpolators
2639  ProSHADE_internal_maths::BicubicInterpolator* biCubInterp = new ProSHADE_internal_maths::BicubicInterpolator ( interpGrid, bestLattitude, bestLongitude - 1.0 );
2640  interpols->emplace_back ( biCubInterp );
2641 
2642  //============================================ Release memory
2643  for ( proshade_unsign iter = 0; iter < 4; iter++ ) { delete[] interpGrid[iter]; }
2644  delete[] interpGrid;
2645  }
2646 
2647  //================================================ Done
2648  return ;
2649 }

◆ 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 2662 of file ProSHADE_maths.cpp.

2663 {
2664  //================================================ Initialise local variables
2665  proshade_signed latHlp, lonHlp;
2666  proshade_signed angDim = static_cast< proshade_signed > ( sphereMappedRotFun->at(0)->getAngularDim() );
2667 
2668  //================================================ Prepare the interpolator objects for interpolation around the position
2669  for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2670  {
2671  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along first dimension)
2672  proshade_double** interpGrid = new proshade_double*[4];
2673  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid, __FILE__, __LINE__, __func__ );
2674 
2675  //============================================ Allocate memory for the value grid on which the interpolation is to be done (along second dimension)
2676  for ( proshade_unsign iter = 0; iter < 4; iter++ )
2677  {
2678  interpGrid[iter] = new proshade_double[4];
2679  ProSHADE_internal_misc::checkMemoryAllocation ( interpGrid[iter], __FILE__, __LINE__, __func__ );
2680  }
2681 
2682  //============================================ Fill in the value grid on which the interpolation is to be done
2683  for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2684  {
2685  for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2686  {
2687  latHlp = static_cast< proshade_signed > ( bestLattitude - 1 + static_cast< proshade_double > ( latIt ) ); if ( latHlp < 0 ) { latHlp += angDim; } if ( latHlp >= angDim ) { latHlp -= angDim; }
2688  lonHlp = static_cast< proshade_signed > ( bestLongitude - 1 + static_cast< proshade_double > ( lonIt ) ); if ( lonHlp < 0 ) { lonHlp += angDim; } if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2689  interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( static_cast< proshade_unsign > ( latHlp ), static_cast< proshade_unsign > ( lonHlp ) );
2690  }
2691  }
2692 
2693  //============================================ Create the interpolators
2694  ProSHADE_internal_maths::BicubicInterpolator* biCubInterp = new ProSHADE_internal_maths::BicubicInterpolator ( interpGrid, bestLattitude, bestLongitude );
2695  interpols->emplace_back ( biCubInterp );
2696 
2697  //============================================ Release memory
2698  for ( proshade_unsign iter = 0; iter < 4; iter++ ) { delete[] interpGrid[iter]; }
2699  delete[] interpGrid;
2700  }
2701 
2702  //================================================ Done
2703  return ;
2704 }

◆ 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 1707 of file ProSHADE_maths.cpp.

1708 {
1709  //================================================ Initialise variables
1710  std::vector < proshade_signed > ret;
1711 
1712  //================================================ Deal with negative numbers
1713  bool changeSign = false;
1714  if ( number < 0 ) { changeSign = true; number = -number; }
1715 
1716  //================================================ Deal with zero and one
1717  if ( number == 0 ) { ProSHADE_internal_misc::addToSignedVector ( &ret, 0 ); return ( ret ); }
1718  if ( number == 1 ) { ProSHADE_internal_misc::addToSignedVector ( &ret, 1 ); return ( ret ); }
1719 
1720  //================================================ Divide by 2 as long as you can
1721  while ( number % 2 == 0 )
1722  {
1724  number = number / 2;
1725  }
1726 
1727  //================================================ Check all odd numbers up to the square root
1728  for ( proshade_double posDiv = 3; posDiv <= sqrt ( static_cast< proshade_double > ( number ) ); posDiv += 2.0 )
1729  {
1730  // If posDiv is a divisor of the number, save the result
1731  while ( number % static_cast< proshade_signed > ( posDiv ) == 0 )
1732  {
1733  ProSHADE_internal_misc::addToSignedVector ( &ret, static_cast< proshade_signed > ( posDiv ) );
1734  number = number / static_cast< proshade_signed > ( posDiv );
1735  }
1736  }
1737 
1738  //================================================ If the number was a large prime number, save it as it is
1739  if ( number > 2 ) { ProSHADE_internal_misc::addToSignedVector ( &ret, number ); }
1740 
1741  //================================================ Finish dealing with negative numbers
1742  if ( changeSign ) { ret.at(0) = -ret.at(0); }
1743 
1744  //================================================ Done
1745  return ( ret );
1746 
1747 }

◆ 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 2295 of file ProSHADE_maths.cpp.

2296 {
2297  //================================================ Initialise variables
2298  bool ret = false;
2299 
2300  //================================================ Compute trace of mat1 * mat2^T
2301  proshade_double trace = ( mat1->at(0) * mat2->at(0) ) + ( mat1->at(1) * mat2->at(1) ) + ( mat1->at(2) * mat2->at(2) );
2302  trace += ( mat1->at(3) * mat2->at(3) ) + ( mat1->at(4) * mat2->at(4) ) + ( mat1->at(5) * mat2->at(5) );
2303  trace += ( mat1->at(6) * mat2->at(6) ) + ( mat1->at(7) * mat2->at(7) ) + ( mat1->at(8) * mat2->at(8) );
2304 
2305  //================================================ Subtract 3 (so that we would have 0 in case of idenity matrix)
2306  trace -= 3.0;
2307 
2308  //================================================ Compare to tolerance
2309  if ( tolerance > std::abs ( trace ) ) { ret = true; }
2310 
2311  //================================================ Done
2312  return ( ret );
2313 
2314 }

◆ smoothen1D()

std::vector< proshade_double > ProSHADE_internal_maths::smoothen1D ( proshade_double  step,
proshade_signed  windowSize,
proshade_double  sigma,
std::vector< proshade_double >  data 
)

This function takes a 1D vector and computes smoothened version based on the parameters.

This function firstly computes the Gaussian weights for each position in a window size accordingly to the parameters amd then proceeds to compute the weighted sum for each position created by sliding this window along the data. This results in smoothening of the data in accordance with the parameters.

Parameters
[in]stepThe size of the step on scale from 0.0 to 1.0 including boarders.
[in]windowSizeThe size of the averaged over window. It is assumed to be odd.
[in]sigmaThe standard deviation of the Gaussian to be used for smoothening.
[in]dataThe data to be smoothened.
[out]smoothenedA vector of smoothened values for the input data with length hist.size() - (windowSize - 1).

Definition at line 2869 of file ProSHADE_maths.cpp.

2870 {
2871  //================================================ Initialise local variables
2872  proshade_signed windowHalf = ( windowSize - 1 ) / 2;
2873  proshade_signed totSize = static_cast< proshade_signed > ( ( 1.0 / step ) + 1 );
2874  std::vector< proshade_double > smoothened ( static_cast< size_t > ( totSize - ( windowSize - 1 ) ), 0.0 );
2875  std::vector< proshade_double > winWeights ( static_cast< size_t > ( windowSize ), 0.0 );
2876 
2877  //================================================ Prepare window weights
2878  for ( proshade_double winIt = 0.0; winIt < static_cast< proshade_double > ( windowSize ); winIt += 1.0 ) { winWeights.at( static_cast< proshade_unsign > ( winIt ) ) = ProSHADE_internal_maths::computeGaussian ( ( winIt - static_cast< proshade_double > ( windowHalf ) ) * step, sigma ); }
2879 
2880  //================================================ Compute smoothened data
2881  for ( proshade_unsign it = 0; it < static_cast< proshade_unsign > ( smoothened.size() ); it++ )
2882  {
2883  //============================================ Compute window weighted average
2884  for ( proshade_signed winIt = 0; winIt < windowSize; winIt++ )
2885  {
2886  smoothened.at(it) += winWeights.at( static_cast< size_t > ( winIt ) ) * data.at( static_cast< size_t > ( static_cast< proshade_signed > ( it ) + winIt ) );
2887  }
2888  }
2889 
2890  //================================================ Done
2891  return ( smoothened );
2892 
2893 }

◆ transpose3x3MatrixInPlace()

void ProSHADE_internal_maths::transpose3x3MatrixInPlace ( proshade_double *  mat)

Transposes 3x3 matrix in place.

Parameters
[in]matThe matrix to be transposed.

Definition at line 1935 of file ProSHADE_maths.cpp.

1936 {
1937  //================================================ Initialise variables
1938  proshade_double tmp;
1939 
1940  //================================================ Transpose the non-diagonal values
1941  tmp = mat[1];
1942  mat[1] = mat[3];
1943  mat[3] = tmp;
1944 
1945  tmp = mat[2];
1946  mat[2] = mat[6];
1947  mat[6] = tmp;
1948 
1949  tmp = mat[5];
1950  mat[5] = mat[7];
1951  mat[7] = tmp;
1952 
1953  //================================================ Done
1954  return ;
1955 
1956 }

◆ 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  const FloatingPoint< proshade_double > lhs1 ( ret[0] );
132  const FloatingPoint< proshade_double > lhs2 ( ret[1] );
133  if ( !lhs1.AlmostEquals ( lhs1 ) ) { ret[0] = 0.0; }
134  if ( !lhs2.AlmostEquals ( lhs2 ) ) { ret[1] = 0.0; }
135 
136  //================================================ Return
137  return ;
138 
139 }

◆ 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 149 of file ProSHADE_maths.cpp.

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

◆ 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 2332 of file ProSHADE_maths.cpp.

2333 {
2334  //================================================ Initialise variables
2335  bool ret = false;
2336 
2337  //================================================ Cosine distance
2338  proshade_double cosDist = ( ( a1 * b1 ) + ( a2 * b2 ) + ( a3 * b3 ) ) /
2339  ( sqrt( pow( a1, 2.0 ) + pow( a2, 2.0 ) + pow( a3, 2.0 ) ) *
2340  sqrt( pow( b1, 2.0 ) + pow( b2, 2.0 ) + pow( b3, 2.0 ) ) );
2341 
2342  //================================================ Compare the absolute value of distance to 1.0 - tolerance
2343  if ( std::abs( cosDist ) > ( 1.0 - tolerance ) ) { ret = true; }
2344 
2345  //================================================ Done
2346  return ( ret );
2347 
2348 }

◆ 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 2366 of file ProSHADE_maths.cpp.

2367 {
2368  //================================================ Initialise variables
2369  bool ret = false;
2370 
2371  //================================================ Cosine distance
2372  proshade_double cosDist = ( ( a1 * b1 ) + ( a2 * b2 ) + ( a3 * b3 ) ) /
2373  ( sqrt( pow( a1, 2.0 ) + pow( a2, 2.0 ) + pow( a3, 2.0 ) ) *
2374  sqrt( pow( b1, 2.0 ) + pow( b2, 2.0 ) + pow( b3, 2.0 ) ) );
2375 
2376  //================================================ Compare the absolute value of distance to 1.0 - tolerance
2377  if ( cosDist > ( 1.0 - tolerance ) ) { ret = true; }
2378 
2379  //================================================ Done
2380  return ( ret );
2381 
2382 }
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:2552
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:1805
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:2607
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:449
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:386
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:2396
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:2332
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:2497
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::computeGaussian
proshade_double computeGaussian(proshade_double val, proshade_double sigma)
This function computes a Gaussian (normal) distribution value given distance from mean and sigma.
Definition: ProSHADE_maths.cpp:2843
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:1773
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:105
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:349
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:479
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:67
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:523
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:2662
ProSHADE_internal_maths::getResolutionOfReflection
proshade_single getResolutionOfReflection(proshade_single h, proshade_single k, proshade_single l, proshade_single xDim, proshade_single yDim, proshade_single zDim)
This function computes the resolution of a particular reflection.
Definition: ProSHADE_maths.cpp:2905
ProSHADE_internal_maths::compute3x3MatrixInverse
proshade_double * compute3x3MatrixInverse(proshade_double *mat)
Function for computing a 3x3 matrix inverse.
Definition: ProSHADE_maths.cpp:1902
ProSHADE_internal_maths::compute3x3MatrixMultiplication
proshade_double * compute3x3MatrixMultiplication(proshade_double *mat1, proshade_double *mat2)
Function for computing a 3x3 matrix multiplication.
Definition: ProSHADE_maths.cpp:1827