#ifndef COMPLEX_H #define COMPLEX_H /* * METRIC --- Mode expansion modeling in integrated optics / photonics * http://metric.computational-photonics.eu/ */ /* * complex.h * Handling complex numbers */ #include #include #include #include #include "inout.h" /* Complex: complex number */ class Complex { public: // real & imaginary part double re; double im; // initialize inline Complex() { re=0.0; im=0.0; } inline Complex(const double& x) { re=x; im=0.0; } inline Complex(const double& x, const double& y) { re=x; im=y; } // copy & assignment inline Complex(const Complex& z) { re = z.re; im = z.im; } inline Complex& operator=(const Complex& z) { if(this != &z) { re = z.re; im = z.im; } return *this; } // the negative inline Complex operator-() const { return Complex(-re, -im); } // complex conjugation inline Complex conj() const { return Complex(re, -im); } // modulus inline double abs() const { return sqrt(re*re+im*im); } // modulus squared inline double sqabs() const { return re*re+im*im; } // argument; double arg() const; // range [0, 2 PI[ double args() const; // range ]-PI, PI] // multiply by (-1) & assignment inline void neg() { re = -re; im = -im; } // multiply by (i) & assignment inline void muli() { double t=re; re = -im; im = t; } // multiply by (-i) & assignment inline void mulmi() { double t=re; re = im; im = -t; } // 1/(z) & assignment inline void inv() { double t = re*re+im*im; re = re/t; im = -im/t; } // addition inline Complex operator+(const Complex& z) const { return Complex(re+z.re, im+z.im); } inline Complex operator+(const double& x) const { return Complex(re+x, im); } friend inline Complex operator+(const double& x, const Complex& z) { return Complex(x+z.re, z.im); } // addition & assignment inline void operator+=(const Complex& z) { re+=z.re; im+=z.im; return;} inline void operator+=(const double& x) { re+=x; return; } // subtraction inline Complex operator-(const Complex& z) const { return Complex(re-z.re, im-z.im); } inline Complex operator-(const double& x) const { return Complex(re-x, im); } friend inline Complex operator-(const double& x, const Complex& z) { return Complex(x-z.re, -z.im); } // subtraction & assignment inline void operator-=(const Complex& z) { re-=z.re; im-=z.im; return; } inline void operator-=(const double& x) { re-=x; return; } // multiplication inline Complex operator*(const Complex& z) const { return Complex(re*z.re-im*z.im, re*z.im+im*z.re); } inline Complex operator*(const double& x) const { return Complex(re*x, im*x); } friend inline Complex operator*(const double& x, const Complex& z) { return Complex(x*z.re, x*z.im); } // multiplication & assignment inline void operator*=(const Complex& z) { double r = re, i = im; re = r*z.re-i*z.im; im = r*z.im+i*z.re; return; } inline void operator*=(const double& x) { re*=x; im*=x; return; } // division inline Complex operator/(const Complex& z) const { double den=z.re*z.re+z.im*z.im; return Complex((re*z.re+im*z.im)/den, (im*z.re-re*z.im)/den); } inline Complex operator/(const double& x) const { return Complex(re/x, im/x); } friend inline Complex operator/(const double& x, const Complex& z) { double den=z.re*z.re+z.im*z.im; return Complex(x*z.re/den, -x*z.im/den); } // division & assignment inline void operator/=(const Complex& z) { double r = re, i = im; double den=z.re*z.re+z.im*z.im; re = (r*z.re+i*z.im)/den; im = (i*z.re-r*z.im)/den; return; } inline void operator/=(const double& x) { re/=x; im/=x; return; } // complex conjugate & multiply inline Complex conjmult(const Complex& z) { return Complex(re*z.re+im*z.im, re*z.im-im*z.re); } inline Complex conjmult(const double& x) { return Complex(re*x, -im*x); } // input void read(FILE *dat); // output void write(FILE *dat) const; void nice(FILE *dat) const; // ckeck validity: re, im in {nan, -nan, inf, -inf} is not OK bool isOK() const; bool isnotOK() const; }; /* some important constants ... */ #define CC0 (Complex(0.0, 0.0)) #define CC1 (Complex(1.0, 0.0)) #define CCm1 (Complex(-1.0, 0.0)) #define CCi (Complex(0.0, 1.0)) #define CCI (Complex(0.0, 1.0)) #define CCmi (Complex(0.0, -1.0)) #define CCmI (Complex(0.0, -1.0)) /* exponential function */ /* exp(c.re)*(cos(c.im) + i sin(c.im)) */ inline Complex exp(const Complex& c) { double e = exp(c.re); return Complex(e*cos(c.im), e*sin(c.im)); } /* cos(phi) + i sin(phi) */ inline Complex expi(double phi) { return Complex(cos(phi), sin(phi)); } /* exp(-c.im)*(cos(c.re) + i sin(c.re)) */ inline Complex expi(const Complex& c) { double e = exp(-c.im); return Complex(e*cos(c.re), e*sin(c.re)); } /* exp(c.im)*(cos(c.re) - i sin(c.re)) */ inline Complex expmi(const Complex& c) { double e = exp(c.im); return Complex(e*cos(c.re), -e*sin(c.re)); } /* trigonometric functions */ Complex sin(const Complex& z); Complex cos(const Complex& z); /* squareroot s = 0: sqrt(|c|)*exp(i arg(c)/2), positive squareroot for real c s = 1: sqrt(|c|)*exp(i arg(c)/2 + pi), negative squareroot for real c */ Complex sqrt(const Complex& c, int s); Complex sqrt_pos(const Complex& c); Complex sqrt_neg(const Complex& c); /* principal value; +/-sqrt(|c|)*exp(i arg(c)/2), sign such that real part of the result is positive; sign of pure imaginary result z.im equal to the sign of c.im */ Complex sqrt(const Complex& c); /* natural logarithm, imaginary part in ]-PI, PI] */ Complex log(const Complex& z); /* powers of complex numbers */ Complex pow(const Complex& z, const int& n); Complex pow(const double& d, const Complex& c); Complex pow(const Complex&z, const double& d); Complex pow(const Complex& z, const Complex& c); #endif // COMPLEX_H