/*
* $Id: libcerf.c,v 1.0 2013/07/12 18:11:05 sfeam Exp $
*/
/*
* Complex error function (cerf) and related special functions from libcerf.
* libcerf itself uses the C99 _Complex mechanism for decribing complex
* numbers. This set of wrapper routines converts back and forth from
* gnuplot's own representation of complex numbers.
* Ethan A Merritt - July 2013
*/
#include "syscfg.h"
#ifdef HAVE_LIBCERF
#include <complex.h> /* C99 _Complex */
#include <cerf.h> /* libcerf library header */
#include "eval.h"
#include "stdfn.h" /* for not_a_number */
#include "util.h" /* for int_error() */
#include "libcerf.h" /* our own prototypes */
/* The libcerf complex error function
* cerf(z) = 2/sqrt(pi) * int[0,z] exp(-t^2) dt
*/
void
f_cerf(union argument *arg)
{
struct value a;
complex double z;
pop(&a);
z = real(&a) + I * imag(&a); /* Convert gnuplot complex to C99 complex */
z = cerf(z); /* libcerf complex -> complex function */
push(Gcomplex(&a, creal(z), cimag(z)));
}
/* The libcerf cdawson function returns Dawson's integral
* cdawson(z) = exp(-z^2) int[0,z] exp(t^2) dt
* = sqrt(pi)/2 * exp(-z^2) * erfi(z)
* for complex z.
*/
void
f_cdawson(union argument *arg)
{
struct value a;
complex double z;
pop(&a);
z = real(&a) + I * imag(&a); /* Convert gnuplot complex to C99 complex */
z = cdawson(z); /* libcerf complex -> complex function */
push(Gcomplex(&a, creal(z), cimag(z)));
}
/* The libcerf routine w_of_z returns the Faddeeva rescaled complex error function
* w(z) = exp(-z^2) * erfc(-i*z)
* This corresponds to Abramowitz & Stegun Eqs. 7.1.3 and 7.1.4
* This is also known as the plasma dispersion function.
*/
void
f_faddeeva(union argument *arg)
{
struct value a;
complex double z;
pop(&a);
z = real(&a) + I * imag(&a); /* Convert gnuplot complex to C99 complex */
z = w_of_z(z); /* libcerf complex -> complex function */
push(Gcomplex(&a, creal(z), cimag(z)));
}
/* The libcerf voigt(z, sigma, gamma) function returns the Voigt profile
* corresponding to the convolution
* voigt(x,sigma,gamma) = integral G(t,sigma) L(x-t,gamma) dt
* of Gaussian
* G(x,sigma) = 1/sqrt(2*pi)/|sigma| * exp(-x^2/2/sigma^2)
* with Lorentzian
* L(x,gamma) = |gamma| / pi / ( x^2 + gamma^2 )
* over the integral from -infinity to +infinity.
*/
void
f_voigtp(union argument *arg)
{
struct value a;
double z, sigma, gamma;
gamma = real(pop(&a));
sigma = real(pop(&a));
z = real(pop(&a));
z = voigt(z, sigma, gamma); /* libcerf double -> double function */
push(Gcomplex(&a, z, 0.0));
}
/* The libcery routine re_w_of_z( double x, double y )
* is equivalent to the previously implemented gnuplot routine voigt(x,y).
* Use it in preference to the previous one.
*/
void
f_voigt(union argument *arg)
{
struct value a;
double x, y, w;
y = real(pop(&a));
x = real(pop(&a));
w = re_w_of_z(x, y);
push(Gcomplex(&a, w, 0.0));
}
/* erfi(z) = -i * erf(iz) */
void
f_erfi(union argument *arg)
{
struct value a;
double z;
z = real(pop(&a));
push(Gcomplex(&a, erfi(z), 0.0));
}
#endif