00001 #ifndef __CKCOMPLEX_H__
00002 #define __CKCOMPLEX_H__
00003 
00004 #include "pup.h"
00005 
00006 #if USE_FFTW_DECLS
00007 #include "fftw.h"
00008 typedef fftw_real  RealType;
00009 #else
00010 typedef double     RealType;
00011 #endif
00012 
00013 #include <cmath>
00014 #include <ostream>
00015 
00016 struct ckcomplex {
00017     RealType  re;
00018     RealType  im;   
00019 
00020  
00021     inline ckcomplex(RealType _re=0., RealType _im=0.): re(_re), im(_im){}
00022     
00023     
00024     
00025     inline ~ckcomplex() {}
00026 
00027     inline RealType getMagSqr(void) const { 
00028       return re*re+im*im; 
00029     }
00030 
00031     inline ckcomplex operator+(ckcomplex a) { 
00032       return ckcomplex(re+a.re,im+a.im); 
00033     }
00034 
00035     
00036     inline friend ckcomplex operator-(ckcomplex lhs, ckcomplex rhs) {
00037       return ckcomplex(lhs.re - rhs.re, lhs.im - rhs.im);
00038     }
00039 
00040     inline ckcomplex conj(void) { 
00041         return ckcomplex(re, -im); 
00042     }
00043 
00044     inline void operator+=(ckcomplex a) { 
00045       re+=a.re; im+=a.im; 
00046     }
00047     
00048     
00049     inline friend ckcomplex operator*(RealType lhs, ckcomplex rhs) {
00050       return ckcomplex(rhs.re*lhs, rhs.im*lhs);
00051     }
00052 
00053     inline friend ckcomplex operator/(ckcomplex lhs, RealType rhs) {
00054         return ckcomplex(lhs.re/rhs, lhs.im/rhs);
00055     }
00056 
00057     inline ckcomplex operator*(RealType a) { 
00058       return ckcomplex(re*a, im*a); 
00059     } 
00060 
00061     inline bool notzero() const { return( (0.0 != re) ? true : (0.0 != im)); }
00062 
00063     inline void operator*=(ckcomplex a) {        
00064         RealType treal, tim;
00065         treal = re * a.re - im * a.im;
00066         tim = re * a.im + im * a.re;
00067         re = treal;
00068         im = tim;
00069     }
00070 
00071     inline ckcomplex operator*(ckcomplex a) {
00072       return ckcomplex( re * a.re - im * a.im, re * a.im + im * a.re); 
00073     }
00074 
00075 
00076     inline void operator -= (ckcomplex a) {
00077       re -= a.re; im -= a.im;
00078     }
00079 
00080 
00081     inline ckcomplex multiplyByi () {
00082         return ckcomplex(-im, re);
00083     }
00084         
00085     inline void * operator new[] (size_t size){
00086         void *buf = malloc(size);
00087         return buf;
00088     }
00089     
00090     inline void operator delete[] (void *buf){
00091         free(buf);
00092     }
00093 
00094     inline friend std::ostream& operator<< (std::ostream &out, const ckcomplex &aNum) {
00095         out<<"("<<aNum.re<<","<<aNum.im<<")";
00096         return out;
00097     }
00098 };
00099 
00100 typedef ckcomplex complex;
00101 
00102 PUPbytes(ckcomplex)
00103 
00104 
00105 
00106 inline int isfinite(complex aNum) { return ( std::isfinite(aNum.re) && std::isfinite(aNum.im) ); }
00107 
00109 inline RealType norm(complex aNum) { return ( aNum.re*aNum.re + aNum.im*aNum.im ); }
00111 inline RealType abs(complex aNum) { return ( sqrt(norm(aNum)) ); }
00112 
00113 #endif
00114