FSpherical.hpp 5.1 KB
Newer Older
1
// ===================================================================================
2
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
3 4 5 6 7 8 9 10 11 12 13 14
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.  
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info". 
// "http://www.gnu.org/licenses".
15
// ===================================================================================
berenger-bramas's avatar
berenger-bramas committed
16 17
#ifndef FSPHERICAL_HPP
#define FSPHERICAL_HPP
COULAUD Olivier's avatar
COULAUD Olivier committed
18
#include <iostream>
berenger-bramas's avatar
berenger-bramas committed
19 20

#include "FGlobal.hpp"
21
#include "FMath.hpp"
BRAMAS Berenger's avatar
BRAMAS Berenger committed
22
#include "FPoint.hpp"
BRAMAS Berenger's avatar
BRAMAS Berenger committed
23
#include "FLog.hpp"
berenger-bramas's avatar
berenger-bramas committed
24

25 26 27 28 29
/**
* This class is a Spherical position
*
* @brief Spherical coordinate system
*
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
* We consider the spherical coordinate system \f$(r, \theta, \varphi)\f$ commonly used in physics. r is the radial distance, \f$\theta\f$ the polar/inclination angle and \f$\varphi\f$ the azimuthal angle.<br>
* The <b>radial distance</b> is the Euclidean distance from the origin O to P.<br>
* The <b>inclination (or polar angle) </b>is the angle between the zenith direction and the line segment OP.<br>
* The <b>azimuth (or azimuthal angle)</b> is the signed angle measured from the azimuth reference direction to the orthogonal projection of the line segment OP on the reference plane.<br>
*
* The spherical coordinates of a point can be obtained from its Cartesian coordinates (x, y, z) by the formulae
* \f$ \displaystyle r = \sqrt{x^2 + y^2 + z^2}\f$<br>
* \f$ \displaystyle \theta = \displaystyle\arccos\left(\frac{z}{r}\right) \f$<br>
* \f$ \displaystyle \varphi = \displaystyle\arctan\left(\frac{y}{x}\right) \f$<br>
*and \f$\varphi\in[0,2\pi[ \f$ \f$ \theta\in[0,\pi]\f$<br>
*
* The spherical coordinate system  is retrieved from the the spherical coordinates by <br>
* \f$x = r \sin(\theta) \cos(\varphi)\f$ <br>
* \f$y = r \sin(\theta) \sin(\varphi)\f$ <br>
* \f$z = r \cos(\theta) \f$<br>
* with  \f$\varphi\in[-\pi,\pi[ \f$ \f$ \theta\in[0,\pi]\f$<br>
*
47 48 49 50
* 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.
51
*  See http://en.wikipedia.org/wiki/Spherical_coordinate_system
berenger-bramas's avatar
berenger-bramas committed
52 53
*/
class FSpherical {
54
    // The attributes of a sphere
55
    FReal r;         //!< the radial distance
COULAUD Olivier's avatar
COULAUD Olivier committed
56 57
    FReal theta;     //!< the inclination angle [0, pi] - colatitude, polar angle
    FReal phi;       //!< the azimuth angle [-pi,pi] - longitude - around z axis
58 59
    FReal cosTheta;
    FReal sinTheta;
berenger-bramas's avatar
berenger-bramas committed
60
public:
BRAMAS Berenger's avatar
BRAMAS Berenger committed
61 62
    /** Default Constructor, set attributes to 0 */
    FSpherical()
63
        : r(0.0), theta(0.0), phi(0.0), cosTheta(0.0), sinTheta(0.0) {
BRAMAS Berenger's avatar
BRAMAS Berenger committed
64 65
    }

66
    /** From now, we just need a constructor based on a 3D position */
BRAMAS Berenger's avatar
BRAMAS Berenger committed
67
    explicit FSpherical(const FPoint& inVector){
berenger-bramas's avatar
berenger-bramas committed
68
        const FReal x2y2 = (inVector.getX() * inVector.getX()) + (inVector.getY() * inVector.getY());
COULAUD Olivier's avatar
COULAUD Olivier committed
69
        this->r          = FMath::Sqrt( x2y2 + (inVector.getZ() * inVector.getZ()));
70

COULAUD Olivier's avatar
COULAUD Olivier committed
71
        this->phi        = FMath::Atan2(inVector.getY(),inVector.getX());
72 73 74

        this->cosTheta = inVector.getZ() / r;
        this->sinTheta = FMath::Sqrt(x2y2) / r;
COULAUD Olivier's avatar
COULAUD Olivier committed
75
        this->theta    = FMath::ACos(this->cosTheta);
76
        // if r == 0 we cannot divide!
BRAMAS Berenger's avatar
BRAMAS Berenger committed
77
        FLOG(if( r < FMath::Epsilon ) FLog::Controller << "!!! In FSpherical, r == 0!\n"; )
berenger-bramas's avatar
berenger-bramas committed
78 79
    }

COULAUD Olivier's avatar
COULAUD Olivier committed
80
    /** Get the radius */
berenger-bramas's avatar
berenger-bramas committed
81 82 83 84
    FReal getR() const{
        return r;
    }

85 86 87 88 89 90 91 92 93 94 95 96 97
    /** 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;
    }
98
    /** Get the azimuth angle [0,2pi]. You should use this method in order to obtain (x,y,z)*/
BRAMAS Berenger's avatar
BRAMAS Berenger committed
99
    FReal getPhiZero2Pi() const{
COULAUD Olivier's avatar
COULAUD Olivier committed
100
        return (phi < 0 ? FMath::FTwoPi + phi : phi);
101 102 103
    }

    /** Get the cos of theta = z / r */
berenger-bramas's avatar
berenger-bramas committed
104 105 106 107
    FReal getCosTheta() const{
        return cosTheta;
    }

108
    /** Get the sin of theta = sqrt(x2y2) / r */
berenger-bramas's avatar
berenger-bramas committed
109 110 111
    FReal getSinTheta() const{
        return sinTheta;
    }
COULAUD Olivier's avatar
COULAUD Olivier committed
112 113 114 115 116 117 118 119 120 121

    /**
     * Operator stream FPoint to std::ostream
     * This can be used to simpldata[1] write out a position
     * @param[in,out] output where to write the position
     * @param[in] inPosition the position to write out
     * @return the output for multiple << operators
     */
    template <class StreamClass>
    friend StreamClass& operator<<(StreamClass& output, const FSpherical& inPosition){
122
        output << "(" <<  inPosition.getR() << ", " << inPosition.getTheta() << ", " << inPosition.getPhi() <<")";
COULAUD Olivier's avatar
COULAUD Olivier committed
123 124 125 126
        return output;  // for multiple << operators.
    }


berenger-bramas's avatar
berenger-bramas committed
127 128 129
};

#endif // FSPHERICAL_HPP