Computes the complex error function $ W(z) $ for a given complex number $ z $ . More...

#include "mdb.h"
#include "f2c.h"

Macros
#define	abs(x)

Functions
double	pow_di (doublereal ap, integer bp)

integer	i_dnnt (doublereal *x)

int	wofz (doublereal xi, doublereal yi, doublereal u, doublereal v, logical *flag__)
	Computes the complex error function for a given complex number .

Detailed Description

Computes the complex error function $ W(z) $ for a given complex number $ z $ .

Copyright

(c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
(c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License: This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author: H. Shang, M. Borland, R. Soliday

Definition in file wofz.c.

Macro Definition Documentation

◆ abs

#define abs ( x )

Value:

((x) >= 0 ? (x) : -(x))

Definition at line 46 of file wofz.c.

Function Documentation

◆ i_dnnt()

integer i_dnnt ( doublereal * x )

Definition at line 48 of file wofz.c.

                              {
  return (integer)(*x >= 0. ? floor(*x + .5) : -floor(.5 - *x));
}

◆ pow_di()

double pow_di	(	doublereal *	ap,
		integer *	bp )

Definition at line 19 of file wofz.c.

                                           {
  double pow, x;
  integer n;
  unsigned long u;
 
  pow = 1;
  x = *ap;
  n = *bp;
 
  if (n != 0) {
    if (n < 0) {
      n = -n;
      x = 1 / x;
    }
    for (u = n;;) {
      if (u & 01)
        pow *= x;
      if (u >>= 1)
        x *= x;
      else
        break;
    }
  }
  return (pow);
}

◆ wofz()

int wofz	(	doublereal *	xi,
		doublereal *	yi,
		doublereal *	u,
		doublereal *	v,
		logical *	flag__ )

Computes the complex error function $ W(z) $ for a given complex number $ z $ .

Given a complex number $z = (\text{xi}, \text{yi})$ , this function computes the value of the Faddeeva function $W(z) = \exp(-z^2) \cdot \text{ERFC}(-i z)$ , where $\text{ERFC}$ is the complex complementary error function and $i = \sqrt{-1}$ .

The accuracy of the algorithm is:

14 significant digits for in the 1st and 2nd quadrants.
13 significant digits for in the 3rd and 4th quadrants, outside a circular region with a radius of 0.126 around a zero of the function.

Parameters

xi	Pointer to the real part of the input complex number .
yi	Pointer to the imaginary part of the input complex number .
u	Pointer to store the real part of the computed .
v	Pointer to store the imaginary part of the computed .
flag__	Pointer to a logical flag indicating whether an overflow has occurred. `false:` No error condition. `true:` Overflow has occurred; the routine becomes inactive.

Returns: Returns 0 on successful computation. If an overflow is detected, flag__ is set to true.

Definition at line 79 of file wofz.c.

{
  /* System generated locals */
  doublereal d__1, d__2;
 
  /* Builtin functions */
  /*
    double sqrt();
    integer i_dnnt();
    double exp(), cos(), sin(), pow_di();
    */
  /* Local variables */
  static integer kapn;
  static doublereal xabs, yabs, daux, qrho, xaux, xsum, ysum;
  static logical a, b;
  static doublereal c__, h__;
  static integer i__, j, n;
  static doublereal x, y, xabsq, xquad, yquad, h2, u1, v1, u2, v2, w1;
  static integer nu;
  static doublereal rx, ry, sx, sy, tx, ty;
  static integer np1;
  static doublereal qlambda;
 
  /*  GIVEN A COMPLEX NUMBER Z = (XI,YI), THIS SUBROUTINE COMPUTES */
  /*  THE VALUE OF THE FADDEEVA-FUNCTION W(Z) = EXP(-Z**2)*ERFC(-I*Z), */
  /*  WHERE ERFC IS THE COMPLEX COMPLEMENTARY ERROR-FUNCTION AND I */
  /*  MEANS SQRT(-1). */
  /*  THE ACCURACY OF THE ALGORITHM FOR Z IN THE 1ST AND 2ND QUADRANT */
  /*  IS 14 SIGNIFICANT DIGITS; IN THE 3RD AND 4TH IT IS 13 SIGNIFICANT */
  /*  DIGITS OUTSIDE A CIRCULAR REGION WITH RADIUS 0.126 AROUND A ZERO */
  /*  OF THE FUNCTION. */
  /*  ALL REAL VARIABLES IN THE PROGRAM ARE DOUBLE PRECISION. */
 
  /*  THE CODE CONTAINS A FEW COMPILER-DEPENDENT PARAMETERS : */
  /*     RMAXREAL = THE MAXIMUM VALUE OF RMAXREAL EQUALS THE ROOT OF */
  /*                RMAX = THE LARGEST NUMBER WHICH CAN STILL BE */
  /*                IMPLEMENTED ON THE COMPUTER IN DOUBLE PRECISION */
  /*                FLOATING-POINT ARITHMETIC */
  /*     RMAXEXP  = LN(RMAX) - LN(2) */
  /*     RMAXGONI = THE LARGEST POSSIBLE ARGUMENT OF A DOUBLE PRECISION */
  /*                GONIOMETRIC FUNCTION (DCOS, DSIN, ...) */
  /*  THE REASON WHY THESE PARAMETERS ARE NEEDED AS THEY ARE DEFINED WILL */
  /*  BE EXPLAINED IN THE CODE BY MEANS OF COMMENTS */
 
  /*  PARAMETER LIST */
  /*     XI     = REAL      PART OF Z */
  /*     YI     = IMAGINARY PART OF Z */
  /*     U      = REAL      PART OF W(Z) */
  /*     V      = IMAGINARY PART OF W(Z) */
  /*     FLAG   = AN ERROR FLAG INDICATING WHETHER OVERFLOW WILL */
  /*              OCCUR OR NOT; TYPE LOGICAL; */
  /*              THE VALUES OF THIS VARIABLE HAVE THE FOLLOWING */
  /*              MEANING : */
  /*              FLAG=.FALSE. : NO ERROR CONDITION */
  /*              FLAG=.TRUE.  : OVERFLOW WILL OCCUR, THE ROUTINE */
  /*                             BECOMES INACTIVE */
  /*  XI, YI      ARE THE INPUT-PARAMETERS */
  /*  U, V, FLAG  ARE THE OUTPUT-PARAMETERS */
 
  /*  FURTHERMORE THE PARAMETER FACTOR EQUALS 2/SQRT(PI) */
 
  /*  THE ROUTINE IS NOT UNDERFLOW-PROTECTED BUT ANY VARIABLE CAN BE */
  /*  PUT TO 0 UPON UNDERFLOW; */
 
  /*  REFERENCE - GPM POPPE, CMJ WIJERS; MORE EFFICIENT COMPUTATION OF */
  /*  THE COMPLEX ERROR-FUNCTION, ACM TRANS. MATH. SOFTWARE. */
 
  *flag__ = FALSE_;
 
  xabs = abs(*xi);
  yabs = abs(*yi);
  x = xabs / (float)6.3;
  y = yabs / (float)4.4;
 
  /*     THE FOLLOWING IF-STATEMENT PROTECTS */
  /*     QRHO = (X**2 + Y**2) AGAINST OVERFLOW */
 
  if (xabs > 5e153 || yabs > 5e153) {
    goto L100;
  }
 
  /* Computing 2nd power */
  d__1 = x;
  /* Computing 2nd power */
  d__2 = y;
  qrho = d__1 * d__1 + d__2 * d__2;
 
  /* Computing 2nd power */
  d__1 = xabs;
  xabsq = d__1 * d__1;
  /* Computing 2nd power */
  d__1 = yabs;
  xquad = xabsq - d__1 * d__1;
  yquad = xabs * 2 * yabs;
 
  a = qrho < .085264;
 
  if (a) {
 
    /*  IF (QRHO.LT.0.085264D0) THEN THE FADDEEVA-FUNCTION IS EVALUATED */
    /*  USING A POWER-SERIES (ABRAMOWITZ/STEGUN, EQUATION (7.1.5), P.297) 
*/
    /*  N IS THE MINIMUM NUMBER OF TERMS NEEDED TO OBTAIN THE REQUIRED */
    /*  ACCURACY */
 
    qrho = (1 - y * (float).85) * sqrt(qrho);
    d__1 = qrho * 72 + 6;
    n = i_dnnt(&d__1);
    j = (n << 1) + 1;
    xsum = (float)1. / j;
    ysum = 0.;
    for (i__ = n; i__ >= 1; --i__) {
      j += -2;
      xaux = (xsum * xquad - ysum * yquad) / i__;
      ysum = (xsum * yquad + ysum * xquad) / i__;
      xsum = xaux + (float)1. / j;
      /* L10: */
    }
    u1 = (xsum * yabs + ysum * xabs) * -1.12837916709551257388 + (float)1.;
    v1 = (xsum * xabs - ysum * yabs) * 1.12837916709551257388;
    daux = exp(-xquad);
    u2 = daux * cos(yquad);
    v2 = -daux * sin(yquad);
 
    *u = u1 * u2 - v1 * v2;
    *v = u1 * v2 + v1 * u2;
 
  } else {
 
    /*  IF (QRHO.GT.1.O) THEN W(Z) IS EVALUATED USING THE LAPLACE */
    /*  CONTINUED FRACTION */
    /*  NU IS THE MINIMUM NUMBER OF TERMS NEEDED TO OBTAIN THE REQUIRED */
    /*  ACCURACY */
 
    /*  IF ((QRHO.GT.0.085264D0).AND.(QRHO.LT.1.0)) THEN W(Z) IS EVALUATED
 */
    /*  BY A TRUNCATED TAYLOR EXPANSION, WHERE THE LAPLACE CONTINUED FRACT
ION */
    /*  IS USED TO CALCULATE THE DERIVATIVES OF W(Z) */
    /*  KAPN IS THE MINIMUM NUMBER OF TERMS IN THE TAYLOR EXPANSION NEEDED
 */
    /*  TO OBTAIN THE REQUIRED ACCURACY */
    /*  NU IS THE MINIMUM NUMBER OF TERMS OF THE CONTINUED FRACTION NEEDED
 */
    /*  TO CALCULATE THE DERIVATIVES WITH THE REQUIRED ACCURACY */
 
    if (qrho > (float)1.) {
      h__ = 0.;
      kapn = 0;
      qrho = sqrt(qrho);
      nu = (integer)(1442 / (qrho * 26 + 77) + 3);
    } else {
      qrho = (1 - y) * sqrt(1 - qrho);
      h__ = qrho * (float)1.88;
      h2 = h__ * 2;
      d__1 = qrho * 34 + 7;
      kapn = i_dnnt(&d__1);
      d__1 = qrho * 26 + 16;
      nu = i_dnnt(&d__1);
    }
 
    b = h__ > (float)0.;
 
    if (b) {
      qlambda = pow_di(&h2, &kapn);
    }
 
    rx = (float)0.;
    ry = (float)0.;
    sx = (float)0.;
    sy = (float)0.;
 
    for (n = nu; n >= 0; --n) {
      np1 = n + 1;
      tx = yabs + h__ + np1 * rx;
      ty = xabs - np1 * ry;
      /* Computing 2nd power */
      d__1 = tx;
      /* Computing 2nd power */
      d__2 = ty;
      c__ = (float).5 / (d__1 * d__1 + d__2 * d__2);
      rx = c__ * tx;
      ry = c__ * ty;
      if (b && n <= kapn) {
        tx = qlambda + sx;
        sx = rx * tx - ry * sy;
        sy = ry * tx + rx * sy;
        qlambda /= h2;
      }
      /* L11: */
    }
 
    if (h__ == (float)0.) {
      *u = rx * 1.12837916709551257388;
      *v = ry * 1.12837916709551257388;
    } else {
      *u = sx * 1.12837916709551257388;
      *v = sy * 1.12837916709551257388;
    }
 
    if (yabs == (float)0.) {
      /* Computing 2nd power */
      d__1 = xabs;
      *u = exp(-(d__1 * d__1));
    }
  }
 
  /*  EVALUATION OF W(Z) IN THE OTHER QUADRANTS */
 
  if (*yi < (float)0.) {
 
    if (a) {
      u2 *= 2;
      v2 *= 2;
    } else {
      xquad = -xquad;
 
      /*         THE FOLLOWING IF-STATEMENT PROTECTS 2*EXP(-Z**2) */
      /*         AGAINST OVERFLOW */
 
      if (yquad > 3537118876014220. || xquad > 708.503061461606) {
        goto L100;
      }
 
      w1 = exp(xquad) * 2;
      u2 = w1 * cos(yquad);
      v2 = -w1 * sin(yquad);
    }
 
    *u = u2 - *u;
    *v = v2 - *v;
    if (*xi > (float)0.) {
      *v = -(*v);
    }
  } else {
    if (*xi < (float)0.) {
      *v = -(*v);
    }
  }
 
  return 0;
 
L100:
  *flag__ = TRUE_;
  return 0;
 
} /* wofz */

Macros