src/util/ckcomplex.h

   1 #ifndef __CKCOMPLEX_H__
   2 #define __CKCOMPLEX_H__
   3
   4 #include "pup.h"
   5
   6 #if USE_FFTW_DECLS
   7 #include "fftw.h"
   8 typedef fftw_real  RealType;
   9 #else
  10 typedef double     RealType;
  11 #endif
  12
  13 #include <cmath>
  14 #include <ostream>
  15
  16 struct ckcomplex {
  17     RealType  re;
  18     RealType  im;
  19
  20
  21     inline ckcomplex(RealType _re=0., RealType _im=0.): re(_re), im(_im){}
  22     //    inline ckcomplex(RealType r) {re=r; im=0;}
  23     //    inline ckcomplex(RealType r,RealType i) {re=r; im=i;}
  24
  25     inline ~ckcomplex() {}
  26
  27     inline RealType getMagSqr(void) const {
  28       return re*re+im*im;
  29     }
  30
  31     inline ckcomplex operator+(ckcomplex a) {
  32       return ckcomplex(re+a.re,im+a.im);
  33     }
  34
  35     // note: not a member
  36     inline friend ckcomplex operator-(ckcomplex lhs, ckcomplex rhs) {
  37       return ckcomplex(lhs.re - rhs.re, lhs.im - rhs.im);
  38     }
  39
  40     inline ckcomplex conj(void) {
  41         return ckcomplex(re, -im);
  42     }
  43
  44     inline void operator+=(ckcomplex a) {
  45       re+=a.re; im+=a.im;
  46     }
  47
  48     // note: not a member
  49     inline friend ckcomplex operator*(RealType lhs, ckcomplex rhs) {
  50       return ckcomplex(rhs.re*lhs, rhs.im*lhs);
  51     }
  52
  53     inline friend ckcomplex operator/(ckcomplex lhs, RealType rhs) {
  54         return ckcomplex(lhs.re/rhs, lhs.im/rhs);
  55     }
  56
  57     inline ckcomplex operator*(RealType a) {
  58       return ckcomplex(re*a, im*a);
  59     }
  60
  61     inline bool notzero() const { return( (0.0 != re) ? true : (0.0 != im)); }
  62
  63     inline void operator*=(ckcomplex a) {
  64         RealType treal, tim;
  65         treal = re * a.re - im * a.im;
  66         tim = re * a.im + im * a.re;
  67         re = treal;
  68         im = tim;
  69     }
  70
  71     inline ckcomplex operator*(ckcomplex a) {
  72       return ckcomplex( re * a.re - im * a.im, re * a.im + im * a.re);
  73     }
  74
  75
  76     inline void operator -= (ckcomplex a) {
  77       re -= a.re; im -= a.im;
  78     }
  79
  80
  81     inline ckcomplex multiplyByi () {
  82         return ckcomplex(-im, re);
  83     }
  84
  85     inline void * operator new[] (size_t size){
  86         void *buf = malloc(size);
  87         return buf;
  88     }
  89
  90     inline void operator delete[] (void *buf){
  91         free(buf);
  92     }
  93
  94     inline friend std::ostream& operator<< (std::ostream &out, const ckcomplex &aNum) {
  95         out<<"("<<aNum.re<<","<<aNum.im<<")";
  96         return out;
  97     }
  98 };
  99
 100 typedef ckcomplex complex;
 101
 102 PUPbytes(ckcomplex)
 103
 104
 105 /// Overload std::isfinite for complex numbers. @note: Doesn't seem to be part of the standard
 106 inline int isfinite(complex aNum) { return ( std::isfinite(aNum.re) && std::isfinite(aNum.im) ); }
 107
 108 /// Like std::norm, return the square of the distance from (0,0) in the complex number plane
 109 inline RealType norm(complex aNum) { return ( aNum.re*aNum.re + aNum.im*aNum.im ); }
 110 /// Return the distance from (0,0) in the complex plane
 111 inline RealType abs(complex aNum) { return ( sqrt(norm(aNum)) ); }
 112
 113 #endif
 114