MATSGN_Complex.h

Go to the documentation of this file.

#ifndef MATSGN_Complex_H
#define MATSGN_Complex_H

#include "CORE_Object.h"
#include "CORE_Exception.h"
#include "CORE_Real.h"
#include "fftw3.h"

DEFINE_SPTR(MATSGN_Complex);
class MATSGN_Complex : public virtual CORE_Object {
    SP_OBJECT(MATSGN_Complex);
    // ATTRIBUTES


private:
    fftw_complex mValue;
    mutable tDouble mWork;
    // ASSOCIATIONS
    
    
public:
    // METHODS
    
    // CONSTRUCTORS 
    
    MATSGN_Complex(void);  
    MATSGN_Complex(const tDouble& a,const tDouble& b);
    MATSGN_Complex(const fftw_complex& f);
    
    MATSGN_Complex(const MATSGN_Complex& f);
    

    
    
    // DESTRUCTORS 
    
    
    virtual ~MATSGN_Complex(void);   
    
public:
    // operators
    // =========
    inline tDouble& operator[](const tUSInt& i)  {
        ASSERT_IN(i<2);
        return mValue[i];
    };
    inline const tDouble& operator[](const tUSInt& i)  const {
        ASSERT_IN(i<2);
        return mValue[i];
    };
  
    inline MATSGN_Complex& operator=(const MATSGN_Complex&  s)  {
        mValue[0]=s[0];
        mValue[1]=s[1];
        return (*this);
    };
    inline MATSGN_Complex& operator=(const fftw_complex&  s)  {
        mValue[0]=s[0];
        mValue[1]=s[1];
        return (*this);
    };
   
    inline MATSGN_Complex& operator=(const tReal&  s)  {
        mValue[0]=s;
        mValue[1]=0;
        return (*this);
    };
   
    inline MATSGN_Complex& operator*=(const tReal&  s)  {
        mValue[0]*=s;
        mValue[1]*=s;
        return (*this);
    };
    inline MATSGN_Complex& operator*=(const MATSGN_Complex&  s)  {
        mWork=mValue[0]*s[0]-mValue[1]*s[1];
        mValue[1]=mValue[0]*s[1]+mValue[1]*s[0];
        mValue[0]=mWork;
        return (*this);
    };
    inline MATSGN_Complex& operator*=(const fftw_complex&  s)  {
        mWork=mValue[0]*s[0]-mValue[1]*s[1];
        mValue[1]=mValue[0]*s[1]+mValue[1]*s[0];
        mValue[0]=mWork;
        return (*this);
    };
    
    inline MATSGN_Complex& operator/=(const tReal&  s)  {
        if (s==0) throw CORE_Exception("math/signam","MATSGN_Complex/=","division by 0");
        mValue[0]/=s;
        mValue[1]/=s;
        return (*this);
    };
    inline MATSGN_Complex& operator/=(const MATSGN_Complex&  s)  {
        mWork=mValue[0]*s[0]+mValue[1]*s[1];
        mValue[1]=(-mValue[0]*s[1]+mValue[1]*s[0]);
        mValue[0]=mWork;
        mWork=s.module2();
        if (mWork==0) throw CORE_Exception("math/signam","MATSGN_Complex/=","division by 0");
        mValue[0]/=mWork;
        mValue[1]/=mWork;
        return (*this);
    };
    inline MATSGN_Complex& operator+=(const MATSGN_Complex&  s)  {
        mValue[0]+=s[0];
        mValue[1]+=s[1];
        return (*this);
    };
    inline MATSGN_Complex& operator+=(const fftw_complex&  s)  {
        mValue[0]+=s[0];
        mValue[1]+=s[1];
        return (*this);
    };
    inline MATSGN_Complex& operator-=(const MATSGN_Complex&  s)  {
        mValue[0]-=s[0];
        mValue[1]-=s[1];
        return (*this);
    };
    inline MATSGN_Complex& operator-=(const fftw_complex&  s)  {
        mValue[0]-=s[0];
        mValue[1]-=s[1];
        return (*this);
    };
    // SET methods

 


    inline void setValue(const tDouble& x,const tDouble& y) {
        mValue[0]=x;
        mValue[1]=y;
    }
    
    // GET methods

    inline const fftw_complex& getValue() const {
        return mValue;
    }
    inline fftw_complex& getValue()  {
        return mValue;
    }

    inline void conjugate() {
        mValue[1]*=-1;
    }
    inline tDouble module2() const {
        return module2(mValue);
    }
    inline tDouble module() const {
        return module(mValue);
    }
    inline static tDouble module2(const fftw_complex& f)  {
        return f[0]*f[0]+f[1]*f[1];
    }
    inline static tDouble module(const fftw_complex& f) {
        return sqrt(module2(f));
    }
    
    /* \brief  complex product a*=b;
     * @param a: the a value
     * @param b: the b value  
     */
    inline static void product(fftw_complex&  a,const fftw_complex& b) {
        tDouble temp;
        product(a,b,temp);
    }
    /* \brief  complex product a*=b;
     * @param a: the a value
     * @param b: the b value
     */
    inline static void product(fftw_complex&  a,const fftw_complex& b,tDouble& temp) {
        temp=a[0]*b[0]-a[1]*b[1];
        a[1]=a[0]*b[1]+a[1]*b[0];
        a[0]=temp;
    }
    
    /* \brief  complex product r=a*b;
     * @param a: the a value
     * @param b: the b value
     * @param r: the result value
     * @param temp: work value
     */
    inline static void product(const fftw_complex&  a,const fftw_complex& b,fftw_complex& r,tDouble& temp) {
        temp=a[0]*b[0]-a[1]*b[1];
        r[1]=a[0]*b[1]+a[1]*b[0];
        r[0]=temp;
    }
    /* \brief  complex product r=a*b
     * @param a: the a value
     * @param b: the b value
     * @param r: the result value
     */
    inline static void product(const fftw_complex&  a,const fftw_complex& b,fftw_complex& r) {
        ASSERT_IN(&b!=&r);
        ASSERT_IN(&a!=&r);
        r[0]=a[0]*b[0]-a[1]*b[1];
        r[1]=a[0]*b[1]+a[1]*b[0];
    }
    /* \brief  complex product r=a*b;
     * @param a: the a value
     * @param b: the b value
     * @param r: the result value
     */
    inline static void product(const MATSGN_Complex&  a,const MATSGN_Complex& b,MATSGN_Complex& r) {
        const tDouble &ra=a[0];
        const tDouble &ia=a[1];
        const tDouble &rb=b[0];
        const tDouble &ib=b[1];
        tDouble tmp=ra*rb-ia*ib;
        r[1]=ia*rb+ra*ib;
        r[0]=tmp;
    }
    /* \brief  complex product r=(a+ib)*(x+i*y)
     * @param a: real part
     * @param b: imaginary part
     * @param x: real part
     * @param y: imaginary part
     * @param r: result 
     */
    inline static void product(const tDouble&  a,const tDouble& b,const tDouble& x,const tDouble& y,fftw_complex& r) {
        r[0]=a*x-b*y;
        r[1]=a*y+b*x;
    }
    /* \brief  complex add a+=b
     * @param a: the a value
     * @param b: the b value
     */
    inline static void add(fftw_complex&  a,const fftw_complex& b) {
        a[0]+=b[0];
        a[1]+=b[1];
    }
    /* \brief  complex add  r=a+b;
     * @param a: the a value
     * @param b: the b value
     * @param r: the result value
     */
    inline static void add(const fftw_complex&  a,const fftw_complex& b,fftw_complex& r) {
        r[0]=a[0]+b[0];
        r[1]=a[0]+b[1];
    }
    /* \brief  complex sub a-=b;
     * @param a: the a value
     * @param b: the b value
     */
    inline static void sub(fftw_complex&  a,const fftw_complex& b) {
        a[0]-=b[0];
        a[1]-=b[1];
    }
    /* \brief  complex sub : r=a-b;
     * @param a: the a value
     * @param b : the b value
     * @param r: the result value
     */
    inline static void sub(const fftw_complex&  a,const fftw_complex& b,fftw_complex& r) {
        r[0]=a[0]-b[0];
        r[1]=a[1]-b[1];
    }
    /* \brief  swap a & b
     * @param a: the value to swap
     * @param b: the value to swap
     */
    inline static void swap(fftw_complex&  a,fftw_complex& b) {
        std::swap(a[0],b[0]);
        std::swap(a[1],b[1]);
    }

    inline static tDouble norm2(const fftw_complex&  a,const fftw_complex& b) {
       
        return (a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1]);
    }
    inline static tDouble norm(const fftw_complex&  a,const fftw_complex& b) {
       
        return sqrt(norm2(a,b));
    }



    virtual tString toString() const {
        return toString(mValue);
    }
    inline static tString toString(const fftw_complex& f) {
        tString ret="";
        ret+=CORE_Real::toString(f[0]);
        if (f[1]>0) {
            ret+="+"+CORE_Real::toString(f[1])+"i";
        } else if (f[1]<0) {
            ret+=CORE_Real::toString(f[1])+"i";
        }
        return ret;
    }
    inline static tString toString(const fftw_complex& f,const tUSInt& nDigits) {
        tString ret="";
        ret+=CORE_Real::toString(f[0],nDigits);
        if (f[1]>0) {
            ret+="+"+CORE_Real::toString(f[1],nDigits)+"i";
        } else if (f[1]<0) {
            ret+=CORE_Real::toString(f[1],nDigits)+"i";
        }
        return ret;
    }
    
    
    
    
    
    
};

#endif