41 *retReal = (*r1)*(*r2) - (*i1)*(*i2);
42 *retImag = (*r1)*(*i2) + (*i1)*(*r2);
65 *retReal = (*r1)*(*r2) + (*i1)*(*i2);
66 *retImag = -(*r1)*(*i2) + (*i1)*(*r2);
86 proshade_double ret = (*r1)*(*r2) - (*i1)*(*i2);
106 proshade_double ret = (*r1)*(*r2) + (*i1)*(*i2);
124 ret[0] = std::accumulate ( vec->begin(), vec->end(), 0.0 ) /
static_cast<proshade_double
> ( vec->size() );
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 ) );
131 if ( ret[0] != ret[0] ) { ret[0] = 0.0; }
132 if ( ret[1] != ret[1] ) { ret[1] = 0.0; }
150 if ( vec->size() < 3 ) { ret[0] = 0.0; ret[1] = 0.0;
return; }
153 std::sort ( vec->begin(), vec->end() );
156 if (
static_cast<proshade_unsign
> ( vec->size() ) % 2 == 0)
158 ret[0] = ( vec->at( (
static_cast<proshade_unsign
> ( vec->size() ) / 2 ) - 1 ) +
159 vec->at(
static_cast<proshade_unsign
> ( vec->size() ) / 2 ) ) / 2.0;
163 ret[0] = vec->at(
static_cast<proshade_unsign
> ( vec->size() ) / 2 );
167 proshade_double Q1, Q3;
168 if (
static_cast<proshade_unsign
> ( vec->size() ) % 2 == 0)
170 Q1 = ( vec->at( (
static_cast<proshade_unsign
> ( vec->size() ) / 4 ) - 1 ) +
171 vec->at(
static_cast<proshade_unsign
> ( vec->size() ) / 4 ) ) / 2.0;
172 Q3 = ( vec->at( ( (
static_cast<proshade_unsign
> ( vec->size() ) / 4 ) * 3 ) - 1 ) +
173 vec->at( (
static_cast<proshade_unsign
> ( vec->size() ) / 4 ) * 3 ) ) / 2.0;
177 Q1 = vec->at(
static_cast<proshade_unsign
> ( vec->size() ) / 4 );
178 Q3 = vec->at( (
static_cast<proshade_unsign
> ( vec->size() ) / 4 ) * 3 );
201 std::sort ( vec, vec + vecSize );
204 if ( vecSize % 2 == 0)
206 ret[0] = ( vec[ ( vecSize / 2 ) - 1 ] + vec[ vecSize / 2 ] ) / 2.0;
210 ret[0] = vec[ vecSize / 2 ];
214 proshade_double Q1, Q3;
215 if ( vecSize % 2 == 0)
217 Q1 = ( vec[ ( vecSize / 4 ) - 1 ] + vec[ vecSize / 4 ] ) / 2.0;
218 Q3 = ( vec[ ( ( vecSize / 4 ) * 3 ) - 1 ] + vec[ ( vecSize / 4 ) * 3 ] ) / 2.0;
222 Q1 = vec[ vecSize / 4 ];
223 Q3 = vec[ ( vecSize / 4 ) * 3 ];
247 proshade_double xMean = 0.0;
248 proshade_double yMean = 0.0;
249 proshade_double zeroCount = 0.0;
250 for ( proshade_unsign iter = 0; iter < length; iter++ )
252 xMean += valSet1[iter];
253 yMean += valSet2[iter];
255 xMean /=
static_cast<proshade_double
> ( length - zeroCount );
256 yMean /=
static_cast<proshade_double
> ( length - zeroCount );
259 proshade_double xmmymm = 0.0;
260 proshade_double xmmsq = 0.0;
261 proshade_double ymmsq = 0.0;
262 for ( proshade_unsign iter = 0; iter < length; iter++ )
264 xmmymm += ( valSet1[iter] - xMean ) * ( valSet2[iter] - yMean );
265 xmmsq += pow( valSet1[iter] - xMean, 2.0 );
266 ymmsq += pow( valSet2[iter] - yMean, 2.0 );
269 proshade_double ret = xmmymm / ( sqrt(xmmsq) * sqrt(ymmsq) );
272 if ( std::isnan ( ret ) ) {
return ( 0.0 ); }
292 throw ProSHADE_exception (
"The integration order is too low.",
"EI00019", __FILE__, __LINE__, __func__,
"The Gauss-Legendre integration order is less than 2. This\n : seems very low; if you have a very small structure or very\n : low resolution, please manually increase the integration\n : order. Otherwise, please report this as a bug." );
296 proshade_double polyValue = 0.0;
297 proshade_double deriValue = 0.0;
298 proshade_double weightSum = 0.0;
306 if ( order % 2 == 1 )
308 abscissas[((order-1)/2)] = polyValue;
309 weights[((order-1)/2)] = deriValue;
314 getGLFirstEvenRoot ( polyValue, order, &abscissas[(order/2)], &weights[(order/2)], taylorSeriesCap );
321 for ( proshade_unsign iter = 0; iter < order; iter++ )
323 weights[iter] = 2.0 / ( 1.0 - abscissas[iter] ) / ( 1.0 + abscissas[iter] ) / weights[iter] / weights[iter];
324 weightSum = weightSum + weights[iter];
328 for ( proshade_unsign iter = 0; iter < order; iter++ )
330 weights[iter] = 2.0 * weights[iter] / weightSum;
350 proshade_double hlpVal = 0.0;
351 proshade_double prevPoly = 1.0;
352 proshade_double prevPrevPoly = 0.0;
353 proshade_double prevDeri = 0.0;
354 proshade_double prevPrevDeri = 0.0;
356 for ( proshade_unsign ordIt = 0; ordIt < order; ordIt++ )
358 hlpVal =
static_cast<proshade_double
> ( ordIt );
359 *polyValue = -hlpVal * prevPrevPoly / ( hlpVal + 1.0 );
360 *deriValue = ( ( 2.0 * hlpVal + 1.0 ) * prevPoly - hlpVal * prevPrevDeri ) / ( hlpVal + 1.0 );
361 prevPrevPoly = prevPoly;
362 prevPoly = *polyValue;
363 prevPrevDeri = prevDeri;
364 prevDeri = *deriValue;
387 if ( taylorSeriesCap < 2 )
389 throw ProSHADE_exception (
"The Taylor series cap is too low.",
"EI00020", __FILE__, __LINE__, __func__,
"The Taylor series expansion limit is less than 2. This\n : seems very low; if you have a very small structure or very\n : low resolution, please manually increase the integration\n : order. Otherwise, please report this as a bug." );
394 proshade_double hlp = 0.0;
395 proshade_double hlpVal =
static_cast<proshade_double
> ( order );
396 proshade_double *abscSteps;
397 proshade_double *weightSteps;
400 abscSteps =
new proshade_double [taylorSeriesCap+2];
401 weightSteps =
new proshade_double [taylorSeriesCap+1];
405 abscSteps[1] = polyAtZero;
406 weightSteps[0] = 0.0;
409 for ( proshade_unsign iter = 0; iter <= taylorSeriesCap - 2; iter = iter + 2 )
411 hlp =
static_cast<proshade_double
> ( iter );
413 abscSteps[iter+2] = 0.0;
414 abscSteps[iter+3] = ( hlp * ( hlp + 1.0 ) - hlpVal * ( hlpVal + 1.0 ) ) * abscSteps[iter+1] / (hlp + 1.0) / (hlp + 2.0 );
416 weightSteps[iter+1] = 0.0;
417 weightSteps[iter+2] = ( hlp + 2.0 ) * abscSteps[iter+3];
421 for ( proshade_double iter = 0; iter < 5; iter++ )
423 *abscAtZero = *abscAtZero -
evaluateGLSeries ( abscSteps, *abscAtZero, taylorSeriesCap ) /
evaluateGLSeries ( weightSteps, *abscAtZero, taylorSeriesCap-1 );
425 *weighAtZero =
evaluateGLSeries ( weightSteps, *abscAtZero, taylorSeriesCap-1 );
450 proshade_double factorialValue = 1.0;
451 proshade_double value = 0.0;
454 for ( proshade_unsign iter = 1; iter <= terms; iter++ )
456 value = value + series[iter] * factorialValue;
457 factorialValue = factorialValue * target;
480 proshade_double hlpVal = 0.0;
481 proshade_double stepSize = 0.0;
482 proshade_double valChange = 0.0;
483 proshade_double valSecChange = 0.0;
484 proshade_double squareSteps = 0.0;
485 proshade_double curVal = 0.0;
488 stepSize = ( to - from ) /
static_cast<proshade_double
> ( taylorSeriesCap );
489 squareSteps = sqrt (
static_cast<proshade_double
> ( noSteps * ( noSteps + 1 ) ) );
493 for ( proshade_unsign iter = 0; iter < taylorSeriesCap; iter++ )
495 hlpVal = ( 1.0 - valAtFrom ) * ( 1.0 + valAtFrom );
496 valChange = - stepSize * hlpVal / ( squareSteps * sqrt ( hlpVal ) - 0.5 * valAtFrom * sin ( 2.0 * curVal ) );
497 valAtFrom = valAtFrom + valChange;
499 curVal = curVal + stepSize;
501 hlpVal = ( 1.0 - valAtFrom ) * ( 1.0 + valAtFrom );
502 valSecChange = - stepSize * hlpVal / ( squareSteps * sqrt ( hlpVal ) - 0.5 * valAtFrom * sin ( 2.0 * curVal ) );
503 valAtFrom = valAtFrom + 0.5 * ( valSecChange - valChange );
524 proshade_double hlpTaylorVal = 0.0;
525 proshade_double hlpOrderVal =
static_cast<proshade_double
> ( order );
526 proshade_double abscValueChange = 0.0;
527 proshade_double prevAbsc = 0.0;
528 proshade_double *hlpAbscSeries;
529 proshade_double *hlpWeightSeries;
530 proshade_unsign noSeriesElems = 0;
531 proshade_unsign oddEvenSwitch = 0;
534 if ( order % 2 == 1 )
536 noSeriesElems = ( order - 1 ) / 2 - 1;
541 noSeriesElems = order / 2 - 1;
546 hlpAbscSeries =
new proshade_double[taylorSeriesCap+2];
547 hlpWeightSeries =
new proshade_double[taylorSeriesCap+1];
550 for ( proshade_unsign serIt = noSeriesElems + 1; serIt < order - 1; serIt++ )
553 prevAbsc = abscissas[serIt];
554 abscValueChange =
advanceGLPolyValue ( M_PI/2.0, -M_PI/2.0, prevAbsc, order, taylorSeriesCap ) - prevAbsc;
557 hlpAbscSeries[0] = 0.0;
558 hlpAbscSeries[1] = 0.0;
559 hlpAbscSeries[2] = weights[serIt];
562 hlpWeightSeries[0] = 0.0;
563 hlpWeightSeries[1] = hlpAbscSeries[2];
566 for ( proshade_unsign tayIt = 0; tayIt <= taylorSeriesCap - 2; tayIt++ )
568 hlpTaylorVal =
static_cast<proshade_double
> ( tayIt );
570 hlpAbscSeries[tayIt+3] = ( 2.0 * prevAbsc * ( hlpTaylorVal + 1.0 ) * hlpAbscSeries[tayIt+2] + ( hlpTaylorVal * ( hlpTaylorVal + 1.0 ) - hlpOrderVal *
571 ( hlpOrderVal + 1.0 ) ) * hlpAbscSeries[tayIt+1] / ( hlpTaylorVal + 1.0 ) ) / ( 1.0 - prevAbsc ) / ( 1.0 + prevAbsc ) /
572 ( hlpTaylorVal + 2.0 );
574 hlpWeightSeries[tayIt+2] = ( hlpTaylorVal + 2.0 ) * hlpAbscSeries[tayIt+3];
578 for ( proshade_unsign iter = 0; iter < 5; iter++ )
580 abscValueChange = abscValueChange -
evaluateGLSeries ( hlpAbscSeries, abscValueChange, taylorSeriesCap ) /
585 abscissas[serIt+1] = prevAbsc + abscValueChange;
586 weights[serIt+1] =
evaluateGLSeries ( hlpWeightSeries, abscValueChange, taylorSeriesCap - 1 );
589 for ( proshade_unsign serIt = 0; serIt <= noSeriesElems + oddEvenSwitch; serIt++ )
591 abscissas[serIt] = -abscissas[order-serIt-1];
592 weights[serIt] = weights[order-serIt-1];
596 delete hlpAbscSeries;
597 delete hlpWeightSeries;
622 proshade_double ret = 0.0;
623 proshade_complex* intData =
new proshade_complex[order];
625 proshade_complex posVals;
626 proshade_unsign lesserPos = 0;
627 proshade_unsign upperPos = 0;
628 proshade_double lesserWeight = 0.0;
629 proshade_double upperWeight = 0.0;
632 for ( proshade_unsign absIter = 0; absIter < order; absIter++ )
639 posVals[0] = ( ( abscissas[absIter] + 1.0 ) / 2.0 ) * integralOverRange;
643 for ( proshade_unsign valIt = 0; valIt < valsSize; valIt++ )
645 if ( ( ( valIt * maxSphereDists ) <= posVals[0] ) && ( ( ( valIt + 1 ) * maxSphereDists ) > posVals[0] ) )
647 lesserPos =
static_cast<proshade_unsign
> ( valIt );
648 upperPos =
static_cast<proshade_unsign
> ( valIt + 1 );
656 if ( lesserPos != 0 )
659 lesserWeight = upperPos - ( posVals[0] / maxSphereDists );
660 upperWeight = 1.0 - lesserWeight;
662 posVals[1] = ( lesserWeight * vals[lesserPos-1] ) + ( upperWeight * vals[upperPos-1] );
667 upperWeight = 1.0 - ( upperPos - ( posVals[0] / maxSphereDists ) );
669 posVals[1] = ( upperWeight * vals[upperPos-1] );
672 intData[absIter][0] = posVals[0];
673 intData[absIter][1] = posVals[1];
677 for ( proshade_unsign absPoint = 0; absPoint < order; absPoint++ )
679 ret += ( weights[absPoint] * intData[absPoint][1] );
683 ret *= ( integralOverRange / 2.0 );
709 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 )
712 proshade_triplet* intData =
new proshade_triplet [order];
714 proshade_triplet posVals;
715 proshade_unsign lesserPos = 0;
716 proshade_unsign upperPos = 0;
717 proshade_double lesserWeight = 0.0;
718 proshade_double upperWeight = 0.0;
721 for ( proshade_unsign absIter = 0; absIter < order; absIter++ )
729 posVals[0] = ( ( abscissas[absIter] + 1.0 ) / 2.0 ) * integralOverRange;
733 for ( proshade_unsign valIt = 0; valIt < valsSize; valIt++ )
735 if ( ( ( valIt * maxSphereDists ) <= posVals[0] ) && ( ( ( valIt + 1 ) * maxSphereDists ) > posVals[0] ) )
737 lesserPos =
static_cast<proshade_unsign
> ( valIt );
738 upperPos =
static_cast<proshade_unsign
> ( valIt + 1 );
746 if ( lesserPos != 0 )
749 lesserWeight = upperPos - ( posVals[0] / maxSphereDists );
750 upperWeight = 1.0 - lesserWeight;
752 posVals[1] = ( lesserWeight * vals[lesserPos-1][0] ) + ( upperWeight * vals[upperPos-1][0] );
753 posVals[2] = ( lesserWeight * vals[lesserPos-1][1] ) + ( upperWeight * vals[upperPos-1][1] );
758 upperWeight = 1.0 - ( upperPos - ( posVals[0] / maxSphereDists ) );
760 posVals[1] = ( upperWeight * vals[upperPos-1][0] );
761 posVals[2] = ( upperWeight * vals[upperPos-1][1] );
764 intData[absIter][0] = posVals[0];
765 intData[absIter][1] = posVals[1];
766 intData[absIter][2] = posVals[2];
772 for ( proshade_unsign absPoint = 0; absPoint < order; absPoint++ )
774 *retReal += ( weights[absPoint] * intData[absPoint][1] );
775 *retImag += ( weights[absPoint] * intData[absPoint][2] );
779 *retReal *= ( integralOverRange / 2.0 );
780 *retImag *= ( integralOverRange / 2.0 );
806 std::complex<double> *rotMatU =
new std::complex<double> [dim*dim];
807 std::complex<double> *rotMatV =
new std::complex<double> [dim*dim];
808 std::complex<double> *work =
new std::complex<double> [( 4 * dim)];
809 int workDim = ( 4 * dim);
810 double* rwork =
new double[(7 * dim)];
811 int* iwork =
new int[(8 * dim)];
822 std::complex<double> *matrixToDecompose =
new std::complex<double>[dim*dim];
824 for (
int rowIt = 0; rowIt < dim; rowIt++ )
826 for (
int colIt = 0; colIt < dim; colIt++ )
828 matrixToDecompose[(colIt*dim)+rowIt] = std::complex<double> ( mat[rowIt][colIt][0], mat[rowIt][colIt][1] );
833 zgesdd_ ( &job, &dim, &dim, matrixToDecompose, &dim, singularValues, rotMatU, &dim, rotMatV, &dim,
834 work, &workDim, rwork, iwork, &returnValue );
842 delete[] matrixToDecompose;
845 if ( returnValue != 0 )
847 throw ProSHADE_exception (
"The LAPACK complex SVD algorithm did not converge!",
"EL00021", __FILE__, __LINE__, __func__,
"LAPACK algorithm for computing the singular value\n : decomposition of complex matrices did not converge and\n : therefore it was not possible to combine SH coefficients\n : from multiple shells. Changing the resolution may help,\n : contact me if this error persists." );
873 double* singularValues =
new double[dim];
874 std::complex<double> *rotMatU =
new std::complex<double> [dim*dim];
875 std::complex<double> *rotMatV =
new std::complex<double> [dim*dim];
876 std::complex<double> *work =
new std::complex<double> [
static_cast<proshade_unsign
>( ( 3 * dim) + pow( dim, 2 ) * dim)];
877 int workDim = ( 3 * dim) + pow( dim, 2 );
878 double* rwork =
new double[
static_cast<proshade_unsign
>((5 * dim) + 5 * pow(dim,2))];
879 int* iwork =
new int[(8 * dim)];
889 std::complex<double> *matrixToDecompose =
new std::complex<double>[dim*dim];
891 for (
int rowIt = 0; rowIt < dim; rowIt++ )
893 for (
int colIt = 0; colIt < dim; colIt++ )
895 matrixToDecompose[(colIt*dim)+rowIt] = std::complex<double> ( mat[(rowIt*dim)+colIt], 0.0 );
900 zgesdd_ ( &job, &dim, &dim, matrixToDecompose, &dim, singularValues, rotMatU, &dim, rotMatV, &dim,
901 work, &workDim, rwork, iwork, &returnValue );
907 delete[] matrixToDecompose;
908 delete[] singularValues;
911 if ( ( returnValue != 0 ) && ( fail ) )
913 throw ProSHADE_exception (
"The LAPACK complex SVD algorithm did not converge!",
"EL00022", __FILE__, __LINE__, __func__,
"LAPACK algorithm for computing the singular value\n : decomposition of complex matrices did not converge and\n : therefore it was not possible to optimise the peak\n : positions in the (self-)rotation function. Changing the\n : resolution may help, contact me if this error persists." );
915 if ( ( returnValue != 0 ) && ( !fail ) )
922 for ( proshade_signed rowIt = 0; rowIt < dim; rowIt++ )
924 for ( proshade_signed colIt = 0; colIt < dim; colIt++ )
926 uAndV[(rowIt*3)+colIt] = rotMatU[( rowIt * 3 ) + colIt].real();
931 for ( proshade_signed rowIt = 0; rowIt < dim; rowIt++ )
933 for ( proshade_signed colIt = 0; colIt < dim; colIt++ )
935 uAndV[(rowIt*3)+colIt+9] = rotMatV[( rowIt * 3 ) + colIt].real();
964 *eulerGamma = ( M_PI * y /
static_cast<proshade_double
> ( 1.0 * band ) );
965 *eulerBeta = ( M_PI * x /
static_cast<proshade_double
> ( 2.0 * band ) );
966 *eulerAlpha = ( M_PI * z /
static_cast<proshade_double
> ( 1.0 * band ) );
989 *x =
static_cast<proshade_double
> ( ( eulerBeta *
static_cast<proshade_double
> ( 2.0 * band ) ) / M_PI );
990 *y =
static_cast<proshade_double
> ( ( eulerGamma *
static_cast<proshade_double
> ( band ) ) / M_PI );
991 *z =
static_cast<proshade_double
> ( ( eulerAlpha *
static_cast<proshade_double
> ( band ) ) / M_PI );
1008 matrix[0] = cos ( eulerAlpha ) * cos ( eulerBeta ) * cos ( eulerGamma ) - sin ( eulerAlpha ) * sin ( eulerGamma );
1009 matrix[1] = sin ( eulerAlpha ) * cos ( eulerBeta ) * cos ( eulerGamma ) + cos ( eulerAlpha ) * sin ( eulerGamma );
1010 matrix[2] = -sin ( eulerBeta ) * cos ( eulerGamma );
1013 matrix[3] = -cos ( eulerAlpha ) * cos ( eulerBeta ) * sin ( eulerGamma ) - sin ( eulerAlpha ) * cos ( eulerGamma );
1014 matrix[4] = -sin ( eulerAlpha ) * cos ( eulerBeta ) * sin ( eulerGamma ) + cos ( eulerAlpha ) * cos ( eulerGamma );
1015 matrix[5] = sin ( eulerBeta ) * sin ( eulerGamma );
1018 matrix[6] = cos ( eulerAlpha ) * sin ( eulerBeta );
1019 matrix[7] = sin ( eulerAlpha ) * sin ( eulerBeta );
1020 matrix[8] = cos ( eulerBeta );
1042 proshade_double singAtPiCheck = 0.01;
1043 proshade_double singAtIdentity = 0.05;
1046 if ( ( std::abs ( rotMat[1] - rotMat[3] ) < singAtPiCheck ) &&
1047 ( std::abs ( rotMat[2] - rotMat[6] ) < singAtPiCheck ) &&
1048 ( std::abs ( rotMat[5] - rotMat[7] ) < singAtPiCheck ) )
1051 if ( ( std::abs ( rotMat[1] + rotMat[3] ) < singAtIdentity ) &&
1052 ( std::abs ( rotMat[2] + rotMat[6] ) < singAtIdentity ) &&
1053 ( std::abs ( rotMat[5] + rotMat[7] ) < singAtIdentity ) &&
1054 ( std::abs ( rotMat[0] + rotMat[4] + rotMat[8] - 3.0 ) < singAtIdentity ) )
1069 proshade_double xx = ( rotMat[0] + 1.0 ) / 2.0;
1070 proshade_double yy = ( rotMat[4] + 1.0 ) / 2.0;
1071 proshade_double zz = ( rotMat[8] + 1.0 ) / 2.0;
1072 proshade_double xy = ( rotMat[1] + rotMat[3] ) / 4.0;
1073 proshade_double xz = ( rotMat[2] + rotMat[6] ) / 4.0;
1074 proshade_double yz = ( rotMat[5] + rotMat[7] ) / 4.0;
1076 if ( ( xx > yy ) && ( xx > zz ) )
1078 if ( xx < singAtPiCheck )
1087 *y = xy / sqrt ( xx );
1088 *z = xz / sqrt ( xx );
1094 if ( yy < singAtPiCheck )
1103 *x = xy / sqrt ( yy );
1104 *z = yz / sqrt ( yy );
1110 if ( zz < singAtPiCheck )
1119 *x = xz / sqrt ( zz );
1120 *y = yz / sqrt ( zz );
1125 if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1126 ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1127 ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1134 if ( *ang < 0.0 ) { *ang = ( 2.0 * M_PI ) + *ang; }
1141 *ang = std::acos ( ( std::max ( -1.0, std::min ( 3.0, rotMat[0] + rotMat[4] + rotMat[8] ) ) - 1.0 ) / 2.0 );
1149 if ( std::abs ( *ang ) < singAtPiCheck )
1156 *x = rotMat[7] - rotMat[5];
1157 *y = rotMat[2] - rotMat[6];
1158 *z = rotMat[3] - rotMat[1];
1160 proshade_double normFactor = std::sqrt ( pow ( *x, 2.0 ) + pow ( *y, 2.0 ) + pow ( *z, 2.0 ) );
1166 if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1167 ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1168 ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1175 if ( *ang < 0.0 ) { *ang = ( 2.0 * M_PI ) + *ang; }
1197 proshade_double singAtPiCheck = 0.01;
1198 proshade_double singAtIdentity = 0.05;
1201 if ( ( std::abs ( rotMat->at(1) - rotMat->at(3) ) < singAtPiCheck ) &&
1202 ( std::abs ( rotMat->at(2) - rotMat->at(6) ) < singAtPiCheck ) &&
1203 ( std::abs ( rotMat->at(5) - rotMat->at(7) ) < singAtPiCheck ) )
1206 if ( ( std::abs ( rotMat->at(1) + rotMat->at(3) ) < singAtIdentity ) &&
1207 ( std::abs ( rotMat->at(2) + rotMat->at(6) ) < singAtIdentity ) &&
1208 ( std::abs ( rotMat->at(5) + rotMat->at(7) ) < singAtIdentity ) &&
1209 ( std::abs ( rotMat->at(0) + rotMat->at(4) + rotMat->at(8) - 3.0 ) < singAtIdentity ) )
1224 proshade_double xx = ( rotMat->at(0) + 1.0 ) / 2.0;
1225 proshade_double yy = ( rotMat->at(4) + 1.0 ) / 2.0;
1226 proshade_double zz = ( rotMat->at(8) + 1.0 ) / 2.0;
1227 proshade_double xy = ( rotMat->at(1) + rotMat->at(3) ) / 4.0;
1228 proshade_double xz = ( rotMat->at(2) + rotMat->at(6) ) / 4.0;
1229 proshade_double yz = ( rotMat->at(5) + rotMat->at(7) ) / 4.0;
1231 if ( ( xx > yy ) && ( xx > zz ) )
1233 if ( xx < singAtPiCheck )
1242 *y = xy / sqrt ( xx );
1243 *z = xz / sqrt ( xx );
1249 if ( yy < singAtPiCheck )
1258 *x = xy / sqrt ( yy );
1259 *z = yz / sqrt ( yy );
1265 if ( zz < singAtPiCheck )
1274 *x = xz / sqrt ( zz );
1275 *y = yz / sqrt ( zz );
1280 if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1281 ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1282 ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1295 *ang = std::acos ( ( std::max ( -1.0, std::min ( 3.0, rotMat->at(0) + rotMat->at(4) + rotMat->at(8) ) ) - 1.0 ) / 2.0 );
1303 if ( std::abs ( *ang ) < singAtPiCheck )
1310 *x = rotMat->at(7) - rotMat->at(5);
1311 *y = rotMat->at(2) - rotMat->at(6);
1312 *z = rotMat->at(3) - rotMat->at(1);
1314 proshade_double normFactor = std::sqrt ( pow ( *x, 2.0 ) + pow ( *y, 2.0 ) + pow ( *z, 2.0 ) );
1320 if ( ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *x ) ) && ( *x < 0.0 ) ) ||
1321 ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *y ) ) && ( *y < 0.0 ) ) ||
1322 ( ( std::max ( std::abs ( *x ), std::max ( std::abs ( *y ), std::abs ( *z ) ) ) == std::abs ( *z ) ) && ( *z < 0.0 ) ) )
1346 if ( ( ang == 0.0 ) || ( std::isinf ( ang ) ) )
1349 for ( proshade_unsign i = 0; i < 9; i++ ) { rotMat[i] = 0.0; }
1359 proshade_double cAng = cos ( ang );
1360 proshade_double sAng = sin ( ang );
1361 proshade_double tAng = 1.0 - cAng;
1363 rotMat[0] = cAng + x * x * tAng;
1364 rotMat[4] = cAng + y * y * tAng;
1365 rotMat[8] = cAng + z * z * tAng;
1367 proshade_double tmp1 = x * y * tAng;
1368 proshade_double tmp2 = z * sAng;
1369 rotMat[3] = tmp1 + tmp2;
1370 rotMat[1] = tmp1 - tmp2;
1372 tmp1 = x * z * tAng;
1374 rotMat[6] = tmp1 - tmp2;
1375 rotMat[2] = tmp1 + tmp2;
1377 tmp1 = y * z * tAng;
1379 rotMat[7] = tmp1 + tmp2;
1380 rotMat[5] = tmp1 - tmp2;
1397 *eA = atan2 ( rotMat[7], rotMat[6] );
1398 *eB = acos ( rotMat[8] );
1399 *eG = atan2 ( rotMat[5], -rotMat[2] );
1402 proshade_double errLimit = 0.001;
1403 if ( ( ( rotMat[7] < errLimit ) && ( rotMat[7] > -errLimit ) ) && ( ( rotMat[6] < errLimit ) && ( rotMat[6] > -errLimit ) ) )
1409 if ( ( ( rotMat[5] < errLimit ) && ( rotMat[5] > -errLimit ) ) && ( ( rotMat[2] < errLimit ) && ( rotMat[2] > -errLimit ) ) )
1416 if ( *eA < 0.0 ) { *eA = 2.0 * M_PI + *eA; }
1417 if ( *eB < 0.0 ) { *eB = M_PI + *eB; }
1418 if ( *eG < 0.0 ) { *eG = 2.0 * M_PI + *eG; }
1440 if ( ( axAng == 0.0 ) || ( std::isinf ( axAng ) ) )
1452 proshade_double cAng = std::cos ( axAng );
1453 proshade_double sAng = std::sin ( axAng );
1454 proshade_double tAng = 1.0 - cAng;
1456 proshade_double element22 = cAng + axZ * axZ * tAng;
1458 proshade_double tmp1 = axX * axZ * tAng;
1459 proshade_double tmp2 = axY * sAng;
1460 proshade_double element20 = tmp1 - tmp2;
1461 proshade_double element02 = tmp1 + tmp2;
1463 tmp1 = axY * axZ * tAng;
1465 proshade_double element21 = tmp1 + tmp2;
1466 proshade_double element12 = tmp1 - tmp2;
1469 if ( std::abs( element22 ) <= 0.99999 )
1472 *eA = std::atan2 ( element21, element20 );
1473 *eB = std::acos ( element22 );
1474 *eG = std::atan2 ( element12, -element02 );
1479 tmp1 = axX * axY * tAng;
1481 proshade_double element10 = tmp1 + tmp2;
1482 proshade_double element00 = cAng + axX * axX * tAng;
1485 if ( element22 >= 0.99999 )
1488 *eA = std::atan2 ( element10, element00 );
1492 if ( element22 <= -0.99999 )
1495 *eA = std::atan2 ( element10, element00 );
1502 if ( *eA < 0.0 ) { *eA = 2.0 * M_PI + *eA; }
1503 if ( *eB < 0.0 ) { *eB = M_PI + *eB; }
1504 if ( *eG < 0.0 ) { *eG = 2.0 * M_PI + *eG; }
1529 proshade_double bestDist = 999.9;
1530 proshade_double eAHlp, eBHlp, eGHlp, axXHlp, axYHlp, axZHlp, axAngHlp, axDist;
1533 proshade_double* rMat =
new proshade_double[9];
1537 for ( proshade_signed xIt = 0; xIt < angDim; xIt++ )
1539 for ( proshade_signed yIt = 0; yIt < angDim; yIt++ )
1541 for ( proshade_signed zIt = 0; zIt < angDim; zIt++ )
1544 if ( bestDist < 0.001 ) {
break; }
1554 axAng = ( 2.0 * M_PI ) - axAng;
1559 if ( ( ( std::max( std::abs( axXHlp ), std::max( std::abs( axYHlp ), std::abs( axZHlp ) ) ) == std::abs( axXHlp ) ) && ( axXHlp < 0.0 ) ) ||
1560 ( ( std::max( std::abs( axXHlp ), std::max( std::abs( axYHlp ), std::abs( axZHlp ) ) ) == std::abs( axYHlp ) ) && ( axYHlp < 0.0 ) ) ||
1561 ( ( std::max( std::abs( axXHlp ), std::max( std::abs( axYHlp ), std::abs( axZHlp ) ) ) == std::abs( axZHlp ) ) && ( axZHlp < 0.0 ) ) )
1569 if ( ( ( std::max( std::abs( axX ), std::max( std::abs( axY ), std::abs( axZ ) ) ) == std::abs( axX ) ) && ( axX < 0.0 ) ) ||
1570 ( ( std::max( std::abs( axX ), std::max( std::abs( axY ), std::abs( axZ ) ) ) == std::abs( axY ) ) && ( axY < 0.0 ) ) ||
1571 ( ( std::max( std::abs( axX ), std::max( std::abs( axY ), std::abs( axZ ) ) ) == std::abs( axZ ) ) && ( axZ < 0.0 ) ) )
1580 axDist = std::abs( axAng - axAngHlp ) + ( 1.0 - std::abs ( ( ( axX * axXHlp ) + ( axY * axYHlp ) + ( axZ * axZHlp ) ) /
1581 ( sqrt( pow( axX, 2.0 ) + pow( axY, 2.0 ) + pow( axZ, 2.0 ) ) * sqrt( pow( axXHlp, 2.0 ) + pow( axYHlp, 2.0 ) + pow( axZHlp, 2.0 ) ) ) ) );
1584 if ( std::abs ( axDist ) < bestDist )
1587 bestDist = std::abs ( axDist );
1616 for ( proshade_unsign iter = 0; iter < 9; iter++ ) { res[iter] = 0.0; }
1619 for ( proshade_unsign row = 0; row < dim; row++ )
1621 for ( proshade_unsign col = 0; col < dim; col++ )
1623 for ( proshade_unsign inner = 0; inner < dim; inner++ )
1625 res[(row*dim)+col] += A[(inner*dim)+row] * B[(col*dim)+inner];
1642 std::vector < proshade_signed > ret;
1645 bool changeSign =
false;
1646 if ( number < 0 ) { changeSign =
true; number = -number; }
1653 while ( number % 2 == 0 )
1656 number = number / 2;
1660 for ( proshade_unsign posDiv = 3; posDiv <= sqrt ( number ); posDiv += 2)
1663 while ( number % posDiv == 0 )
1666 number = number / posDiv;
1674 if ( changeSign ) { ret.at(0) = -ret.at(0); }
1691 return ( ( 1.0 / sqrt ( 2.0 * M_PI * pow(standardDev,2.0) ) ) * std::exp ( - pow( value - mean, 2.0 ) / 2.0 * pow(standardDev,2.0) ) );
1708 return ( (*x1 * *x2) + (*y1 * *y2) + (*z1 * *z2) );
1724 return ( (x1 * x2) + (y1 * y2) + (z1 * z2) );
1740 proshade_double* crossProd =
new proshade_double[3];
1744 crossProd[0] = ( (*y1) * (*z2) ) - ( (*z1) * (*y2) );
1745 crossProd[1] = ( (*z1) * (*x2) ) - ( (*x1) * (*z2) );
1746 crossProd[2] = ( (*x1) * (*y2) ) - ( (*y1) * (*x2) );
1749 return ( crossProd );
1762 proshade_double* ret =
new proshade_double[9];
1766 ret[0] = ( mat1[0] * mat2[0] ) + ( mat1[1] * mat2[3] ) + ( mat1[2] * mat2[6] );
1767 ret[1] = ( mat1[0] * mat2[1] ) + ( mat1[1] * mat2[4] ) + ( mat1[2] * mat2[7] );
1768 ret[2] = ( mat1[0] * mat2[2] ) + ( mat1[1] * mat2[5] ) + ( mat1[2] * mat2[8] );
1769 ret[3] = ( mat1[3] * mat2[0] ) + ( mat1[4] * mat2[3] ) + ( mat1[5] * mat2[6] );
1770 ret[4] = ( mat1[3] * mat2[1] ) + ( mat1[4] * mat2[4] ) + ( mat1[5] * mat2[7] );
1771 ret[5] = ( mat1[3] * mat2[2] ) + ( mat1[4] * mat2[5] ) + ( mat1[5] * mat2[8] );
1772 ret[6] = ( mat1[6] * mat2[0] ) + ( mat1[7] * mat2[3] ) + ( mat1[8] * mat2[6] );
1773 ret[7] = ( mat1[6] * mat2[1] ) + ( mat1[7] * mat2[4] ) + ( mat1[8] * mat2[7] );
1774 ret[8] = ( mat1[6] * mat2[2] ) + ( mat1[7] * mat2[5] ) + ( mat1[8] * mat2[8] );
1792 proshade_double* ret =
new proshade_double[3];
1796 ret[0] = ( x * mat[0] ) + ( y * mat[1] ) + ( z * mat[2] );
1797 ret[1] = ( x * mat[3] ) + ( y * mat[4] ) + ( z * mat[5] );
1798 ret[2] = ( x * mat[6] ) + ( y * mat[7] ) + ( z * mat[8] );
1813 proshade_double* inverse =
new proshade_double[9];
1817 proshade_double matDet = ( mat[0] * mat[4] * mat[8] ) +
1818 ( mat[1] * mat[5] * mat[6] ) +
1819 ( mat[2] * mat[3] * mat[7] ) -
1820 ( mat[0] * mat[5] * mat[7] ) -
1821 ( mat[1] * mat[3] * mat[8] ) -
1822 ( mat[2] * mat[4] * mat[6] );
1825 inverse[0] = ( mat[4] * mat[8] - mat[5] * mat[7] ) / matDet;
1826 inverse[1] = ( mat[2] * mat[7] - mat[1] * mat[8] ) / matDet;
1827 inverse[2] = ( mat[1] * mat[5] - mat[2] * mat[4] ) / matDet;
1828 inverse[3] = ( mat[5] * mat[6] - mat[3] * mat[8] ) / matDet;
1829 inverse[4] = ( mat[0] * mat[8] - mat[2] * mat[6] ) / matDet;
1830 inverse[5] = ( mat[2] * mat[3] - mat[0] * mat[5] ) / matDet;
1831 inverse[6] = ( mat[3] * mat[7] - mat[4] * mat[6] ) / matDet;
1832 inverse[7] = ( mat[1] * mat[6] - mat[0] * mat[7] ) / matDet;
1833 inverse[8] = ( mat[0] * mat[4] - mat[1] * mat[3] ) / matDet;
1846 proshade_double tmp;
1886 proshade_double* inPlaneRotation =
new proshade_double[9];
1887 proshade_double* basisChangeMat =
new proshade_double[9];
1892 proshade_double normF = std::sqrt( std::pow( x1, 2.0 ) + std::pow ( y1, 2.0 ) + std::pow ( z1, 2.0 ) );
1893 x1 /= normF; y1 /= normF; z1 /= normF;
1895 normF = std::sqrt( std::pow( x2, 2.0 ) + std::pow ( y2, 2.0 ) + std::pow ( z2, 2.0 ) );
1896 x2 /= normF; y2 /= normF; z2 /= normF;
1900 proshade_double crossProdMag = std::sqrt( std::pow( crossProd[0], 2.0 ) + std::pow ( crossProd[1], 2.0 ) + std::pow ( crossProd[2], 2.0 ) );
1907 inPlaneRotation[0] = dotProd; inPlaneRotation[1] = -crossProdMag; inPlaneRotation[2] = 0.0;
1908 inPlaneRotation[3] = crossProdMag; inPlaneRotation[4] = dotProd; inPlaneRotation[5] = 0.0;
1909 inPlaneRotation[6] = 0.0; inPlaneRotation[7] = 0.0; inPlaneRotation[8] = 1.0;
1913 normF = std::sqrt ( std::pow ( x2 - ( dotProd * x1 ), 2.0 ) + std::pow ( y2 - ( dotProd * y1 ), 2.0 ) + std::pow ( z2 - ( dotProd * z1 ), 2.0 ) );
1915 basisChangeMat[0] = x1; basisChangeMat[1] = ( x2 - ( dotProd * x1 ) ) / normF; basisChangeMat[2] = crossProd[0];
1916 basisChangeMat[3] = y1; basisChangeMat[4] = ( y2 - ( dotProd * y1 ) ) / normF; basisChangeMat[5] = crossProd[1];
1917 basisChangeMat[6] = z1; basisChangeMat[7] = ( z2 - ( dotProd * z1 ) ) / normF; basisChangeMat[8] = crossProd[2];
1928 delete[] inPlaneRotation;
1929 delete[] basisChangeMat;
1930 delete[] basisChangeMatInverse;
1962 std::vector < proshade_double > ret;
1965 proshade_double solX = ( -sqrt ( pow ( 2.0 * x1 * y1 * dot2 * y2 + 2.0 * x1 * z1 * dot2 * z2 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( y1, 2.0 ) * dot2 * x2 + 2.0 * y1 * dot1 * x2 * y2 - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
1966 4.0 * ( pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * x1 * y1 * x2 * y2 - 2.0 * x1 * z1 * x2 * z2 + pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) ) *
1967 ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 - 2.0 * y1 * dot1 * dot2 * y2 + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) - 2.0 * z1 * dot1 * dot2 * z2 + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) ) ) -
1968 2.0 * x1 * y1 * dot2 * y2 - 2.0 * x1 * z1 * dot2 * z2 + 2.0 * x1 * dot1 * pow ( y2, 2.0 ) + 2.0 * x1 * dot1 * pow ( z2, 2.0 ) + 2.0 * pow ( y1, 2.0 ) * dot2 * x2 - 2.0 * y1 * dot1 * x2 * y2 + 2.0 * pow ( z1, 2.0 ) * dot2 * x2 - 2.0 * z1 * dot1 * x2 * z2 ) /
1969 ( 2.0 * ( pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * x1 * y1 * x2 * y2 - 2.0 * x1 * z1 * x2 * z2 + pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) ) );
1970 proshade_double solY = ( ( dot2 * pow ( x2, 2.0 ) * pow ( z1, 3.0 ) ) /
1971 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1972 ( dot1 * pow ( x2, 2.0 ) * z2 * pow ( z1, 2.0 ) ) /
1973 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1974 ( 2.0 * x1 * dot2 * x2 * z2 * pow ( z1, 2.0 ) ) /
1975 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) - dot2 * z1 -
1976 ( x2 * sqrt ( pow ( - 2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
1977 4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
1978 ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) * z1 ) /
1979 ( 2.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
1980 ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) * z1 ) /
1981 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
1982 ( x1 * dot1 * x2 * pow ( y2, 2.0 ) * z1 ) /
1983 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
1984 ( pow ( x1, 2.0 ) * dot2 * pow ( z2, 2.0 ) * z1 ) /
1985 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
1986 ( 2.0 * x1 * dot1 * x2 * pow ( z2, 2.0 ) * z1 ) /
1987 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1988 ( y1 * dot1 * pow ( x2, 2.0 ) * y2 * z1 ) /
1989 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1990 ( x1 * y1 * dot2 * x2 * y2 * z1 ) /
1991 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) + dot1 * z2 +
1992 ( x1 * z2 * sqrt ( pow ( -2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
1993 4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
1994 ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) ) /
1995 ( 2.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
1996 ( pow ( x1, 2.0 ) * dot1 * pow ( z2, 3.0 ) ) /
1997 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
1998 ( pow ( x1, 2.0 ) * dot1 * pow ( y2, 2.0 ) * z2 ) /
1999 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2000 ( x1 * pow ( y1, 2.0 ) * dot2 * x2 * z2 ) /
2001 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2002 ( pow ( x1, 2.0 ) * y1 * dot2 * y2 * z2 ) /
2003 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2004 ( x1 * y1 * dot1 * x2 * y2 * z2 ) /
2005 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) / ( y1 * z2 - z1 * y2 );
2006 proshade_double solZ = ( - ( dot2 * pow ( x2, 2.0 ) * y2 * pow ( z1, 3.0 ) ) /
2007 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2008 ( dot2 * pow ( x2, 2.0 ) * pow ( z1, 2.0 ) ) /
2009 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2010 ( dot1 * pow ( x2, 2.0 ) * y2 * z2 * pow ( z1, 2.0 ) ) /
2011 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2012 ( 2.0 * x1 * dot2 * x2 * y2 * z2 * pow ( z1, 2.0 ) ) /
2013 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2014 ( x2 * y2 * sqrt ( pow ( -2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
2015 4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
2016 ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) * z1 ) /
2017 ( 2.0 * ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2018 ( dot2 * y2 * z1 ) / ( y1 * z2 - z1 * y2 ) +
2019 ( dot1 * pow ( x2, 2.0 ) * z2 * z1 ) /
2020 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2021 ( x1 * dot2 * x2 * z2 * z1 ) /
2022 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2023 ( x1 * dot1 * x2 * pow ( y2, 3.0 ) * z1 ) /
2024 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2025 ( y1 * dot1 * pow ( x2, 2.0 ) * pow ( y2, 2.0 ) * z1 ) /
2026 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2027 ( x1 * y1 * dot2 * x2 * pow ( y2, 2.0 ) * z1 ) /
2028 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2029 ( pow ( x1, 2.0 ) * dot2 * y2 * pow ( z2, 2.0 ) * z1 ) /
2030 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2031 ( 2.0 * x1 * dot1 * x2 * y2 * pow ( z2, 2.0 ) * z1 ) /
2032 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2033 ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) * y2 * z1 ) /
2034 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) + dot2 +
2035 ( x2 * sqrt ( pow ( - 2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
2036 4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
2037 ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) ) /
2038 ( 2.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2039 ( x1 * y2 * z2 * sqrt ( pow ( - 2.0 * dot2 * x2 * pow ( y1, 2.0 ) + 2.0 * x1 * dot2 * y2 * y1 + 2.0 * dot1 * x2 * y2 * y1 - 2.0 * x1 * dot1 * pow ( y2, 2.0 ) - 2.0 * x1 * dot1 * pow ( z2, 2.0 ) - 2.0 * pow ( z1, 2.0 ) * dot2 * x2 + 2.0 * x1 * z1 * dot2 * z2 + 2.0 * z1 * dot1 * x2 * z2, 2.0 ) -
2040 4.0 * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) *
2041 ( pow ( y1, 2.0 ) * pow ( dot2, 2.0 ) + pow ( z1, 2.0 ) * pow ( dot2, 2.0 ) - 2.0 * y1 * dot1 * y2 * dot2 - 2.0 * z1 * dot1 * z2 * dot2 - pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( y2, 2.0 ) - pow ( y1, 2.0 ) * pow ( z2, 2.0 ) + pow ( dot1, 2.0 ) * pow ( z2, 2.0 ) + 2.0 * y1 * z1 * y2 * z2 ) ) ) /
2042 ( 2.0 * ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2043 ( dot1 * y2 * z2 ) / ( y1 * z2 - z1 * y2 ) -
2044 ( pow ( y1, 2.0 ) * dot2 * pow ( x2, 2.0 ) ) /
2045 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2046 ( x1 * dot1 * x2 * pow ( y2, 2.0 ) ) /
2047 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) -
2048 ( x1 * dot1 * x2 * pow ( z2, 2.0 ) ) /
2049 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2050 ( y1 * dot1 * pow ( x2, 2.0 ) * y2 ) /
2051 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2052 ( x1 * y1 * dot2 * x2 * y2 ) /
2053 ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) +
2054 ( pow ( x1, 2.0 ) * dot1 * y2 * pow ( z2, 3.0 ) ) /
2055 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2056 ( pow ( x1, 2.0 ) * dot1 * pow ( y2, 3.0 ) * z2 ) /
2057 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2058 ( pow ( x1, 2.0 ) * y1 * dot2 * pow ( y2, 2.0 ) * z2 ) /
2059 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) -
2060 ( x1 * y1 * dot1 * x2 * pow ( y2, 2.0 ) * z2 ) /
2061 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) +
2062 ( x1 * pow ( y1, 2.0 ) * dot2 * x2 * y2 * z2 ) /
2063 ( ( y1 * z2 - z1 * y2 ) * ( pow ( y1, 2.0 ) * pow ( x2, 2.0 ) + pow ( z1, 2.0 ) * pow ( x2, 2.0 ) - 2.0 * x1 * y1 * y2 * x2 - 2.0 * x1 * z1 * z2 * x2 + pow ( x1, 2.0 ) * pow ( y2, 2.0 ) + pow ( z1, 2.0 ) * pow ( y2, 2.0 ) + pow ( x1, 2.0 ) * pow ( z2, 2.0 ) + pow ( y1, 2.0 ) * pow ( z2, 2.0 ) - 2.0 * y1 * z1 * y2 * z2 ) ) ) / z2;
2066 if ( ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solX ) ) && ( solX < 0.0 ) ) ||
2067 ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solY ) ) && ( solY < 0.0 ) ) ||
2068 ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solZ ) ) && ( solZ < 0.0 ) ) ) { solX *= -1.0; solY *= -1.0; solZ *= -1.0; }
2101 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 )
2104 std::vector < proshade_double > ret;
2107 proshade_double solX = - ( y1 * dot2 * z3 - y1 * dot3 * z2 - z1 * dot2 * y3 + z1 * dot3 * y2 + dot1 * y3 * z2 - dot1 * z3 * y2 ) /
2108 ( -x1 * y3 * z2 + x1 * z3 * y2 + y1 * x3 * z2 - y1 * z3 * x2 - z1 * x3 * y2 + z1 * y3 * x2 );
2109 proshade_double solY = - ( x1 * dot2 * z3 - x1 * dot3 * z2 - z1 * dot2 * x3 + z1 * dot3 * x2 + dot1 * x3 * z2 - dot1 * z3 * x2 ) /
2110 ( x1 * y3 * z2 - x1 * z3 * y2 - y1 * x3 * z2 + y1 * z3 * x2 + z1 * x3 * y2 - z1 * y3 * x2 );
2111 proshade_double solZ = - ( x1 * dot2 * y3 - x1 * dot3 * y2 - y1 * dot2 * x3 + y1 * dot3 * x2 + dot1 * x3 * y2 - dot1 * y3 * x2 ) /
2112 ( -x1 * y3 * z2 + x1 * z3 * y2 + y1 * x3 * z2 - y1 * z3 * x2 - z1 * x3 * y2 + z1 * y3 * x2 );
2115 proshade_double normFactor = sqrt ( pow ( solX, 2.0 ) + pow ( solY, 2.0 ) + pow ( solZ, 2.0 ) );
2121 if ( ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solX ) ) && ( solX < 0.0 ) ) ||
2122 ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solY ) ) && ( solY < 0.0 ) ) ||
2123 ( ( std::max ( std::abs ( solX ), std::max( std::abs ( solY ), std::abs ( solZ ) ) ) == std::abs ( solZ ) ) && ( solZ < 0.0 ) ) ) { solX *= -1.0; solY *= -1.0; solZ *= -1.0; }
2144 std::vector< proshade_double > ret;
2148 ( el1->at(1) * el2->at(3) ) +
2149 ( el1->at(2) * el2->at(6) ) );
2151 ( el1->at(1) * el2->at(4) ) +
2152 ( el1->at(2) * el2->at(7) ) );
2154 ( el1->at(1) * el2->at(5) ) +
2155 ( el1->at(2) * el2->at(8) ) );
2158 ( el1->at(4) * el2->at(3) ) +
2159 ( el1->at(5) * el2->at(6) ) );
2161 ( el1->at(4) * el2->at(4) ) +
2162 ( el1->at(5) * el2->at(7) ) );
2164 ( el1->at(4) * el2->at(5) ) +
2165 ( el1->at(5) * el2->at(8) ) );
2168 ( el1->at(7) * el2->at(3) ) +
2169 ( el1->at(8) * el2->at(6) ) );
2171 ( el1->at(7) * el2->at(4) ) +
2172 ( el1->at(8) * el2->at(7) ) );
2174 ( el1->at(7) * el2->at(5) ) +
2175 ( el1->at(8) * el2->at(8) ) );
2201 proshade_double trace = ( mat1->at(0) * mat2->at(0) ) + ( mat1->at(1) * mat2->at(1) ) + ( mat1->at(2) * mat2->at(2) );
2202 trace += ( mat1->at(3) * mat2->at(3) ) + ( mat1->at(4) * mat2->at(4) ) + ( mat1->at(5) * mat2->at(5) );
2203 trace += ( mat1->at(6) * mat2->at(6) ) + ( mat1->at(7) * mat2->at(7) ) + ( mat1->at(8) * mat2->at(8) );
2209 if ( tolerance > std::abs ( trace ) ) { ret =
true; }
2238 proshade_double cosDist = ( ( a1 * b1 ) + ( a2 * b2 ) + ( a3 * b3 ) ) /
2239 ( sqrt( pow( a1, 2.0 ) + pow( a2, 2.0 ) + pow( a3, 2.0 ) ) *
2240 sqrt( pow( b1, 2.0 ) + pow( b2, 2.0 ) + pow( b3, 2.0 ) ) );
2243 if ( std::abs( cosDist ) > ( 1.0 - tolerance ) ) { ret =
true; }
2272 proshade_double cosDist = ( ( a1 * b1 ) + ( a2 * b2 ) + ( a3 * b3 ) ) /
2273 ( sqrt( pow( a1, 2.0 ) + pow( a2, 2.0 ) + pow( a3, 2.0 ) ) *
2274 sqrt( pow( b1, 2.0 ) + pow( b2, 2.0 ) + pow( b3, 2.0 ) ) );
2277 if ( cosDist > ( 1.0 - tolerance ) ) { ret =
true; }
2296 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 )
2299 proshade_double lonM, lonP, latM, latP, movSum;
2300 std::vector<proshade_double> latVals ( 3 );
2301 std::vector<proshade_double> lonVals ( 3 );
2302 proshade_double learningRate = 0.1;
2303 proshade_double prevVal = *bestSum;
2304 proshade_double valChange = 999.9;
2305 proshade_double origBestLat = std::round ( *bestLattitude );
2306 proshade_double origBestLon = std::round ( *bestLongitude );
2307 proshade_double tmpVal;
2310 std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsMinusMinus;
2311 std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsMinusPlus;
2312 std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsPlusMinus;
2313 std::vector<ProSHADE_internal_maths::BicubicInterpolator*> interpolsPlusPlus;
2320 while ( valChange > 0.0001 )
2323 lonM = *bestLongitude - step;
2324 lonP = *bestLongitude + step;
2325 latM = *bestLattitude - step;
2326 latP = *bestLattitude + step;
2329 if ( latM < ( origBestLat - 1.0 ) ) { tmpVal = *bestLattitude; *bestLattitude = origBestLat - 1.0;
optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step );
if ( *bestLattitude == origBestLat - 1.0 ) { *bestLattitude = tmpVal; }
break; }
2330 if ( latP > ( origBestLat + 1.0 ) ) { tmpVal = *bestLattitude; *bestLattitude = origBestLat + 1.0;
optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step );
if ( *bestLattitude == origBestLat + 1.0 ) { *bestLattitude = tmpVal; }
break; }
2331 if ( lonM < ( origBestLon - 1.0 ) ) { tmpVal = *bestLongitude; *bestLongitude = origBestLon - 1.0;
optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step );
if ( *bestLongitude == origBestLon - 1.0 ) { *bestLongitude = tmpVal; }
break; }
2332 if ( lonP > ( origBestLon + 1.0 ) ) { tmpVal = *bestLongitude; *bestLongitude = origBestLon + 1.0;
optimiseAxisBiCubicInterpolation ( bestLattitude, bestLongitude, bestSum, sphereList, sphereMappedRotFun, step );
if ( *bestLongitude == origBestLon + 1.0 ) { *bestLongitude = tmpVal; }
break; }
2335 latVals.at(0) = latM; latVals.at(1) = *bestLattitude; latVals.at(2) = latP;
2336 lonVals.at(0) = lonM; lonVals.at(1) = *bestLongitude; lonVals.at(2) = lonP;
2339 for ( proshade_unsign laIt = 0; laIt < static_cast<proshade_unsign> ( latVals.size() ); laIt++ )
2341 for ( proshade_unsign loIt = 0; loIt < static_cast<proshade_unsign> ( lonVals.size() ); loIt++ )
2345 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sphereList->size() ); iter++ )
2348 if ( ( latVals.at(laIt) <= origBestLat ) && ( lonVals.at(loIt) <= origBestLon ) ) { movSum += interpolsMinusMinus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2349 if ( ( latVals.at(laIt) <= origBestLat ) && ( lonVals.at(loIt) > origBestLon ) ) { movSum += interpolsMinusPlus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2350 if ( ( latVals.at(laIt) > origBestLat ) && ( lonVals.at(loIt) <= origBestLon ) ) { movSum += interpolsPlusMinus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2351 if ( ( latVals.at(laIt) > origBestLat ) && ( lonVals.at(loIt) > origBestLon ) ) { movSum += interpolsPlusPlus.at(iter)->getValue ( latVals.at(laIt), lonVals.at(loIt) ); }
2355 if ( *bestSum < movSum )
2358 *bestLongitude = lonVals.at(loIt);
2359 *bestLattitude = latVals.at(laIt);
2365 valChange = std::floor ( 100000.0 * ( *bestSum - prevVal ) ) / 100000.0;
2366 prevVal = std::floor ( 100000.0 * ( *bestSum ) ) / 100000.0;
2367 step = std::max ( ( valChange / step ) * learningRate, 0.01 );
2368 if ( learningRate >= 0.02 ) { learningRate -= 0.01; }
2373 for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsMinusMinus.size() ); intIt++ ) {
delete interpolsMinusMinus.at(intIt); }
2374 for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsMinusPlus.size() ); intIt++ ) {
delete interpolsMinusPlus.at(intIt); }
2375 for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsPlusMinus.size() ); intIt++ ) {
delete interpolsPlusMinus.at(intIt); }
2376 for ( proshade_unsign intIt = 0; intIt < static_cast<proshade_unsign> ( interpolsPlusPlus.size() ); intIt++ ) {
delete interpolsPlusPlus.at(intIt); }
2393 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 )
2396 proshade_signed latHlp, lonHlp;
2397 proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2400 for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2403 proshade_double** interpGrid =
new proshade_double*[4];
2407 for ( proshade_unsign iter = 0; iter < 4; iter++ )
2409 interpGrid[iter] =
new proshade_double[4];
2414 for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2416 for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2418 latHlp = bestLattitude - 2 + latIt;
if ( latHlp < 0.0 ) { latHlp += angDim; }
if ( latHlp >= angDim ) { latHlp -= angDim; }
2419 lonHlp = bestLongitude - 2 + lonIt;
if ( lonHlp < 0.0 ) { lonHlp += angDim; }
if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2420 interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2426 interpols->emplace_back ( biCubInterp );
2429 for ( proshade_unsign iter = 0; iter < 4; iter++ ) {
delete[] interpGrid[iter]; }
2430 delete[] interpGrid;
2448 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 )
2451 proshade_signed latHlp, lonHlp;
2452 proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2455 for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2458 proshade_double** interpGrid =
new proshade_double*[4];
2462 for ( proshade_unsign iter = 0; iter < 4; iter++ )
2464 interpGrid[iter] =
new proshade_double[4];
2469 for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2471 for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2473 latHlp = bestLattitude - 2 + latIt;
if ( latHlp < 0.0 ) { latHlp += angDim; }
if ( latHlp >= angDim ) { latHlp -= angDim; }
2474 lonHlp = bestLongitude - 1 + lonIt;
if ( lonHlp < 0.0 ) { lonHlp += angDim; }
if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2475 interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2481 interpols->emplace_back ( biCubInterp );
2484 for ( proshade_unsign iter = 0; iter < 4; iter++ ) {
delete[] interpGrid[iter]; }
2485 delete[] interpGrid;
2503 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 )
2506 proshade_signed latHlp, lonHlp;
2507 proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2510 for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2513 proshade_double** interpGrid =
new proshade_double*[4];
2517 for ( proshade_unsign iter = 0; iter < 4; iter++ )
2519 interpGrid[iter] =
new proshade_double[4];
2524 for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2526 for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2528 latHlp = bestLattitude - 1 + latIt;
if ( latHlp < 0.0 ) { latHlp += angDim; }
if ( latHlp >= angDim ) { latHlp -= angDim; }
2529 lonHlp = bestLongitude - 2 + lonIt;
if ( lonHlp < 0.0 ) { lonHlp += angDim; }
if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2530 interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2536 interpols->emplace_back ( biCubInterp );
2539 for ( proshade_unsign iter = 0; iter < 4; iter++ ) {
delete[] interpGrid[iter]; }
2540 delete[] interpGrid;
2558 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 )
2561 proshade_signed latHlp, lonHlp;
2562 proshade_signed angDim = sphereMappedRotFun->at(0)->getAngularDim();
2565 for ( proshade_unsign sphereIt = 0; sphereIt < static_cast<proshade_unsign> ( sphereList->size() ); sphereIt++ )
2568 proshade_double** interpGrid =
new proshade_double*[4];
2572 for ( proshade_unsign iter = 0; iter < 4; iter++ )
2574 interpGrid[iter] =
new proshade_double[4];
2579 for ( proshade_unsign latIt = 0; latIt < 4; latIt++ )
2581 for ( proshade_unsign lonIt = 0; lonIt < 4; lonIt++ )
2583 latHlp = bestLattitude - 1 + latIt;
if ( latHlp < 0.0 ) { latHlp += angDim; }
if ( latHlp >= angDim ) { latHlp -= angDim; }
2584 lonHlp = bestLongitude - 1 + lonIt;
if ( lonHlp < 0.0 ) { lonHlp += angDim; }
if ( lonHlp >= angDim ) { lonHlp -= angDim; }
2585 interpGrid[latIt][lonIt] = sphereMappedRotFun->at(sphereList->at(sphereIt))->getSphereLatLonPosition ( latHlp, lonHlp );
2591 interpols->emplace_back ( biCubInterp );
2594 for ( proshade_unsign iter = 0; iter < 4; iter++ ) {
delete[] interpGrid[iter]; }
2595 delete[] interpGrid;
2617 proshade_unsign whichImprove;
2620 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( CSymList->size() ); grIt++ )
2623 if ( CSymList->at(grIt)[0] == axis[0] )
2628 whichImprove = grIt;
2635 if ( improve && !ret )
2637 CSymList->at(whichImprove)[1] = axis[1];
2638 CSymList->at(whichImprove)[2] = axis[2];
2639 CSymList->at(whichImprove)[3] = axis[3];
2640 CSymList->at(whichImprove)[4] = axis[4];
2641 CSymList->at(whichImprove)[5] = axis[5];
2667 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( CSymList->size() ); grIt++ )
2669 if ( fold == CSymList->at(grIt)[0] )
2678 std::cout << std::endl;
2696 std::vector< proshade_unsign > ret;
2697 std::vector< std::pair< proshade_unsign, bool > > sieveOfEratosthenesArray;
2700 if ( upTo < 2 ) {
return ( ret ); }
2703 for ( proshade_unsign iter = 2; iter <= upTo; iter++ ) { sieveOfEratosthenesArray.emplace_back ( std::pair< proshade_unsign, bool > ( iter,
true ) ); }
2706 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); iter++ )
2709 if ( sieveOfEratosthenesArray.at(iter).second )
2712 for ( proshade_unsign it = iter + sieveOfEratosthenesArray.at(iter).first; it < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); it += sieveOfEratosthenesArray.at(iter).first )
2714 sieveOfEratosthenesArray.at(it).second =
false;
2720 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( sieveOfEratosthenesArray.size() ); iter++ )