71 proshade_double xSamplRate = ( this->xDimSizeOriginal /
static_cast<proshade_double
> ( this->xDimIndicesOriginal ) );
72 proshade_double ySamplRate = ( this->yDimSizeOriginal /
static_cast<proshade_double
> ( this->yDimIndicesOriginal ) );
73 proshade_double zSamplRate = ( this->zDimSizeOriginal /
static_cast<proshade_double
> ( this->zDimIndicesOriginal ) );
76 proshade_double xRotPos = (
static_cast<proshade_double
> ( this->xFrom - this->mapMovFromsChangeX ) * xSamplRate ) +
77 ( ( (
static_cast<proshade_double
> ( this->xTo ) -
static_cast<proshade_double
> ( this->xFrom ) ) / 2.0 ) * xSamplRate );
79 proshade_double yRotPos = (
static_cast<proshade_double
> ( this->yFrom - this->mapMovFromsChangeY ) * ySamplRate ) +
80 ( ( (
static_cast<proshade_double
> ( this->yTo ) -
static_cast<proshade_double
> ( this->yFrom ) ) / 2.0 ) * ySamplRate );
82 proshade_double zRotPos = (
static_cast<proshade_double
> ( this->zFrom - this->mapMovFromsChangeZ ) * zSamplRate ) +
83 ( ( (
static_cast<proshade_double
> ( this->zTo ) -
static_cast<proshade_double
> ( this->zFrom ) ) / 2.0 ) * zSamplRate );
86 this->originalPdbRotCenX = xRotPos - ( this->mapCOMProcessChangeX / 2.0 );
87 this->originalPdbRotCenY = yRotPos - ( this->mapCOMProcessChangeY / 2.0 );
88 this->originalPdbRotCenZ = zRotPos - ( this->mapCOMProcessChangeZ / 2.0 );
112 this->originalPdbTransX = 0.0;
113 this->originalPdbTransY = 0.0;
114 this->originalPdbTransZ = 0.0;
117 if ( ( eulA != 0.0 ) || ( eulB != 0.0 ) || ( eulG != 0.0 ) )
125 this->originalPdbTransX = this->mapCOMProcessChangeX;
126 this->originalPdbTransY = this->mapCOMProcessChangeY;
127 this->originalPdbTransZ = this->mapCOMProcessChangeZ;
131 this->originalPdbTransX += trsX;
132 this->originalPdbTransY += trsY;
133 this->originalPdbTransZ += trsZ;
178 eulA, eulB, eulG, settings );
218 movingStructure->
rotateMap ( settings, eulA, eulB, eulG );
222 std::max ( staticStructure->
getYDim(), movingStructure->
getYDim() ),
223 std::max ( staticStructure->
getZDim(), movingStructure->
getZDim() ) );
225 std::max ( staticStructure->
getYDim(), movingStructure->
getYDim() ),
226 std::max ( staticStructure->
getZDim(), movingStructure->
getZDim() ) );
235 proshade_double mapPeak = 0.0;
236 proshade_unsign xDimS = staticStructure->
getXDim();
237 proshade_unsign yDimS = staticStructure->
getYDim();
238 proshade_unsign zDimS = staticStructure->
getZDim();
242 if ( *trsX > ( xDimS / 2 ) ) { *trsX = *trsX - xDimS; }
243 if ( *trsY > ( yDimS / 2 ) ) { *trsY = *trsY - yDimS; }
244 if ( *trsZ > ( zDimS / 2 ) ) { *trsZ = *trsZ - zDimS; }
247 proshade_single xCor = ( staticStructure->
xFrom - movingStructure->
xFrom ) *
248 (
static_cast<proshade_double
> ( staticStructure->
getXDimSize() ) / staticStructure->
getXDim() );
249 proshade_single yCor = ( staticStructure->
yFrom - movingStructure->
yFrom ) *
250 (
static_cast<proshade_double
> ( staticStructure->
getYDimSize() ) / staticStructure->
getYDim() );
251 proshade_single zCor = ( staticStructure->
zFrom - movingStructure->
zFrom ) *
252 (
static_cast<proshade_double
> ( staticStructure->
getZDimSize() ) / staticStructure->
getZDim() );
254 ( *trsX *
static_cast<proshade_double
> ( staticStructure->
getXDimSize() ) / staticStructure->
getXDim() );
256 ( *trsY *
static_cast<proshade_double
> ( staticStructure->
getYDimSize() ) / staticStructure->
getYDim() );
258 ( *trsZ *
static_cast<proshade_double
> ( staticStructure->
getZDimSize() ) / staticStructure->
getZDim() );
266 std::stringstream hlpSS;
267 hlpSS <<
"Optimal map translation distances are " << *trsX <<
" ; " << *trsY <<
" ; " << *trsZ <<
" Angstroms with peak height " << mapPeak / ( xDimS * yDimS * zDimS );
311 std::vector< proshade_double > ret;
312 proshade_double mapPeak = 0.0;
313 proshade_double trsX = 0.0, trsY = 0.0, trsZ = 0.0;
324 if ( trsX > ( staticStructure->
getXDim() / 2 ) ) { trsX = trsX - this->getXDim(); }
325 if ( trsY > ( staticStructure->
getYDim() / 2 ) ) { trsY = trsY - this->getYDim(); }
326 if ( trsZ > ( staticStructure->
getZDim() / 2 ) ) { trsZ = trsZ - this->getZDim(); }
329 proshade_single xCor = ( staticStructure->
xFrom - this->xFrom ) *
330 (
static_cast<proshade_double
> ( staticStructure->
getXDimSize() ) / staticStructure->
getXDim() );
331 proshade_single yCor = ( staticStructure->
yFrom - this->yFrom ) *
332 (
static_cast<proshade_double
> ( staticStructure->
getYDimSize() ) / staticStructure->
getYDim() );
333 proshade_single zCor = ( staticStructure->
zFrom - this->zFrom ) *
334 (
static_cast<proshade_double
> ( staticStructure->
getZDimSize() ) / staticStructure->
getZDim() );
336 ( trsX *
static_cast<proshade_double
> ( staticStructure->
getXDimSize() ) / staticStructure->
getXDim() );
338 ( trsY *
static_cast<proshade_double
> ( staticStructure->
getYDimSize() ) / staticStructure->
getYDim() );
340 ( trsZ *
static_cast<proshade_double
> ( staticStructure->
getZDimSize() ) / staticStructure->
getZDim() );
341 proshade_single modXMov = xMov;
342 proshade_single modYMov = yMov;
343 proshade_single modZMov = zMov;
351 this->mapMovFromsChangeX = this->xFrom;
352 this->mapMovFromsChangeY = this->yFrom;
353 this->mapMovFromsChangeZ = this->zFrom;
357 this->getXFromPtr(), this->getXToPtr(),
358 this->getYFromPtr(), this->getYToPtr(),
359 this->getZFromPtr(), this->getZToPtr(),
360 this->getXAxisOrigin(), this->getYAxisOrigin(), this->getZAxisOrigin() );
363 this->getXDimSize(), this->getYDimSize(), this->getZDimSize(),
364 this->getXDim(), this->getYDim(), this->getZDim() );
367 this->mapMovFromsChangeX = this->xFrom - this->mapMovFromsChangeX;
368 this->mapMovFromsChangeY = this->yFrom - this->mapMovFromsChangeY;
369 this->mapMovFromsChangeZ = this->zFrom - this->mapMovFromsChangeZ;
372 this->computePdbRotationCentre ( );
373 this->computeOptimalTranslation ( eulA, eulB, eulG, -xMov, -yMov, -zMov );
392 fftw_complex *tmpIn1 = NULL, *tmpOut1 = NULL, *tmpIn2 = NULL, *tmpOut2 = NULL, *resOut = NULL;
393 fftw_plan forwardFourierObj1, forwardFourierObj2, inverseFourierCombo;
394 proshade_unsign dimMult = staticStructure->
getXDim() * staticStructure->
getYDim() * staticStructure->
getZDim();
395 ProSHADE_internal_overlay::allocateTranslationFunctionMemory ( tmpIn1, tmpOut1, tmpIn2, tmpOut2, this->translationMap, resOut, forwardFourierObj1, forwardFourierObj2, inverseFourierCombo, staticStructure->
getXDim(), staticStructure->
getYDim(), staticStructure->
getZDim() );
398 for ( proshade_unsign iter = 0; iter < dimMult; iter++ ) { tmpIn1[iter][0] = staticStructure->
getMapValue ( iter ); tmpIn1[iter][1] = 0.0; }
399 for ( proshade_unsign iter = 0; iter < dimMult; iter++ ) { tmpIn2[iter][0] = this->getMapValue ( iter ); tmpIn2[iter][1] = 0.0; }
402 fftw_execute ( forwardFourierObj1 );
403 fftw_execute ( forwardFourierObj2 );
407 fftw_execute ( inverseFourierCombo );
432 void ProSHADE_internal_overlay::allocateTranslationFunctionMemory ( fftw_complex*& tmpIn1, fftw_complex*& tmpOut1, fftw_complex*& tmpIn2, fftw_complex*& tmpOut2, fftw_complex*& resIn, fftw_complex*& resOut, fftw_plan& forwardFourierObj1, fftw_plan& forwardFourierObj2, fftw_plan& inverseFourierCombo, proshade_unsign xD, proshade_unsign yD, proshade_unsign zD )
435 tmpIn1 =
new fftw_complex[xD * yD * zD];
436 tmpOut1 =
new fftw_complex[xD * yD * zD];
437 tmpIn2 =
new fftw_complex[xD * yD * zD];
438 tmpOut2 =
new fftw_complex[xD * yD * zD];
439 resIn =
new fftw_complex[xD * yD * zD];
440 resOut =
new fftw_complex[xD * yD * zD];
451 forwardFourierObj1 = fftw_plan_dft_3d ( xD, yD, zD, tmpIn1, tmpOut1, FFTW_FORWARD , FFTW_ESTIMATE );
452 forwardFourierObj2 = fftw_plan_dft_3d ( xD, yD, zD, tmpIn2, tmpOut2, FFTW_FORWARD , FFTW_ESTIMATE );
453 inverseFourierCombo = fftw_plan_dft_3d ( xD, yD, zD, resOut, resIn , FFTW_BACKWARD, FFTW_ESTIMATE );
474 fftw_destroy_plan ( forwardFourierObj1 );
475 fftw_destroy_plan ( forwardFourierObj2 );
476 fftw_destroy_plan ( inverseFourierCombo );
500 double normFactor =
static_cast<double> ( xD * yD * zD );
501 proshade_signed h1, k1, l1, hlpPos, arrPos;
502 proshade_signed uo2 = xD / 2;
503 proshade_signed vo2 = yD / 2;
504 proshade_signed wo2 = zD / 2;
507 for ( proshade_signed uIt = 0; uIt < static_cast<proshade_signed> ( xD ); uIt++ )
509 for ( proshade_signed vIt = 0; vIt < static_cast<proshade_signed> ( yD ); vIt++ )
511 for ( proshade_signed wIt = 0; wIt < static_cast<proshade_signed> ( zD ); wIt++ )
514 if ( uIt > uo2 ) { h1 = uIt - xD; }
else { h1 = uIt; }
515 if ( vIt > vo2 ) { k1 = vIt - yD; }
else { k1 = vIt; }
516 if ( wIt > wo2 ) { l1 = wIt - zD; }
else { l1 = wIt; }
519 if ( h1 < 0 ) { h1 += xD; }
520 if ( k1 < 0 ) { k1 += yD; }
521 if ( l1 < 0 ) { l1 += zD; }
524 hlpPos = l1 + zD * ( k1 + yD * h1 );
525 arrPos = wIt + zD * ( vIt + yD * uIt );
533 &resOut[hlpPos][1] );
536 resOut[hlpPos][0] /= normFactor;
537 resOut[hlpPos][1] /= normFactor;
561 proshade_signed arrPos;
565 for ( proshade_signed uIt = 0; uIt < static_cast<proshade_signed> ( xD ); uIt++ )
567 for ( proshade_signed vIt = 0; vIt < static_cast<proshade_signed> ( yD ); vIt++ )
569 for ( proshade_signed wIt = 0; wIt < static_cast<proshade_signed> ( zD ); wIt++ )
571 arrPos = wIt + zD * ( vIt + yD * uIt );
572 if ( resIn[arrPos][0] > *mapPeak )
574 *mapPeak = resIn[arrPos][0];
601 if ( ( this->xDimIndices > xDim ) || ( this->yDimIndices > yDim ) || ( this->zDimIndices > zDim ) )
603 throw ProSHADE_exception (
"Cannot zero-pad in negative direction.",
"EO00034", __FILE__, __LINE__, __func__,
"The requested padded size of a structure is smaller than\n : the current size. If the user sees this error, there is\n : likely a considerable bug. Please report this error." );
607 if ( ( this->xDimIndices == xDim ) && ( this->yDimIndices == yDim ) && ( this->zDimIndices == zDim ) ) { return ; }
610 proshade_unsign addXPre, addYPre, addZPre, addXPost, addYPost, addZPost;
611 ProSHADE_internal_overlay::computeBeforeAfterZeroCounts ( &addXPre, &addYPre, &addZPre, &addXPost, &addYPost, &addZPost, xDim, yDim, zDim, this->xDimIndices, this->yDimIndices, this->zDimIndices );
614 proshade_double* newMap =
new proshade_double [xDim * yDim * zDim];
617 ProSHADE_internal_overlay::paddMapWithZeroes ( this->internalMap, newMap, xDim, yDim, zDim, this->xDimIndices, this->yDimIndices, this->zDimIndices, addXPre, addYPre, addZPre );
620 delete[] this->internalMap;
621 this->internalMap =
new proshade_double [xDim * yDim * zDim];
622 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { this->internalMap[iter] = newMap[iter]; }
628 this->xDimSize = xDim * ( this->xDimSize / this->xDimIndices );
629 this->yDimSize = yDim * ( this->yDimSize / this->yDimIndices );
630 this->zDimSize = zDim * ( this->zDimSize / this->zDimIndices );
631 this->xDimIndices = xDim ; this->yDimIndices = yDim ; this->zDimIndices = zDim;
632 this->xGridIndices = xDim ; this->yGridIndices = yDim ; this->zGridIndices = zDim;
633 this->xFrom -= addXPre ; this->yFrom -= addYPre ; this->zFrom -= addZPre;
634 this->xTo += addXPost; this->yTo += addYPost; this->zTo += addZPost;
635 this->xAxisOrigin -= addXPre ; this->yAxisOrigin -= addYPre ; this->zAxisOrigin -= addZPre ;
657 void ProSHADE_internal_overlay::computeBeforeAfterZeroCounts ( proshade_unsign* addXPre, proshade_unsign* addYPre, proshade_unsign* addZPre, proshade_unsign* addXPost, proshade_unsign* addYPost, proshade_unsign* addZPost, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_unsign xDimIndices, proshade_unsign yDimIndices, proshade_unsign zDimIndices )
660 *addXPre = ( xDim - xDimIndices ) / 2;
661 *addYPre = ( yDim - yDimIndices ) / 2;
662 *addZPre = ( zDim - zDimIndices ) / 2;
663 *addXPost = *addXPre;
if ( ( xDim - xDimIndices ) % 2 == 1 ) { *addXPost += 1; }
664 *addYPost = *addYPre;
if ( ( yDim - yDimIndices ) % 2 == 1 ) { *addYPost += 1; }
665 *addZPost = *addZPre;
if ( ( zDim - zDimIndices ) % 2 == 1 ) { *addZPost += 1; }
686 void ProSHADE_internal_overlay::paddMapWithZeroes ( proshade_double* oldMap, proshade_double*& newMap, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_unsign xDimIndices, proshade_unsign yDimIndices, proshade_unsign zDimIndices, proshade_unsign addXPre, proshade_unsign addYPre, proshade_unsign addZPre )
689 proshade_unsign newMapIndex = 0;
690 proshade_unsign oldMapIndex = 0;
693 for ( proshade_unsign xIt = 0; xIt < xDim; xIt++ )
695 for ( proshade_unsign yIt = 0; yIt < yDim; yIt++ )
697 for ( proshade_unsign zIt = 0; zIt < zDim; zIt++ )
700 newMapIndex = zIt + zDim * ( yIt + yDim * xIt );
703 if ( xIt < addXPre ) { newMap[newMapIndex] = 0.0;
continue; }
704 if ( yIt < addYPre ) { newMap[newMapIndex] = 0.0;
continue; }
705 if ( zIt < addZPre ) { newMap[newMapIndex] = 0.0;
continue; }
708 if ( xIt >= ( addXPre + xDimIndices ) ) { newMap[newMapIndex] = 0.0;
continue; }
709 if ( yIt >= ( addYPre + yDimIndices ) ) { newMap[newMapIndex] = 0.0;
continue; }
710 if ( zIt >= ( addZPre + zDimIndices ) ) { newMap[newMapIndex] = 0.0;
continue; }
713 oldMapIndex = (zIt-addZPre) + zDimIndices * ( (yIt-addYPre) + yDimIndices * (xIt-addXPre) );
714 newMap[newMapIndex] = oldMap[oldMapIndex];
739 this->maxCompBand = this->spheres[this->noSpheres-1]->getLocalBandwidth();
742 this->findMapCOM ( );
743 this->mapCOMProcessChangeX = this->xCom - this->originalMapXCom;
744 this->mapCOMProcessChangeY = this->yCom - this->originalMapYCom;
745 this->mapCOMProcessChangeZ = this->zCom - this->originalMapZCom;
751 this->allocateRotatedSHMemory ( settings );
754 this->computeRotatedSH ( settings );
757 this->invertSHCoefficients ( );
760 std::vector<proshade_double> lonCO, latCO;
764 proshade_double *densityMapRotated =
new proshade_double [this->xDimIndices * this->yDimIndices * this->zDimIndices];
766 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { densityMapRotated[iter] = 0.0; }
769 this->interpolateMapFromSpheres ( settings, densityMapRotated );
772 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
774 this->internalMap[iter] = densityMapRotated[iter];
778 delete[] densityMapRotated;
798 proshade_single xMov = -trsX, yMov = -trsY, zMov = -trsZ;
802 this->getXFromPtr(), this->getXToPtr(), this->getYFromPtr(), this->getYToPtr(),
803 this->getZFromPtr(), this->getZToPtr(), this->getXAxisOrigin(), this->getYAxisOrigin(), this->getZAxisOrigin() );
807 this->getXDim(), this->getYDim(), this->getZDim() );
821 this->rotSphericalHarmonics =
new proshade_complex* [this->noSpheres];
825 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
828 this->rotSphericalHarmonics[iter] =
new proshade_complex [
static_cast<proshade_unsign
> ( pow ( (this->spheres[iter]->getLocalBandwidth() * 2), 2) )];
832 for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( pow ( (this->spheres[iter]->getLocalBandwidth() * 2), 2) ); it++ )
834 this->rotSphericalHarmonics[iter][it][0] = 0.0;
835 this->rotSphericalHarmonics[iter][it][1] = 0.0;
851 proshade_double WigDR, WigDI, *ShR, *ShI, retR, retI;
852 proshade_unsign arrPos;
855 for ( proshade_signed shell = 0; shell < static_cast<proshade_signed> ( this->noSpheres ); shell++ )
858 for ( proshade_signed bandIter = 0; bandIter < static_cast<proshade_signed> ( this->spheres[shell]->getLocalBandwidth() ); bandIter++ )
861 for ( proshade_signed order1 = 0; order1 < ( ( bandIter * 2 ) + 1 ); order1++ )
864 ShR = getRealSphHarmValue ( bandIter, order1, shell );
865 ShI = getImagSphHarmValue ( bandIter, order1, shell );
868 for ( proshade_signed order2 = 0; order2 < ( ( bandIter * 2 ) + 1 ); order2++ )
871 this->getWignerMatrixValue ( bandIter, order1, order2, &WigDR, &WigDI );
877 arrPos =
static_cast<proshade_unsign
> ( seanindex ( order2-bandIter, bandIter, this->spheres[shell]->getLocalBandwidth() ) );
878 this->rotSphericalHarmonics[shell][arrPos][0] += retR;
879 this->rotSphericalHarmonics[shell][arrPos][1] += retI;
905 proshade_unsign oneDim = shBand * 2;
908 sigR =
new double [(oneDim * oneDim * 4)];
909 sigI = sigR + (oneDim * oneDim * 2);
910 rcoeffs =
new double [(oneDim * oneDim * 2)];
911 icoeffs = rcoeffs + (oneDim * oneDim);
912 weights =
new double [4 * shBand];
913 workspace =
new double [( 20 * shBand * shBand ) + ( 24 * shBand )];
922 idctPlan = fftw_plan_r2r_1d ( oneDim, weights, workspace, FFTW_REDFT01, FFTW_ESTIMATE );
925 int rank, howmany_rank;
926 fftw_iodim dims[1], howmany_dims[1];
929 dims[0].n = 2 * shBand;
930 dims[0].is = 2 * shBand;
933 howmany_dims[0].n = 2 * shBand;
934 howmany_dims[0].is = 1;
935 howmany_dims[0].os = 2 * shBand;
938 ifftPlan = fftw_plan_guru_split_dft ( rank, dims, howmany_rank, howmany_dims, sigR, sigI, rcoeffs, icoeffs, FFTW_ESTIMATE );
951 double *sigR = NULL, *sigI = NULL, *rcoeffs = NULL, *icoeffs = NULL, *weights = NULL, *workspace = NULL;
952 fftw_plan idctPlan, ifftPlan;
955 for (
int shell = 0; shell < static_cast<int> ( this->noSpheres ); shell++ )
958 proshade_unsign oneDim = this->spheres[shell]->getLocalBandwidth() * 2;
964 makeweights ( this->spheres[shell]->getLocalBandwidth(), weights );
967 this->spheres[shell]->allocateRotatedMap ( );
970 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( oneDim * oneDim ); iter++ )
972 rcoeffs[iter] = this->rotSphericalHarmonics[shell][iter][0];
973 icoeffs[iter] = this->rotSphericalHarmonics[shell][iter][1];
979 InvFST_semi_fly ( rcoeffs,
983 this->spheres[shell]->getLocalBandwidth(),
986 this->spheres[shell]->getLocalBandwidth(),
991 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( oneDim * oneDim ); iter++ )
993 this->spheres[shell]->setRotatedMappedData ( iter, sigR[iter] );
997 fftw_destroy_plan ( idctPlan );
998 fftw_destroy_plan ( ifftPlan );
1021 for ( proshade_unsign iter = 0; iter <= angRes; iter++ )
1023 ProSHADE_internal_misc::addToDoubleVector ( lonCO,
static_cast<proshade_double
> ( iter ) * ( (
static_cast<proshade_double
> ( M_PI ) * 2.0 ) /
static_cast<proshade_double
> ( angRes ) ) - (
static_cast<double> ( M_PI ) ) );
1027 for ( proshade_unsign iter = 0; iter <= angRes; iter++ )
1029 ProSHADE_internal_misc::addToDoubleVector ( latCO,
static_cast<proshade_double
> ( iter ) * (
static_cast<proshade_double
> ( M_PI ) /
static_cast<proshade_double
> ( angRes ) ) - (
static_cast<proshade_double
> ( M_PI ) / 2.0 ) );
1045 proshade_double rad = 0.0, lon = 0.0, lat = 0.0, newU = 0.0, newV = 0.0, newW = 0.0;
1046 proshade_unsign lowerLonL = 0, upperLonL = 0, lowerLonU = 0, upperLonU = 0, lowerLatL = 0, upperLatL = 0, lowerLatU = 0, upperLatU = 0, lowerShell = 0, upperShell = 0;
1047 proshade_double x00 = 0.0, x01 = 0.0, x10 = 0.0, x11 = 0.0, distLLon = 0.0, distLLat = 0.0, distLRad = 0.0, valLLon = 0.0, valULon = 0.0;
1048 proshade_double lowerShellValue = 0.0, upperShellValue = 0.0;
1049 proshade_double xSamplingRate =
static_cast<proshade_double
> ( this->xDimSize ) / this->xDimIndices;
1050 proshade_double ySamplingRate =
static_cast<proshade_double
> ( this->yDimSize ) / this->yDimIndices;
1051 proshade_double zSamplingRate =
static_cast<proshade_double
> ( this->zDimSize ) / this->zDimIndices;
1052 proshade_signed arrPos;
1053 std::vector<proshade_double> lonCOU, latCOU, lonCOL, latCOL;
1055 for ( proshade_signed uIt = 0; uIt < static_cast<proshade_signed> (this->xDimIndices); uIt++ )
1057 for ( proshade_signed vIt = 0; vIt < static_cast<proshade_signed> (this->yDimIndices); vIt++ )
1059 for ( proshade_signed wIt = 0; wIt < static_cast<proshade_signed> (this->zDimIndices); wIt++ )
1062 newU =
static_cast<proshade_double
> ( uIt - (
static_cast<proshade_signed
> (this->xDimIndices) / 2 ) );
1063 newV =
static_cast<proshade_double
> ( vIt - (
static_cast<proshade_signed
> (this->yDimIndices) / 2 ) );
1064 newW =
static_cast<proshade_double
> ( wIt - (
static_cast<proshade_signed
> (this->zDimIndices) / 2 ) );
1067 if ( ( newU == 0.0 ) && ( newV == 0.0 ) && ( newW == 0.0 ) )
1069 arrPos = wIt + this->zDimIndices * ( vIt + this->yDimIndices * uIt );
1070 densityMapRotated[arrPos] = this->internalMap[arrPos];
1075 rad = sqrt ( pow( ( newU * xSamplingRate ), 2.0 ) +
1076 pow( ( newV * ySamplingRate ), 2.0 ) +
1077 pow( ( newW * zSamplingRate ), 2.0 ) );
1078 lon = atan2 ( ( newV * ySamplingRate ), ( newU * xSamplingRate ) );
1079 lat = asin ( ( newW * zSamplingRate ) / rad );
1082 if ( rad != rad ) { rad = 0.0; }
1083 if ( lon != lon ) { lon = 0.0; }
1084 if ( lat != lat ) { lat = 0.0; }
1089 for ( proshade_unsign iter = 0; iter < (this->noSpheres-1); iter++ )
1091 if ( ( this->spherePos.at(iter) <= rad ) && ( this->spherePos.at(iter+1) > rad ) )
1094 upperShell = iter+1;
1099 if ( upperShell == 0 )
1101 arrPos = wIt + this->zDimIndices * ( vIt + this->yDimIndices * uIt );
1102 densityMapRotated[arrPos] = 0.0;
1107 lonCOL.clear(); latCOL.clear(); lonCOU.clear(); latCOU.clear();
1112 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( lonCOL.size() ); iter++ )
1114 if ( iter == (
static_cast<proshade_unsign
> ( lonCOL.size() ) - 1 ) )
1120 if ( ( std::floor(10000. * lonCOL.at(iter)) <= std::floor(10000. * lon) ) && ( std::floor(10000. * lonCOL.at(iter+1)) > std::floor(10000. * lon) ) )
1127 if ( upperLonL == this->spheres[lowerShell]->getLocalAngRes() ) { upperLonL = 0; }
1129 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( lonCOU.size() ); iter++ )
1131 if ( iter == (
static_cast<proshade_unsign
> ( lonCOU.size() ) - 1 ) )
1137 if ( ( std::floor(10000. * lonCOU.at(iter)) <= std::floor(10000. * lon) ) && ( std::floor(10000. * lonCOU.at(iter+1)) > std::floor(10000. * lon) ) )
1144 if ( upperLonU == this->spheres[upperShell]->getLocalAngRes() ) { upperLonU = 0; }
1146 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( latCOL.size() ); iter++ )
1148 if ( iter == (
static_cast<proshade_unsign
> ( latCOL.size() ) - 1 ) )
1154 if ( ( std::floor(10000. * latCOL.at(iter)) <= std::floor(10000. * lat) ) && ( std::floor(10000. * latCOL.at(iter+1)) > std::floor(10000. * lat) ) )
1161 if ( upperLatL == this->spheres[lowerShell]->getLocalAngRes() ) { upperLatL = 0; }
1163 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( latCOU.size() ); iter++ )
1165 if ( iter == (
static_cast<proshade_unsign
> ( latCOU.size() ) - 1 ) )
1171 if ( ( std::floor(10000. * latCOU.at(iter)) <= std::floor(10000. * lat) ) && ( std::floor(10000. * latCOU.at(iter+1)) > std::floor(10000. * lat) ) )
1178 if ( upperLatU == this->spheres[upperShell]->getLocalAngRes() ) { upperLatU = 0; }
1181 x00 = this->spheres[lowerShell]->getRotatedMappedData ( lowerLatL * this->spheres[lowerShell]->getLocalAngRes() + lowerLonL );
1182 x01 = this->spheres[lowerShell]->getRotatedMappedData ( lowerLatL * this->spheres[lowerShell]->getLocalAngRes() + upperLonL );
1183 x10 = this->spheres[lowerShell]->getRotatedMappedData ( upperLatL * this->spheres[lowerShell]->getLocalAngRes() + lowerLonL );
1184 x11 = this->spheres[lowerShell]->getRotatedMappedData ( upperLatL * this->spheres[lowerShell]->getLocalAngRes() + upperLonL );
1186 distLLon = std::abs ( lon - lonCOL.at(lowerLonL) ) / ( std::abs( lon - lonCOL.at(lowerLonL) ) + std::abs( lon - lonCOL.at(upperLonL) ) );
1187 valLLon = ( ( 1.0 - distLLon ) * x00 ) + ( distLLon * x01 );
1188 valULon = ( ( 1.0 - distLLon ) * x10 ) + ( distLLon * x11 );
1190 distLLat = std::abs ( lat - latCOL.at(lowerLatL) ) / ( std::abs( lat - latCOL.at(lowerLatL) ) + std::abs( lat - latCOL.at(upperLatL) ) );
1191 lowerShellValue = ( ( 1.0 - distLLat ) * valLLon ) + ( distLLat * valULon );
1194 x00 = this->spheres[upperShell]->getRotatedMappedData ( lowerLatU * this->spheres[upperShell]->getLocalAngRes() + lowerLonU );
1195 x01 = this->spheres[upperShell]->getRotatedMappedData ( lowerLatU * this->spheres[upperShell]->getLocalAngRes() + upperLonU );
1196 x10 = this->spheres[upperShell]->getRotatedMappedData ( upperLatU * this->spheres[upperShell]->getLocalAngRes() + lowerLonU );
1197 x11 = this->spheres[upperShell]->getRotatedMappedData ( upperLatU * this->spheres[upperShell]->getLocalAngRes() + upperLonU );
1199 distLLon = std::abs ( lon - lonCOU.at(lowerLonU) ) / ( std::abs( lon - lonCOU.at(lowerLonU) ) + std::abs( lon - lonCOU.at(upperLonU) ) );
1200 valLLon = ( ( 1.0 - distLLon ) * x00 ) + ( distLLon * x01 );
1201 valULon = ( ( 1.0 - distLLon ) * x10 ) + ( distLLon * x11 );
1203 distLLat = std::abs ( lat - latCOU.at(lowerLatU) ) / ( std::abs( lat - latCOU.at(lowerLatU) ) + std::abs( lat - latCOU.at(upperLatU) ) );
1204 upperShellValue = ( ( 1.0 - distLLat ) * valLLon ) + ( distLLat * valULon );
1207 distLRad = std::abs ( rad - this->spherePos.at(lowerShell) ) / ( std::abs( rad - this->spherePos.at(lowerShell) ) +
1208 std::abs( rad - this->spherePos.at(upperShell) ) );
1210 arrPos = wIt + this->zDimIndices * ( vIt + this->yDimIndices * uIt );
1211 densityMapRotated[arrPos] = ( ( 1.0 - distLRad ) * lowerShellValue ) + ( distLRad * upperShellValue );
1231 std::vector < proshade_double > ret;
1232 proshade_double eulA, eulB, eulG;
1236 this->getMaxBand() * 2,
1237 &eulA, &eulB, &eulG, settings );