98 eulA, eulB, eulG, settings );
142 std::max ( staticStructure->
getYDim(), movingStructure->
getYDim() ),
143 std::max ( staticStructure->
getZDim(), movingStructure->
getZDim() ) );
145 std::max ( staticStructure->
getYDim(), movingStructure->
getYDim() ),
146 std::max ( staticStructure->
getZDim(), movingStructure->
getZDim() ) );
155 proshade_double mapPeak = 0.0;
156 proshade_unsign xDimS = staticStructure->
getXDim();
157 proshade_unsign yDimS = staticStructure->
getYDim();
158 proshade_unsign zDimS = staticStructure->
getZDim();
162 if ( *trsX > (
static_cast< proshade_double
> ( xDimS ) / 2.0 ) ) { *trsX = *trsX -
static_cast< proshade_double
> ( xDimS ); }
163 if ( *trsY > (
static_cast< proshade_double
> ( yDimS ) / 2.0 ) ) { *trsY = *trsY -
static_cast< proshade_double
> ( yDimS ); }
164 if ( *trsZ > (
static_cast< proshade_double
> ( zDimS ) / 2.0 ) ) { *trsZ = *trsZ -
static_cast< proshade_double
> ( zDimS ); }
173 proshade_single xMov =
static_cast< proshade_single
> ( -(*trsX) );
174 proshade_single yMov =
static_cast< proshade_single
> ( -(*trsY) );
175 proshade_single zMov =
static_cast< proshade_single
> ( -(*trsZ) );
178 std::stringstream hlpSS;
179 hlpSS <<
"Optimal map translation distances are " << *trsX <<
" ; " << *trsY <<
" ; " << *trsZ <<
" Angstroms with peak height " << mapPeak / (
static_cast< proshade_double
> ( xDimS ) *
static_cast< proshade_double
> ( yDimS ) *
static_cast< proshade_double
> ( zDimS ) );
195 static_cast< proshade_signed
> ( movingStructure->
getXDim() ),
static_cast< proshade_signed
> ( movingStructure->
getYDim() ),
196 static_cast< proshade_signed
> ( movingStructure->
getZDim() ) );
223 proshade_single xSamplRate = (
static_cast< proshade_single
> ( staticStructure->
getXDimSize() ) /
static_cast< proshade_single
> ( staticStructure->
getXDim() ) );
224 proshade_single ySamplRate = (
static_cast< proshade_single
> ( staticStructure->
getYDimSize() ) /
static_cast< proshade_single
> ( staticStructure->
getYDim() ) );
225 proshade_single zSamplRate = (
static_cast< proshade_single
> ( staticStructure->
getZDimSize() ) /
static_cast< proshade_single
> ( staticStructure->
getZDim() ) );
227 proshade_single xCor =
static_cast< proshade_single
> (
static_cast< proshade_single
> ( staticStructure->
xFrom - movingStructure->
xFrom ) * xSamplRate );
228 proshade_single yCor =
static_cast< proshade_single
> (
static_cast< proshade_single
> ( staticStructure->
yFrom - movingStructure->
yFrom ) * ySamplRate );
229 proshade_single zCor =
static_cast< proshade_single
> (
static_cast< proshade_single
> ( staticStructure->
zFrom - movingStructure->
zFrom ) * zSamplRate );
231 proshade_single xMov =
static_cast< proshade_single
> ( 0.0 -
static_cast< proshade_double
> ( xCor ) - ( *trsX *
static_cast< proshade_double
> ( xSamplRate ) ) );
232 proshade_single yMov =
static_cast< proshade_single
> ( 0.0 -
static_cast< proshade_double
> ( yCor ) - ( *trsY *
static_cast< proshade_double
> ( ySamplRate ) ) );
233 proshade_single zMov =
static_cast< proshade_single
> ( 0.0 -
static_cast< proshade_double
> ( zCor ) - ( *trsZ *
static_cast< proshade_double
> ( zSamplRate ) ) );
236 *trsX =
static_cast< proshade_double
> ( -xMov );
237 *trsY =
static_cast< proshade_double
> ( -yMov );
238 *trsZ =
static_cast< proshade_double
> ( -zMov );
253 std::vector< proshade_double > ret;
254 proshade_double mapPeak = 0.0;
255 proshade_unsign xDimS = staticStructure->
getXDim();
256 proshade_unsign yDimS = staticStructure->
getYDim();
257 proshade_unsign zDimS = staticStructure->
getZDim();
258 proshade_double trsX = 0.0, trsY = 0.0, trsZ = 0.0;
269 if ( trsX > (
static_cast< proshade_double
> ( xDimS ) / 2.0 ) ) { trsX = trsX -
static_cast< proshade_double
> ( xDimS ); }
270 if ( trsY > (
static_cast< proshade_double
> ( yDimS ) / 2.0 ) ) { trsY = trsY -
static_cast< proshade_double
> ( yDimS ); }
271 if ( trsZ > (
static_cast< proshade_double
> ( zDimS ) / 2.0 ) ) { trsZ = trsZ -
static_cast< proshade_double
> ( zDimS ); }
276 this->originalPdbTransX = trsX;
277 this->originalPdbTransY = trsY;
278 this->originalPdbTransZ = trsZ;
280 proshade_single xMov =
static_cast< proshade_single
> ( -trsX );
281 proshade_single yMov =
static_cast< proshade_single
> ( -trsY );
282 proshade_single zMov =
static_cast< proshade_single
> ( -trsZ );
290 this->mapMovFromsChangeX =
static_cast<proshade_double
> ( this->xFrom );
291 this->mapMovFromsChangeY =
static_cast<proshade_double
> ( this->yFrom );
292 this->mapMovFromsChangeZ =
static_cast<proshade_double
> ( this->zFrom );
296 this->getXFromPtr(), this->getXToPtr(),
297 this->getYFromPtr(), this->getYToPtr(),
298 this->getZFromPtr(), this->getZToPtr(),
299 this->getXAxisOrigin(), this->getYAxisOrigin(), this->getZAxisOrigin() );
302 this->getXDimSize(), this->getYDimSize(), this->getZDimSize(),
303 static_cast< proshade_signed
> ( this->getXDim() ),
static_cast< proshade_signed
> ( this->getYDim() ),
static_cast< proshade_signed
> ( this->getZDim() ) );
306 this->mapMovFromsChangeX =
static_cast<proshade_double
> ( this->xFrom ) - this->mapMovFromsChangeX;
307 this->mapMovFromsChangeY =
static_cast<proshade_double
> ( this->yFrom ) - this->mapMovFromsChangeY;
308 this->mapMovFromsChangeZ =
static_cast<proshade_double
> ( this->zFrom ) - this->mapMovFromsChangeZ;
327 fftw_complex *tmpIn1 =
nullptr, *tmpOut1 =
nullptr, *tmpIn2 =
nullptr, *tmpOut2 =
nullptr, *resOut =
nullptr;
328 fftw_plan forwardFourierObj1, forwardFourierObj2, inverseFourierCombo;
329 proshade_unsign dimMult = staticStructure->
getXDim() * staticStructure->
getYDim() * staticStructure->
getZDim();
330 ProSHADE_internal_overlay::allocateTranslationFunctionMemory ( tmpIn1, tmpOut1, tmpIn2, tmpOut2, this->translationMap, resOut, forwardFourierObj1, forwardFourierObj2, inverseFourierCombo, staticStructure->
getXDim(), staticStructure->
getYDim(), staticStructure->
getZDim() );
333 for ( proshade_unsign iter = 0; iter < dimMult; iter++ ) { tmpIn1[iter][0] = staticStructure->
getMapValue ( iter ); tmpIn1[iter][1] = 0.0; }
334 for ( proshade_unsign iter = 0; iter < dimMult; iter++ ) { tmpIn2[iter][0] = this->getMapValue ( iter ); tmpIn2[iter][1] = 0.0; }
337 fftw_execute ( forwardFourierObj1 );
338 fftw_execute ( forwardFourierObj2 );
342 fftw_execute ( inverseFourierCombo );
367 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 )
370 tmpIn1 =
new fftw_complex[xD * yD * zD];
371 tmpOut1 =
new fftw_complex[xD * yD * zD];
372 tmpIn2 =
new fftw_complex[xD * yD * zD];
373 tmpOut2 =
new fftw_complex[xD * yD * zD];
374 resIn =
new fftw_complex[xD * yD * zD];
375 resOut =
new fftw_complex[xD * yD * zD];
386 forwardFourierObj1 = fftw_plan_dft_3d (
static_cast< int > ( xD ),
static_cast< int > ( yD ),
static_cast< int > ( zD ), tmpIn1, tmpOut1, FFTW_FORWARD , FFTW_ESTIMATE );
387 forwardFourierObj2 = fftw_plan_dft_3d (
static_cast< int > ( xD ),
static_cast< int > ( yD ),
static_cast< int > ( zD ), tmpIn2, tmpOut2, FFTW_FORWARD , FFTW_ESTIMATE );
388 inverseFourierCombo = fftw_plan_dft_3d (
static_cast< int > ( xD ),
static_cast< int > ( yD ),
static_cast< int > ( zD ), resOut, resIn , FFTW_BACKWARD, FFTW_ESTIMATE );
409 fftw_destroy_plan ( forwardFourierObj1 );
410 fftw_destroy_plan ( forwardFourierObj2 );
411 fftw_destroy_plan ( inverseFourierCombo );
435 double normFactor =
static_cast<double> ( xD * yD * zD );
436 proshade_signed arrPos;
439 for ( proshade_signed xIt = 0; xIt < static_cast< proshade_signed > ( xD ); xIt++ )
441 for ( proshade_signed yIt = 0; yIt < static_cast< proshade_signed > ( yD ); yIt++ )
443 for ( proshade_signed zIt = 0; zIt < static_cast< proshade_signed > ( zD ); zIt++ )
446 arrPos = zIt +
static_cast< proshade_signed
> ( zD ) * ( yIt +
static_cast< proshade_signed
> ( yD ) * xIt );
454 &resOut[arrPos][1] );
457 resOut[arrPos][0] /= normFactor;
458 resOut[arrPos][1] /= normFactor;
482 proshade_signed arrPos;
486 for ( proshade_signed uIt = 0; uIt < static_cast<proshade_signed> ( xD ); uIt++ )
488 for ( proshade_signed vIt = 0; vIt < static_cast<proshade_signed> ( yD ); vIt++ )
490 for ( proshade_signed wIt = 0; wIt < static_cast<proshade_signed> ( zD ); wIt++ )
492 arrPos = wIt +
static_cast< proshade_signed
> ( zD ) * ( vIt +
static_cast< proshade_signed
> ( yD ) * uIt );
493 if ( resIn[arrPos][0] > *mapPeak )
495 *mapPeak = resIn[arrPos][0];
496 *trsX =
static_cast< proshade_double
> ( uIt );
497 *trsY =
static_cast< proshade_double
> ( vIt );
498 *trsZ =
static_cast< proshade_double
> ( wIt );
522 if ( ( this->xDimIndices > xDim ) || ( this->yDimIndices > yDim ) || ( this->zDimIndices > zDim ) )
524 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." );
528 if ( ( this->xDimIndices == xDim ) && ( this->yDimIndices == yDim ) && ( this->zDimIndices == zDim ) ) { return ; }
531 proshade_unsign addXPre, addYPre, addZPre, addXPost, addYPost, addZPost;
532 ProSHADE_internal_overlay::computeBeforeAfterZeroCounts ( &addXPre, &addYPre, &addZPre, &addXPost, &addYPost, &addZPost, xDim, yDim, zDim, this->xDimIndices, this->yDimIndices, this->zDimIndices );
535 proshade_double* newMap =
new proshade_double [xDim * yDim * zDim];
538 ProSHADE_internal_overlay::paddMapWithZeroes ( this->internalMap, newMap, xDim, yDim, zDim, this->xDimIndices, this->yDimIndices, this->zDimIndices, addXPre, addYPre, addZPre );
541 delete[] this->internalMap;
542 this->internalMap =
new proshade_double [xDim * yDim * zDim];
543 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { this->internalMap[iter] = newMap[iter]; }
549 this->xDimSize =
static_cast< proshade_single
> ( xDim ) * ( this->xDimSize /
static_cast< proshade_single
> ( this->xDimIndices ) );
550 this->yDimSize =
static_cast< proshade_single
> ( yDim ) * ( this->yDimSize /
static_cast< proshade_single
> ( this->yDimIndices ) );
551 this->zDimSize =
static_cast< proshade_single
> ( zDim ) * ( this->zDimSize /
static_cast< proshade_single
> ( this->zDimIndices ) );
552 this->xDimIndices = xDim ; this->yDimIndices = yDim ; this->zDimIndices = zDim;
553 this->xGridIndices = xDim ; this->yGridIndices = yDim ; this->zGridIndices = zDim;
554 this->xFrom -= addXPre ; this->yFrom -= addYPre ; this->zFrom -= addZPre;
555 this->xTo += addXPost; this->yTo += addYPost; this->zTo += addZPost;
556 this->xAxisOrigin -= addXPre ; this->yAxisOrigin -= addYPre ; this->zAxisOrigin -= addZPre ;
578 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 )
581 *addXPre = ( xDim - xDimIndices ) / 2;
582 *addYPre = ( yDim - yDimIndices ) / 2;
583 *addZPre = ( zDim - zDimIndices ) / 2;
584 *addXPost = *addXPre;
if ( ( xDim - xDimIndices ) % 2 == 1 ) { *addXPost += 1; }
585 *addYPost = *addYPre;
if ( ( yDim - yDimIndices ) % 2 == 1 ) { *addYPost += 1; }
586 *addZPost = *addZPre;
if ( ( zDim - zDimIndices ) % 2 == 1 ) { *addZPost += 1; }
607 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 )
610 proshade_unsign newMapIndex = 0;
611 proshade_unsign oldMapIndex = 0;
614 for ( proshade_unsign xIt = 0; xIt < xDim; xIt++ )
616 for ( proshade_unsign yIt = 0; yIt < yDim; yIt++ )
618 for ( proshade_unsign zIt = 0; zIt < zDim; zIt++ )
621 newMapIndex = zIt + zDim * ( yIt + yDim * xIt );
624 if ( xIt < addXPre ) { newMap[newMapIndex] = 0.0;
continue; }
625 if ( yIt < addYPre ) { newMap[newMapIndex] = 0.0;
continue; }
626 if ( zIt < addZPre ) { newMap[newMapIndex] = 0.0;
continue; }
629 if ( xIt >= ( addXPre + xDimIndices ) ) { newMap[newMapIndex] = 0.0;
continue; }
630 if ( yIt >= ( addYPre + yDimIndices ) ) { newMap[newMapIndex] = 0.0;
continue; }
631 if ( zIt >= ( addZPre + zDimIndices ) ) { newMap[newMapIndex] = 0.0;
continue; }
634 oldMapIndex = (zIt-addZPre) + zDimIndices * ( (yIt-addYPre) + yDimIndices * (xIt-addXPre) );
635 newMap[newMapIndex] = oldMap[oldMapIndex];
660 this->maxCompBand = this->spheres[this->noSpheres-1]->getLocalBandwidth();
663 this->findMapCOM ( );
664 this->mapCOMProcessChangeX += ( this->xCom - this->originalMapXCom );
665 this->mapCOMProcessChangeY += ( this->yCom - this->originalMapYCom );
666 this->mapCOMProcessChangeZ += ( this->zCom - this->originalMapZCom );
672 this->allocateRotatedSHMemory ( );
675 this->computeRotatedSH ( );
678 this->invertSHCoefficients ( );
681 std::vector<proshade_double> lonCO, latCO;
685 proshade_double *densityMapRotated =
new proshade_double [this->xDimIndices * this->yDimIndices * this->zDimIndices];
687 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { densityMapRotated[iter] = 0.0; }
690 this->interpolateMapFromSpheres ( densityMapRotated );
693 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
695 this->internalMap[iter] = densityMapRotated[iter];
699 delete[] densityMapRotated;
721 bool withinBounds =
true;
722 proshade_double c000, c001, c010, c011, c100, c101, c110, c111, c00, c01, c10, c11, c0, c1;
724 proshade_double xCOM, yCOM, zCOM;
725 std::vector< proshade_double > ret;
728 proshade_single xSampRate = this->xDimSize /
static_cast< proshade_single
> ( this->xTo - this->xFrom );
729 proshade_single ySampRate = this->yDimSize /
static_cast< proshade_single
> ( this->yTo - this->yFrom );
730 proshade_single zSampRate = this->zDimSize /
static_cast< proshade_single
> ( this->zTo - this->zFrom );
733 ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xCOM, &yCOM, &zCOM, this->xDimSize, this->yDimSize, this->zDimSize, this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
736 proshade_single *mins =
new proshade_single[3];
737 proshade_single *maxs =
new proshade_single[3];
738 proshade_single *rotMat =
new proshade_single[9];
739 proshade_single *rotVec;
740 proshade_single *interpMins =
new proshade_single[3];
741 proshade_single *interpMaxs =
new proshade_single[3];
742 proshade_single *interpDiff =
new proshade_single[3];
743 proshade_single *movs =
new proshade_single[3];
744 map =
new proshade_double[ this->xDimIndices * this->yDimIndices * this->zDimIndices ];
757 for (
size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { map[iter] = 0.0; }
760 mins[0] = std::floor (
static_cast< proshade_single
> ( this->xDimIndices ) / -2.0f );
761 mins[1] = std::floor (
static_cast< proshade_single
> ( this->yDimIndices ) / -2.0f );
762 mins[2] = std::floor (
static_cast< proshade_single
> ( this->zDimIndices ) / -2.0f );
768 if ( this->xDimIndices % 2 == 0 ) { maxs[0] -= 1.0f; }
769 if ( this->yDimIndices % 2 == 0 ) { maxs[1] -= 1.0f; }
770 if ( this->zDimIndices % 2 == 0 ) { maxs[2] -= 1.0f; }
786 for ( proshade_single xIt = mins[0]; xIt <= maxs[0]; xIt += 1.0f )
788 for ( proshade_single yIt = mins[1]; yIt <= maxs[1]; yIt += 1.0f )
790 for ( proshade_single zIt = mins[2]; zIt <= maxs[2]; zIt += 1.0f )
797 for (
size_t posIt = 0; posIt < 3; posIt++ )
800 interpMins[posIt] = std::floor ( rotVec[posIt] );
801 interpMaxs[posIt] = interpMins[posIt] + 1.0f;
804 if ( ( maxs[posIt] < interpMins[posIt] ) || ( interpMins[posIt] < mins[posIt] ) || ( maxs[posIt] < interpMaxs[posIt] ) || ( interpMaxs[posIt] < mins[posIt] ) )
806 withinBounds =
false;
811 interpDiff[posIt] = rotVec[posIt] - interpMins[posIt];
813 if ( !withinBounds ) {
continue; }
816 for (
size_t posIt = 0; posIt < 3; posIt++ )
818 interpMaxs[posIt] = std::min ( maxs[posIt], std::max ( mins[posIt], interpMaxs[posIt] ) );
825 arrPos =
static_cast< size_t > ( ( interpMins[2] - mins[2] ) +
static_cast< proshade_single
> ( this->zDimIndices ) * ( ( interpMins[1] - mins[1] ) +
static_cast< proshade_single
> ( this->yDimIndices ) * ( interpMins[0] - mins[0] ) ) );
826 c000 = this->internalMap[arrPos];
828 arrPos =
static_cast< size_t > ( ( interpMaxs[2] - mins[2] ) +
static_cast< proshade_single
> ( this->zDimIndices ) * ( ( interpMins[1] - mins[1] ) +
static_cast< proshade_single
> ( this->yDimIndices ) * ( interpMins[0] - mins[0] ) ) );
829 c001 = this->internalMap[arrPos];
831 arrPos =
static_cast< size_t > ( ( interpMins[2] - mins[2] ) +
static_cast< proshade_single
> ( this->zDimIndices ) * ( ( interpMaxs[1] - mins[1] ) +
static_cast< proshade_single
> ( this->yDimIndices ) * ( interpMins[0] - mins[0] ) ) );
832 c010 = this->internalMap[arrPos];
834 arrPos =
static_cast< size_t > ( ( interpMaxs[2] - mins[2] ) +
static_cast< proshade_single
> ( this->zDimIndices ) * ( ( interpMaxs[1] - mins[1] ) +
static_cast< proshade_single
> ( this->yDimIndices ) * ( interpMins[0] - mins[0] ) ) );
835 c011 = this->internalMap[arrPos];
837 arrPos =
static_cast< size_t > ( ( interpMins[2] - mins[2] ) +
static_cast< proshade_single
> ( this->zDimIndices ) * ( ( interpMins[1] - mins[1] ) +
static_cast< proshade_single
> ( this->yDimIndices ) * ( interpMaxs[0] - mins[0] ) ) );
838 c100 = this->internalMap[arrPos];
840 arrPos =
static_cast< size_t > ( ( interpMaxs[2] - mins[2] ) +
static_cast< proshade_single
> ( this->zDimIndices ) * ( ( interpMins[1] - mins[1] ) +
static_cast< proshade_single
> ( this->yDimIndices ) * ( interpMaxs[0] - mins[0] ) ) );
841 c101 = this->internalMap[arrPos];
843 arrPos =
static_cast< size_t > ( ( interpMins[2] - mins[2] ) +
static_cast< proshade_single
> ( this->zDimIndices ) * ( ( interpMaxs[1] - mins[1] ) +
static_cast< proshade_single
> ( this->yDimIndices ) * ( interpMaxs[0] - mins[0] ) ) );
844 c110 = this->internalMap[arrPos];
846 arrPos =
static_cast< size_t > ( ( interpMaxs[2] - mins[2] ) +
static_cast< proshade_single
> ( this->zDimIndices ) * ( ( interpMaxs[1] - mins[1] ) +
static_cast< proshade_single
> ( this->yDimIndices ) * ( interpMaxs[0] - mins[0] ) ) );
847 c111 = this->internalMap[arrPos];
850 c00 = ( c000 * ( 1.0 -
static_cast< proshade_double
> ( interpDiff[0] ) ) ) + ( c100 *
static_cast< proshade_double
> ( interpDiff[0] ) );
851 c01 = ( c001 * ( 1.0 -
static_cast< proshade_double
> ( interpDiff[0] ) ) ) + ( c101 *
static_cast< proshade_double
> ( interpDiff[0] ) );
852 c10 = ( c010 * ( 1.0 -
static_cast< proshade_double
> ( interpDiff[0] ) ) ) + ( c110 *
static_cast< proshade_double
> ( interpDiff[0] ) );
853 c11 = ( c011 * ( 1.0 -
static_cast< proshade_double
> ( interpDiff[0] ) ) ) + ( c111 *
static_cast< proshade_double
> ( interpDiff[0] ) );
856 c0 = ( c00 * ( 1.0 -
static_cast< proshade_double
> ( interpDiff[1] ) ) ) + ( c10 *
static_cast< proshade_double
> ( interpDiff[1] ) );
857 c1 = ( c01 * ( 1.0 -
static_cast< proshade_double
> ( interpDiff[1] ) ) ) + ( c11 *
static_cast< proshade_double
> ( interpDiff[1] ) );
860 arrPos =
static_cast< size_t > ( ( zIt - mins[2] ) +
static_cast< proshade_single
> ( this->zDimIndices ) * ( ( yIt - mins[1] ) +
static_cast< proshade_single
> ( this->yDimIndices ) * ( xIt - mins[0] ) ) );
861 map[arrPos] = ( c0 * ( 1.0 -
static_cast< proshade_double
> ( interpDiff[2] ) ) ) + ( c1 *
static_cast< proshade_double
> ( interpDiff[2] ) );
894 proshade_double axX, axY, axZ, axAng, tmp, *rMat, *map;
897 rMat =
new proshade_double[9];
922 std::vector< proshade_double > ret = this->rotateMapRealSpace ( axX, axY, axZ, axAng, map );
925 for (
size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
927 this->internalMap[iter] = map[iter];
935 this->originalPdbRotCenX = ret.at(0);
936 this->originalPdbRotCenY = ret.at(1);
937 this->originalPdbRotCenZ = ret.at(2);
957 proshade_single xMov =
static_cast< proshade_single
> ( -trsX );
958 proshade_single yMov =
static_cast< proshade_single
> ( -trsY );
959 proshade_single zMov =
static_cast< proshade_single
> ( -trsZ );
963 this->getXFromPtr(), this->getXToPtr(), this->getYFromPtr(), this->getYToPtr(),
964 this->getZFromPtr(), this->getZToPtr(), this->getXAxisOrigin(), this->getYAxisOrigin(), this->getZAxisOrigin() );
968 static_cast< proshade_signed
> ( this->getXDim() ),
static_cast< proshade_signed
> ( this->getYDim() ),
969 static_cast< proshade_signed
> ( this->getZDim() ) );
981 this->rotSphericalHarmonics =
new proshade_complex* [this->noSpheres];
985 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
988 this->rotSphericalHarmonics[iter] =
new proshade_complex [
static_cast<proshade_unsign
> ( pow ( (this->spheres[iter]->getLocalBandwidth() * 2), 2) )];
992 for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( pow ( (this->spheres[iter]->getLocalBandwidth() * 2), 2) ); it++ )
994 this->rotSphericalHarmonics[iter][it][0] = 0.0;
995 this->rotSphericalHarmonics[iter][it][1] = 0.0;
1009 proshade_double WigDR, WigDI, *ShR, *ShI, retR, retI;
1010 proshade_unsign arrPos;
1013 for ( proshade_signed shell = 0; shell < static_cast<proshade_signed> ( this->noSpheres ); shell++ )
1016 for ( proshade_signed bandIter = 0; bandIter < static_cast<proshade_signed> ( this->spheres[shell]->getLocalBandwidth() ); bandIter++ )
1019 for ( proshade_signed order1 = 0; order1 < ( ( bandIter * 2 ) + 1 ); order1++ )
1022 ShR = getRealSphHarmValue (
static_cast< proshade_unsign
> ( bandIter ),
static_cast< proshade_unsign
> ( order1 ),
static_cast< proshade_unsign
> ( shell ) );
1023 ShI = getImagSphHarmValue (
static_cast< proshade_unsign
> ( bandIter ),
static_cast< proshade_unsign
> ( order1 ),
static_cast< proshade_unsign
> ( shell ) );
1026 for ( proshade_signed order2 = 0; order2 < ( ( bandIter * 2 ) + 1 ); order2++ )
1029 this->getWignerMatrixValue (
static_cast< proshade_unsign
> ( bandIter ),
static_cast< proshade_unsign
> ( order1 ),
static_cast< proshade_unsign
> ( order2 ), &WigDR, &WigDI );
1035 arrPos =
static_cast<proshade_unsign
> ( seanindex (
static_cast< int > ( order2-bandIter ),
static_cast< int > ( bandIter ),
1036 static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ) ) );
1037 this->rotSphericalHarmonics[shell][arrPos][0] += retR;
1038 this->rotSphericalHarmonics[shell][arrPos][1] += retI;
1064 proshade_unsign oneDim = shBand * 2;
1067 sigR =
new double [(oneDim * oneDim * 4)];
1068 sigI = sigR + (oneDim * oneDim * 2);
1069 rcoeffs =
new double [(oneDim * oneDim * 2)];
1070 icoeffs = rcoeffs + (oneDim * oneDim);
1071 weights =
new double [4 * shBand];
1072 workspace =
new double [( 20 * shBand * shBand ) + ( 24 * shBand )];
1081 idctPlan = fftw_plan_r2r_1d (
static_cast< int > ( oneDim ), weights, workspace, FFTW_REDFT01, FFTW_ESTIMATE );
1084 int rank, howmany_rank;
1085 fftw_iodim dims[1], howmany_dims[1];
1088 dims[0].n = 2 *
static_cast< int > ( shBand );
1089 dims[0].is = 2 *
static_cast< int > ( shBand );
1092 howmany_dims[0].n = 2 *
static_cast< int > ( shBand );
1093 howmany_dims[0].is = 1;
1094 howmany_dims[0].os = 2 *
static_cast< int > ( shBand );
1097 ifftPlan = fftw_plan_guru_split_dft ( rank, dims, howmany_rank, howmany_dims, sigR, sigI, rcoeffs, icoeffs, FFTW_ESTIMATE );
1110 double *sigR =
nullptr, *sigI =
nullptr, *rcoeffs =
nullptr, *icoeffs =
nullptr, *weights =
nullptr, *workspace =
nullptr;
1111 fftw_plan idctPlan, ifftPlan;
1114 for (
int shell = 0; shell < static_cast<int> ( this->noSpheres ); shell++ )
1117 proshade_unsign oneDim = this->spheres[shell]->getLocalBandwidth() * 2;
1123 makeweights (
static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ), weights );
1126 this->spheres[shell]->allocateRotatedMap ( );
1129 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( oneDim * oneDim ); iter++ )
1131 rcoeffs[iter] = this->rotSphericalHarmonics[shell][iter][0];
1132 icoeffs[iter] = this->rotSphericalHarmonics[shell][iter][1];
1138 InvFST_semi_fly ( rcoeffs,
1142 static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ),
1145 static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ),
1150 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( oneDim * oneDim ); iter++ )
1152 this->spheres[shell]->setRotatedMappedData ( iter, sigR[iter] );
1156 fftw_destroy_plan ( idctPlan );
1157 fftw_destroy_plan ( ifftPlan );
1180 for ( proshade_unsign iter = 0; iter <= angRes; iter++ )
1182 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 ) ) );
1186 for ( proshade_unsign iter = 0; iter <= angRes; iter++ )
1188 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 ) );
1203 proshade_double rad = 0.0, lon = 0.0, lat = 0.0, newU = 0.0, newV = 0.0, newW = 0.0;
1204 proshade_unsign lowerLonL = 0, upperLonL = 0, lowerLonU = 0, upperLonU = 0, lowerLatL = 0, upperLatL = 0, lowerLatU = 0, upperLatU = 0, lowerShell = 0, upperShell = 0;
1205 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;
1206 proshade_double lowerShellValue = 0.0, upperShellValue = 0.0;
1207 proshade_double xSamplingRate =
static_cast<proshade_double
> ( this->xDimSize ) /
static_cast<proshade_double
> ( this->xDimIndices );
1208 proshade_double ySamplingRate =
static_cast<proshade_double
> ( this->yDimSize ) /
static_cast<proshade_double
> ( this->yDimIndices );
1209 proshade_double zSamplingRate =
static_cast<proshade_double
> ( this->zDimSize ) /
static_cast<proshade_double
> ( this->zDimIndices );
1210 proshade_signed arrPos;
1211 std::vector<proshade_double> lonCOU, latCOU, lonCOL, latCOL;
1213 for ( proshade_signed uIt = 0; uIt < static_cast<proshade_signed> (this->xDimIndices); uIt++ )
1215 for ( proshade_signed vIt = 0; vIt < static_cast<proshade_signed> (this->yDimIndices); vIt++ )
1217 for ( proshade_signed wIt = 0; wIt < static_cast<proshade_signed> (this->zDimIndices); wIt++ )
1220 newU =
static_cast<proshade_double
> ( uIt - (
static_cast<proshade_signed
> (this->xDimIndices) / 2 ) );
1221 newV =
static_cast<proshade_double
> ( vIt - (
static_cast<proshade_signed
> (this->yDimIndices) / 2 ) );
1222 newW =
static_cast<proshade_double
> ( wIt - (
static_cast<proshade_signed
> (this->zDimIndices) / 2 ) );
1225 if ( ( newU == 0.0 ) && ( newV == 0.0 ) && ( newW == 0.0 ) )
1227 arrPos = wIt +
static_cast< proshade_signed
> ( this->zDimIndices ) * ( vIt +
static_cast< proshade_signed
> ( this->yDimIndices ) * uIt );
1228 densityMapRotated[arrPos] = this->internalMap[arrPos];
1233 rad = sqrt ( pow( ( newU * xSamplingRate ), 2.0 ) +
1234 pow( ( newV * ySamplingRate ), 2.0 ) +
1235 pow( ( newW * zSamplingRate ), 2.0 ) );
1236 lon = atan2 ( ( newV * ySamplingRate ), ( newU * xSamplingRate ) );
1237 lat = asin ( ( newW * zSamplingRate ) / rad );
1240 if ( rad != rad ) { rad = 0.0; }
1241 if ( lon != lon ) { lon = 0.0; }
1242 if ( lat != lat ) { lat = 0.0; }
1247 for ( proshade_unsign iter = 0; iter < (this->noSpheres-1); iter++ )
1249 if ( (
static_cast< proshade_double
> ( this->spherePos.at(iter) ) <= rad ) && (
static_cast< proshade_double
> ( this->spherePos.at(iter+1) ) > rad ) )
1252 upperShell = iter+1;
1257 if ( upperShell == 0 )
1259 arrPos = wIt +
static_cast< proshade_signed
> ( this->zDimIndices ) * ( vIt +
static_cast< proshade_signed
> ( this->yDimIndices ) * uIt );
1260 densityMapRotated[arrPos] = 0.0;
1265 lonCOL.clear(); latCOL.clear(); lonCOU.clear(); latCOU.clear();
1270 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( lonCOL.size() ); iter++ )
1272 if ( iter == (
static_cast<proshade_unsign
> ( lonCOL.size() ) - 1 ) )
1278 if ( ( std::floor(10000. * lonCOL.at(iter)) <= std::floor(10000. * lon) ) && ( std::floor(10000. * lonCOL.at(iter+1)) > std::floor(10000. * lon) ) )
1285 if ( upperLonL == this->spheres[lowerShell]->getLocalAngRes() ) { upperLonL = 0; }
1287 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( lonCOU.size() ); iter++ )
1289 if ( iter == (
static_cast<proshade_unsign
> ( lonCOU.size() ) - 1 ) )
1295 if ( ( std::floor(10000. * lonCOU.at(iter)) <= std::floor(10000. * lon) ) && ( std::floor(10000. * lonCOU.at(iter+1)) > std::floor(10000. * lon) ) )
1302 if ( upperLonU == this->spheres[upperShell]->getLocalAngRes() ) { upperLonU = 0; }
1304 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( latCOL.size() ); iter++ )
1306 if ( iter == (
static_cast<proshade_unsign
> ( latCOL.size() ) - 1 ) )
1312 if ( ( std::floor(10000. * latCOL.at(iter)) <= std::floor(10000. * lat) ) && ( std::floor(10000. * latCOL.at(iter+1)) > std::floor(10000. * lat) ) )
1319 if ( upperLatL == this->spheres[lowerShell]->getLocalAngRes() ) { upperLatL = 0; }
1321 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( latCOU.size() ); iter++ )
1323 if ( iter == (
static_cast<proshade_unsign
> ( latCOU.size() ) - 1 ) )
1329 if ( ( std::floor(10000. * latCOU.at(iter)) <= std::floor(10000. * lat) ) && ( std::floor(10000. * latCOU.at(iter+1)) > std::floor(10000. * lat) ) )
1336 if ( upperLatU == this->spheres[upperShell]->getLocalAngRes() ) { upperLatU = 0; }
1339 x00 = this->spheres[lowerShell]->getRotatedMappedData ( lowerLatL * this->spheres[lowerShell]->getLocalAngRes() + lowerLonL );
1340 x01 = this->spheres[lowerShell]->getRotatedMappedData ( lowerLatL * this->spheres[lowerShell]->getLocalAngRes() + upperLonL );
1341 x10 = this->spheres[lowerShell]->getRotatedMappedData ( upperLatL * this->spheres[lowerShell]->getLocalAngRes() + lowerLonL );
1342 x11 = this->spheres[lowerShell]->getRotatedMappedData ( upperLatL * this->spheres[lowerShell]->getLocalAngRes() + upperLonL );
1344 distLLon = std::abs ( lon - lonCOL.at(lowerLonL) ) / ( std::abs( lon - lonCOL.at(lowerLonL) ) + std::abs( lon - lonCOL.at(upperLonL) ) );
1345 valLLon = ( ( 1.0 - distLLon ) * x00 ) + ( distLLon * x01 );
1346 valULon = ( ( 1.0 - distLLon ) * x10 ) + ( distLLon * x11 );
1348 distLLat = std::abs ( lat - latCOL.at(lowerLatL) ) / ( std::abs( lat - latCOL.at(lowerLatL) ) + std::abs( lat - latCOL.at(upperLatL) ) );
1349 lowerShellValue = ( ( 1.0 - distLLat ) * valLLon ) + ( distLLat * valULon );
1352 x00 = this->spheres[upperShell]->getRotatedMappedData ( lowerLatU * this->spheres[upperShell]->getLocalAngRes() + lowerLonU );
1353 x01 = this->spheres[upperShell]->getRotatedMappedData ( lowerLatU * this->spheres[upperShell]->getLocalAngRes() + upperLonU );
1354 x10 = this->spheres[upperShell]->getRotatedMappedData ( upperLatU * this->spheres[upperShell]->getLocalAngRes() + lowerLonU );
1355 x11 = this->spheres[upperShell]->getRotatedMappedData ( upperLatU * this->spheres[upperShell]->getLocalAngRes() + upperLonU );
1357 distLLon = std::abs ( lon - lonCOU.at(lowerLonU) ) / ( std::abs( lon - lonCOU.at(lowerLonU) ) + std::abs( lon - lonCOU.at(upperLonU) ) );
1358 valLLon = ( ( 1.0 - distLLon ) * x00 ) + ( distLLon * x01 );
1359 valULon = ( ( 1.0 - distLLon ) * x10 ) + ( distLLon * x11 );
1361 distLLat = std::abs ( lat - latCOU.at(lowerLatU) ) / ( std::abs( lat - latCOU.at(lowerLatU) ) + std::abs( lat - latCOU.at(upperLatU) ) );
1362 upperShellValue = ( ( 1.0 - distLLat ) * valLLon ) + ( distLLat * valULon );
1365 distLRad = std::abs ( rad -
static_cast< proshade_double
> ( this->spherePos.at(lowerShell) ) ) / ( std::abs( rad -
static_cast< proshade_double
> ( this->spherePos.at(lowerShell) ) ) +
1366 std::abs( rad -
static_cast< proshade_double
> ( this->spherePos.at(upperShell) ) ) );
1368 arrPos = wIt +
static_cast< proshade_signed
> ( this->zDimIndices ) * ( vIt +
static_cast< proshade_signed
> ( this->yDimIndices ) * uIt );
1369 densityMapRotated[arrPos] = ( ( 1.0 - distLRad ) * lowerShellValue ) + ( distLRad * upperShellValue );
1389 std::vector < proshade_double > ret;
1390 proshade_double eulA, eulB, eulG;
1394 this->getMaxBand() * 2,
1395 &eulA, &eulB, &eulG, settings );