core/complex.h

Go to the documentation of this file.

#ifndef _COMPLEX_H
#define _COMPLEX_H

typedef enum {I='i',J='j',A='d'} CNOTATION; 
#define CNOTATION_DEFAULT J /* never set this to A */
#define PI 3.1415926535897932384626433832795
#define E 2.71828182845905

/* only cpp code may actually do complex math */
#ifdef __cplusplus
#include <math.h>
#include "platform.h"
class complex { 
private:
    double r; 
    double i; 
    CNOTATION f; 
public:
    inline complex() 
    {
        r = 0;
        i = 0;
        f = CNOTATION_DEFAULT;
    };
    inline complex(double re) 
    {
        r = re;
        i = 0;
        f = CNOTATION_DEFAULT;
    };
    inline complex(double re, double im, CNOTATION nf=CNOTATION_DEFAULT) 
    {
        r = re;
        i = im;
        f = nf;
    };
    
    /* assignment operations */
    inline complex &operator = (complex x) 
    {
        r = x.r; 
        i = x.i; 
        f = x.f; 
        return *this;
    };
    inline complex &operator = (double x) 
    {
        r = x; 
        i = 0; 
        f = CNOTATION_DEFAULT;
        return *this;
    };

    /* access operations */
    inline double & Re(void) 
    {
        return r;
    };
    inline double & Im(void) 
    {
        return i;
    };
    inline CNOTATION & Notation(void) 
    {
        return f;
    };
    inline double Mag(void) const 
    {
        return sqrt(r*r+i*i);
    };
    inline double Mag(double m)  
    {
        double old = sqrt(r*r+i*i);
        r *= m/old;
        i *= m/old;
        return m;
    };
    inline double Arg(void) const 
    {
        if (r==0)
        {
            if (i>0)
                return PI/2;
            else if (i==0)
                return 0;
            else
                return -PI/2;
        }
        else if (r>0)
            return atan(i/r);
        else
            return PI+atan(i/r);
    };
    inline double Arg(double a)  
    {
        SetPolar(Mag(),a,f);
        return a;
    };
    inline complex Log(void) const 
    { 
        return complex(log(Mag()),Arg(),f);
    };
    inline void SetReal(double v) 
    {
        r = v;
    };
    inline void SetImag(double v) 
    {
        i = v;
    };
    inline void SetNotation(CNOTATION nf) 
    {
        f = nf;
    }
    inline void SetRect(double rp, double ip, CNOTATION nf=CNOTATION_DEFAULT) 
    {
        r = rp;
        i = ip;
        f = nf;
    };
    inline void SetPolar(double m, double a, CNOTATION nf=A) 
    {
        r = (m*cos(a)); 
        i = (m*sin(a));
        f = nf;
    };

#if 0
    //inline operator const double (void) const /**< cast real part to double */
    //{
    //  return r;
    //};
#endif

    inline complex operator - (void) 
    {
        return complex(-r,-i,f);
    };
    inline complex operator ~ (void) 
    { 
        return complex(r,-i,f);
    };

    /* reflexive math operations */
    inline complex &operator += (double x) 
    {
        r += x; 
        return *this;
    };
    inline complex &operator -= (double x) 
    {
        r -= x; 
        return *this;
    };
    inline complex &operator *= (double x) 
    {
        r *= x; 
        i *= x; 
        return *this;
    };
    inline complex &operator /= (double x) 
    {
        r /= x; 
        i /= x;
        return *this;
    };
    inline complex &operator ^= (double x) 
    { 
        double lm = log(Mag()), a = Arg(), b = exp(x*lm), c = x*a; 
        r = (b*cos(c)); 
        i = (b*sin(c)); 
        return *this;
    };
    inline complex &operator += (complex x) 
    {
        r += x.r; 
        i += x.i; 
        return *this;
    };
    inline complex &operator -= (complex x)  
    {
        r -= x.r; 
        i -= x.i; 
        return *this;
    };
    inline complex &operator *= (complex x)  
    {
        double pr=r; 
        r = pr * x.r - i * x.i; 
        i = pr * x.i + i * x.r; 
        return *this;
    };
    inline complex &operator /= (complex y)  
    {
        double xr=r;
        double a = y.r*y.r+y.i*y.i;
        r = (xr*y.r+i*y.i)/a; 
        i = (i*y.r-xr*y.i)/a;
        return *this;
    };
    inline complex &operator ^= (complex x) 
    { 
        double lm = log(Mag()), a = Arg(), b = exp(x.r*lm-x.i*a), c = x.r*a+x.i*lm; 
        r = (b*cos(c)); 
        i = (b*sin(c)); 
        return *this;
    };

    /* binary math operations */
    inline complex operator + (double y) 
    {
        complex x(*this);
        return x+=y;
    };
    inline complex operator - (double y) 
    {
        complex x(*this);
        return x-=y;
    };
    inline complex operator * (double y) 
    {
        complex x(*this);
        return x*=y;
    };
    inline complex operator / (double y) 
    {
        complex x(*this);
        return x/=y;
    };
    inline complex operator ^ (double y) 
    {
        complex x(*this);
        return x^=y;
    };
    inline complex operator + (complex y) 
    {
        complex x(*this);
        return x+=y;
    };
    inline complex operator - (complex y) 
    {
        complex x(*this);
        return x-=y;
    };
    inline complex operator * (complex y) 
    {
        complex x(*this);
        return x*=y;
    };
    inline complex operator / (complex y) 
    {
        complex x(*this);
        return x/=y;
    };
    inline complex operator ^ (complex y) 
    {
        complex x(*this);
        return x^=y;
    };

    inline complex SetPowerFactor(double mag, 
        double pf, 
        CNOTATION n=J) 
    {
        SetPolar(mag/pf, acos(pf),n);
        return *this;
    }


    /* comparison */
    inline bool IsZero(double err=0.0) 
    {
        return Mag()<=err;
    }

    /* magnitude comparisons */
    inline bool operator == (double m)  { return Mag()==m; };
    inline bool operator != (double m)  { return Mag()!=m; };
    inline bool operator < (double m)   { return Mag()<m; };
    inline bool operator <= (double m)  { return Mag()<=m; };
    inline bool operator > (double m)   { return Mag()>m; };
    inline bool operator >= (double m)  { return Mag()>=m; };

    /* angle comparisons */
    inline complex operator == (complex y)  { return fmod(y.Arg()-Arg(),2*PI)==0.0;};
    inline complex operator != (complex y)  { return fmod(y.Arg()-Arg(),2*PI)!=0.0;};
    inline complex operator < (complex y)   { return fmod(y.Arg()-Arg(),2*PI)<PI;};
    inline complex operator <= (complex y)  { return fmod(y.Arg()-Arg(),2*PI)<=PI;};
    inline complex operator > (complex y)   { return fmod(y.Arg()-Arg(),2*PI)>PI;};
    inline complex operator >= (complex y)  { return fmod(y.Arg()-Arg(),2*PI)>=PI;};
#ifndef isfinite
#define isfinite(X) (!isinf(X) && !isnan(X))
#endif
    inline bool IsFinite(void) { return isfinite(r) && isfinite(i); };
};
#else
typedef struct { 
    double r; 
    double i; 
    CNOTATION f; 
} complex; 
#endif

#endif

---

GridLAB-D^TM Version 1.0

An open-source project initiated by the US Department of Energy