complexI.H

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2014 OpenFOAM Foundation
    Copyright (C) 2019-2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.
 
    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
inline constexpr Foam::complex::complex() noexcept
:
    re(0),
    im(0)
{}
 
 
inline constexpr Foam::complex::complex(Foam::zero) noexcept
:
    re(0),
    im(0)
{}
 
 
inline constexpr Foam::complex::complex(const scalar r) noexcept
:
    re(r),
    im(0)
{}
 
 
inline constexpr Foam::complex::complex(const scalar r, const scalar i) noexcept
:
    re(r),
    im(i)
{}
 
 
inline Foam::complex::complex(const std::complex<float>& c)
:
    re(c.real()),
    im(c.imag())
{}
 
 
inline Foam::complex::complex(const std::complex<double>& c)
:
    re(c.real()),
    im(c.imag())
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
inline Foam::complex Foam::complex::conjugate() const
{
    return complex(re, -im);
}
 
 
inline Foam::scalar Foam::complex::magnitude() const
{
    return std::hypot(re, im);
}
 
 
inline Foam::scalar Foam::complex::magSqr() const
{
    return (re*re + im*im);
}
 
 
inline Foam::scalar Foam::complex::cmptSum() const noexcept
{
    return (re + im);
}
 
 
// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
inline void Foam::complex::operator=(Foam::zero)
{
    re = 0;
    im = 0;
}
 
 
inline void Foam::complex::operator=(const scalar s)
{
    re = s;
    im = 0;
}
 
 
inline void Foam::complex::operator+=(const complex& c)
{
    re += c.re;
    im += c.im;
}
 
 
inline void Foam::complex::operator+=(const scalar s)
{
    re += s;
}
 
 
inline void Foam::complex::operator-=(const complex& c)
{
    re -= c.re;
    im -= c.im;
}
 
 
inline void Foam::complex::operator-=(const scalar s)
{
    re -= s;
}
 
 
inline void Foam::complex::operator*=(const complex& c)
{
    *this = (*this)*c;
}
 
 
inline void Foam::complex::operator*=(const scalar s)
{
    re *= s;
    im *= s;
}
 
 
inline void Foam::complex::operator/=(const complex& c)
{
    *this = *this/c;
}
 
 
inline void Foam::complex::operator/=(const scalar s)
{
    re /= s;
    im /= s;
}
 
 
inline bool Foam::complex::operator==(const complex& c) const
{
    return (equal(re, c.re) && equal(im, c.im));
}
 
 
inline bool Foam::complex::operator!=(const complex& c) const
{
    return !operator==(c);
}
 
 
// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
inline Foam::complex Foam::operator~(const complex& c)
{
    return c.conjugate();
}
 
 
// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
namespace Foam
{
 
inline scalar mag(const complex& c)
{
    return c.magnitude();
}
 
 
inline scalar magSqr(const complex& c)
{
    return c.magSqr();
}
 
 
inline complex sqr(const complex& c)
{
    return c*c;
}
 
 
//- sgn() https://en.wikipedia.org/wiki/Sign_function#Complex_signum
inline complex sign(const complex& c)
{
    // TBD: Use volatile to avoid aggressive branch optimization
    const scalar s(mag(c));
    return s < ROOTVSMALL ? Zero : c/s;
}
 
 
//- csgn() https://en.wikipedia.org/wiki/Sign_function#Complex_signum
inline scalar csign(const complex& c)
{
    return equal(c.real(), 0) ? sign(c.imag()) : sign(c.real());
}
 
 
inline const complex& min(const complex& c1, const complex& c2)
{
    return (c1.magSqr() < c2.magSqr()) ? c1 : c2;
}
 
 
inline const complex& max(const complex& c1, const complex& c2)
{
    return (c1.magSqr() < c2.magSqr()) ? c2 : c1;
}
 
 
inline complex limit(const complex& c1, const complex& c2)
{
    return complex
    (
        limit(c1.real(), c2.real()),
        limit(c1.imag(), c2.imag())
    );
}
 
 
//- Linear interpolation of complex a and b by factor t
inline complex lerp(const complex& a, const complex& b, const scalar t)
{
    const scalar onet = (1-t);
 
    return complex
    (
        onet*a.real() + t*b.real(),
        onet*a.imag() + t*b.imag()
    );
}
 
 
inline const complex& sum(const complex& c)
{
    return c;
}
 
 
template<class Cmpt> class Tensor;
 
//- No-op rotational transform for complex
inline complex transform(const Tensor<scalar>&, const complex c)
{
    return c;
}
 
 
// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
inline complex operator-(const complex& c)
{
    return complex(-c.re, -c.im);
}
 
 
inline complex operator+(const complex& c1, const complex& c2)
{
    return complex
    (
        c1.re + c2.re,
        c1.im + c2.im
    );
}
 
 
inline complex operator+(const complex& c, const scalar s)
{
    return complex(c.re + s, c.im);
}
 
 
inline complex operator+(const scalar s, const complex& c)
{
    return complex(c.re + s, c.im);
}
 
 
inline complex operator-(const complex& c1, const complex& c2)
{
    return complex
    (
        c1.re - c2.re,
        c1.im - c2.im
    );
}
 
 
inline complex operator-(const complex& c, const scalar s)
{
    return complex(c.re - s, c.im);
}
 
 
inline complex operator-(const scalar s, const complex& c)
{
    return complex(s - c.re, -c.im);
}
 
 
inline complex operator*(const complex& c1, const complex& c2)
{
    return complex
    (
        c1.re*c2.re - c1.im*c2.im,
        c1.im*c2.re + c1.re*c2.im
    );
}
 
 
inline complex operator*(const complex& c, const scalar s)
{
    return complex(s*c.re, s*c.im);
}
 
 
inline complex operator*(const scalar s, const complex& c)
{
    return complex(s*c.re, s*c.im);
}
 
 
inline complex operator/(const complex& c1, const complex& c2)
{
    const scalar sqrC2 = c2.magSqr();
 
    return complex
    (
        (c1.re*c2.re + c1.im*c2.im)/sqrC2,
        (c1.im*c2.re - c1.re*c2.im)/sqrC2
    );
}
 
 
inline complex operator/(const complex& c, const scalar s)
{
    return complex(c.re/s, c.im/s);
}
 
 
inline complex operator/(const scalar s, const complex& c)
{
    return complex(s)/c;
}
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
} // End namespace Foam
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
// Complex transcendental functions
namespace Foam
{
    #define transFunc(func)                                                   \
    inline complex func(const complex& c)                                     \
    {                                                                         \
        return std:: func (static_cast<std::complex<scalar>>(c));             \
    }
 
transFunc(sqrt)
transFunc(exp)
transFunc(log)
transFunc(log10)
transFunc(sin)
transFunc(cos)
transFunc(tan)
transFunc(asin)
transFunc(acos)
transFunc(atan)
transFunc(sinh)
transFunc(cosh)
transFunc(tanh)
transFunc(asinh)
transFunc(acosh)
transFunc(atanh)
 
// Special treatment for pow()
inline complex pow(const complex& x, const complex& y)
{
    return std::pow
    (
        static_cast<std::complex<scalar>>(x),
        static_cast<std::complex<scalar>>(y)
    );
}
 
 
// Combinations of complex and integral/float
#define powFuncs(type2)                                                       \
    inline complex pow(const complex& x, const type2 y)                       \
    {                                                                         \
        return std::pow(static_cast<std::complex<scalar>>(x), y);             \
    }                                                                         \
                                                                              \
    inline complex pow(const type2 x, const complex& y)                       \
    {                                                                         \
        return std::pow                                                       \
        (                                                                     \
            static_cast<std::complex<scalar>>(x),                             \
            static_cast<std::complex<scalar>>(y)                              \
        );                                                                    \
    }
 
powFuncs(int)
powFuncs(long)
#if (WM_LABEL_SIZE == 64) && (defined(__APPLE__) || defined(_WIN32))
powFuncs(int64_t)
#endif
powFuncs(float)
powFuncs(double)
 
 
inline complex pow3(const complex& c)
{
    return c*sqr(c);
}
 
 
inline complex pow4(const complex& c)
{
    return sqr(sqr(c));
}
 
 
inline complex pow5(const complex& c)
{
    return c*pow4(c);
}
 
 
inline complex pow6(const complex& c)
{
    return pow3(sqr(c));
}
 
 
inline complex pow025(const complex& c)
{
    return sqrt(sqrt(c));
}
 
} // End namespace Foam
 
#undef transFunc
#undef powFuncs
 
// ************************************************************************* //