ComplexErrorFun.cpp

Go to the documentation of this file.

// ------------------------------------------------------------------------
// $RCSfile: ComplexErrorFun.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.2 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Function: Werrf(complex<double>)
//   The complex error function.
//
// ------------------------------------------------------------------------
// Class category: Utilities
// ------------------------------------------------------------------------
//
// $Date: 2001/08/24 19:33:44 $
// $Author: jsberg $
//
// ------------------------------------------------------------------------
 
#include "Utilities/ComplexErrorFun.h"
#include <complex>
 
// Complex error function
// The algorithms is based on:
//   Walter Gautschi,
//   Efficient Computation of the Complex Error Function,
//   SIAM J. Numer. Anal., Vol 7, No. 1, March 1970, pp. 187-198.
// ------------------------------------------------------------------------
 
std::complex<double> Werrf(std::complex<double> z) {
    const double cc   = 1.12837916709551;  // 2 / sqrt(pi)
    const double xlim = 5.33;
    const double ylim = 4.29;
 
    const double x = std::abs(std::real(z));
    const double y = std::abs(std::imag(z));
    std::complex<double> w;
 
    if (y < ylim && x < xlim) {
        // Inside limit rectangle: equation (3.8): q = s(z);
        const double q = (1.0 - y / ylim) * sqrt(1.0 - (x * x) / (xlim * xlim));
        // Equations (3.11).
        const int nc = 6 + int(23.0 * q);
        const int nu = 9 + int(21.0 * q);
        const std::complex<double> zz(1.6 * q + y, -x);  // h - i*z
 
        // Equations (3.12) for n = nu, nu - 1, ..., N.
        std::complex<double> r(0.0);             // r_{nu}
        for (int n = nu; n > nc; n--) {          // for n = nu, nu-1, ..., nc+1
            r = 0.5 / (zz + double(n + 1) * r);  // r_{n-1}
        }
 
        // Equations (3.12) for n = N, N - 1, ..., 0.
        double power     = pow(3.2 * q, nc);     // (2*h)^nc
        const double div = 0.3125 / q;           // 1 / (2*h)
        std::complex<double> s(0.0);             // s_{N}
        for (int n = nc; n >= 0; n--) {          // for n = N, N-1, ..., 0
            r = 0.5 / (zz + double(n + 1) * r);  // r_{n-1}
            s = r * (power + s);                 // s_{n-1}
            power *= div;                        // (2*h)^(n-1)
        }
 
        // Equation (3.13).
        w = cc * s;
    } else {
        // Outside limit rectangle: equations (3.12) for h = 0.
        const std::complex<double> zz(y, -x);    // - i*z
        std::complex<double> r(0.0);             // r_{nu}
        for (int n = 8; n >= 0; n--) {           // for n = nu, nu-1, ..., 0
            r = 0.5 / (zz + double(n + 1) * r);  // r_{n-1}
        }
 
        // Equation (3.13).
        w = cc * r;
    }
 
    // Equations (3.1).
    if (std::imag(z) < 0.0) {
        std::complex<double> zz(x, y);
        w = 2.0 * std::exp(-zz * zz) - w;
        if (std::real(z) > 0.0) w = std::conj(w);
    } else {
        if (std::real(z) < 0.0) w = std::conj(w);
    }
 
    return w;
}