diff --git a/Src/Kernels/FElecForcesKernels.hpp b/Src/Kernels/FElecForcesKernels.hpp index 5b165ce95a076cc3b083726b0355faa7a82274e3..d87520dbdbed7d1f10f5c8d963a249e409f42973 100644 --- a/Src/Kernels/FElecForcesKernels.hpp +++ b/Src/Kernels/FElecForcesKernels.hpp @@ -111,7 +111,7 @@ class FElecForcesKernels : public FAbstractKernels<ParticleClass,CellClass,Conta preM2MTransitionsPer = new FComplexe[periodicLevels * 8 * harmonic.getExpSize()]; FReal treeWidthAtLevel = boxWidth; - for(int idxLevel = -1 ; idxLevel <= -periodicLevels ; --idxLevel ){ + for(int idxLevel = -1 ; idxLevel >= -periodicLevels ; --idxLevel ){ const F3DPosition father(treeWidthAtLevel,treeWidthAtLevel,treeWidthAtLevel); treeWidthAtLevel *= 2; @@ -143,7 +143,7 @@ class FElecForcesKernels : public FAbstractKernels<ParticleClass,CellClass,Conta // so 6 in each dimension treeWidthAtLevel = boxWidth*2; preM2LTransitionsPer = new FComplexe[periodicLevels * (7 * 7 * 7) * devM2lP]; - for(int idxLevel = -1 ; idxLevel <= -periodicLevels ; --idxLevel ){ + for(int idxLevel = -1 ; idxLevel >= -periodicLevels ; --idxLevel ){ for(int idxX = -3 ; idxX <= 3 ; ++idxX ){ for(int idxY = -3 ; idxY <= 3 ; ++idxY ){ for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){ @@ -870,10 +870,11 @@ private: // We want: - gradient(POTENTIAL_SIGN potential). // The -(- 1.0) computing is not the most efficient programming ... const FReal signe = 1.0; - force_vector_in_local_base_x = ( force_vector_in_local_base_x * signe / spherical.getR()); - force_vector_in_local_base_y = ( force_vector_in_local_base_y * signe / spherical.getR()); - force_vector_in_local_base_z = ( force_vector_in_local_base_z * signe / (spherical.getR() * spherical.getSinTheta())); - + if( FMath::Epsilon < spherical.getR()){ + force_vector_in_local_base_x = ( force_vector_in_local_base_x * signe / spherical.getR()); + force_vector_in_local_base_y = ( force_vector_in_local_base_y * signe / spherical.getR()); + force_vector_in_local_base_z = ( force_vector_in_local_base_z * signe / (spherical.getR() * spherical.getSinTheta())); + } ///////////////////////////////////////////////////////////////////// //spherical_position_Set_ph diff --git a/Src/Utils/FSpherical.hpp b/Src/Utils/FSpherical.hpp index 4cf0b5bbcbc0288052a09a6878e9ff183b5458a7..cf591c62021b8730fbe87c168c4970c737b211d2 100644 --- a/Src/Utils/FSpherical.hpp +++ b/Src/Utils/FSpherical.hpp @@ -2,6 +2,7 @@ #define FSPHERICAL_HPP #include "FGlobal.hpp" +#include "FMath.hpp" #include "F3DPosition.hpp" /** This class is a Spherical position @@ -20,8 +21,14 @@ public: const FReal x2y2 = (inVector.getX() * inVector.getX()) + (inVector.getY() * inVector.getY()); this->r = FMath::Sqrt( x2y2 + (inVector.getZ() * inVector.getZ())); this->phi = FMath::Atan2(inVector.getY(),inVector.getX()); - this->cosTheta = inVector.getZ() / r; - this->sinTheta = FMath::Sqrt(x2y2) / r; + if( r < FMath::Epsilon ){ // if r == 0 we cannot divide! + this->cosTheta = FReal(1); + this->sinTheta = FReal(1); + } + else { + this->cosTheta = inVector.getZ() / r; + this->sinTheta = FMath::Sqrt(x2y2) / r; + } } /** Get the rayon */