Commit 44b1f2ec authored by BRAMAS Berenger's avatar BRAMAS Berenger

Update Rotation kernel, still not working

parent 73b8cc2c
This diff is collapsed.
......@@ -16,46 +16,45 @@
#include "FPoint.hpp"
#include "FDebug.hpp"
/** This class is a Spherical position
* This class is currently in its minimum version
/**
* This class is a Spherical position
*
* @brief Spherical coordinate system
*
* The spherical coordinate system (r radius, theta inclination, phi azimuth) is the following <br>
* x = r sin(theta) cos(phi) <br>
* y = r sin(theta) sin(phi) <br>
* z = r cos(theta) <br>
* This system is defined in p 872 of the paper of Epton and Dembart, SIAM J Sci Comput 1995.<br>
*
* Even if it can look different from usual expression (where theta and phi are inversed),
* such expression is used to match the SH expression.
*/
//! Brief Spherical coordinate system
//! The spherical coordinate system (r radius, theta inclination, phi azimuth) is the following <br>
//! x = r sin(theta) cos(phi) <br>
//! y = r sin(theta) sin(phi) <br>
//! z = r cos(theta) <br>
//! This system is defined in p 872 of the paper of Epton and Dembart, SIAM J Sci Comput 1995.
//!
class FSpherical {
// The attributes of a sphere
FReal r; //!< the radius
FReal theta; //!< the inclination angle [0, pi] - colatitude
FReal phi; //!< the azimuth angle [-pi,pi] - longitude - around z axis
FReal cosTheta;
FReal sinTheta;
FReal theta; //!< the inclination angle
FReal phi; //!< the azimuth angle
public:
/** Default Constructor, set attributes to 0 */
FSpherical()
: r(0), cosTheta(0), sinTheta(0), theta(0), phi(0) {
: r(0), theta(0), phi(0), cosTheta(0), sinTheta(0) {
}
/** From now, we just need a constructor based on a 3D position */
explicit FSpherical(const FPoint& inVector){
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());
if( r < FMath::Epsilon ){ // if r == 0 we cannot divide!
FDEBUG( FDebug::Controller << "!!! In FSpherical, r == 0!\n"; )
this->cosTheta = FReal(1);
this->sinTheta = FReal(1);
}
else {
this->cosTheta = inVector.getZ() / r;
this->sinTheta = FMath::Sqrt(x2y2) / r;
}
this->cosTheta = inVector.getZ() / r;
this->sinTheta = FMath::Sqrt(x2y2) / r;
this->theta = FMath::ACos(this->cosTheta);
// if r == 0 we cannot divide!
FDEBUG(if( r < FMath::Epsilon ) FDebug::Controller << "!!! In FSpherical, r == 0!\n"; )
}
/** Get the radius */
......@@ -63,23 +62,33 @@ public:
return r;
}
/** Get the cosine of theta = z / r */
/** Get the inclination angle theta = acos(z/r) [0, pi] */
FReal getTheta() const{
return theta;
}
/** Get the azimuth angle phi = atan2(y,x) [-pi,pi] */
FReal getPhi() const{
return phi;
}
/** Get the inclination angle [0, pi] */
FReal getInclination() const{
return theta;
}
/** Get the azimuth angle [0,2pi] */
FReal getAzimuth() const{
return (phi < 0 ? FMath::FPi*2 + phi : phi);
}
/** Get the cos of theta = z / r */
FReal getCosTheta() const{
return cosTheta;
}
/** Get the sin theta = sqrt(x2y2) / r */
/** Get the sin of theta = sqrt(x2y2) / r */
FReal getSinTheta() const{
return sinTheta;
}
/** Get the inclination angle theta = acos(z/r) */
FReal getTheta() const{
return theta;
}
/** Get the azimuth angle phi = atan2(y,x) */
FReal getPhi() const{
return phi;
}
};
#endif // FSPHERICAL_HPP
......@@ -23,6 +23,7 @@
#include "../../Src/Kernels/Spherical/FSphericalParticle.hpp"
#include "../../Src/Kernels/Rotation/FRotationKernel.hpp"
#include "../../Src/Kernels/Rotation/FRotationCell.hpp"
#include "../../Src/Utils/FMath.hpp"
#include "../../Src/Utils/FMemUtils.hpp"
......@@ -52,7 +53,7 @@ public:
int main(int argc, char** argv){
static const int P = 2;
static const int P = 1;
typedef IndexedParticle ParticleClass;
typedef FRotationCell<P> CellClass;
......@@ -179,7 +180,7 @@ int main(int argc, char** argv){
const ParticleClass& other = particles[leafIter.data().getIndex()];
potentialDiff.add(other.getPotential(),leafIter.data().getPotential());
std::cout << "Direct Potential = " << other.getPotential() << " Fmm Potential = " << leafIter.data().getPotential() << std::endl; // Remove Me
std::cout << leafIter.data().getIndex() << " Direct Potential = " << other.getPotential() << " Fmm Potential = " << leafIter.data().getPotential() << std::endl; // Remove Me
fx.add(other.getForces().getX(),leafIter.data().getForces().getX());
fy.add(other.getForces().getY(),leafIter.data().getForces().getY());
......
......@@ -120,7 +120,8 @@ int main(){
//const FReal sinTheta = FMath::Sin(FMath::FPiDiv2/FReal(2.0));
const FPoint relativPos(FReal(0.1),FReal(0.1),FReal(0.1));
const FSpherical relativPosSphere(relativPos);
const FReal Betha = 2.186276;//FMath::FPi/2; //2*FMath::FPi-FMath::FPi/2
//Theta 0.955317 Phi 2.356194
const FReal Betha = 0.955317;//FMath::FPi/2; //2*FMath::FPi-FMath::FPi/2
const FReal cosTheta = FMath::Cos(Betha); // relativPosSphere.getCosTheta();
const FReal sinTheta = FMath::Sin(Betha); // relativPosSphere.getSinTheta();
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment