Yesterday I needed code to calculate the electric fields of an elliptical two-dimensional Gaussian charge distribution. I decided to implement it in c++ using libcerf for the error function. I used the formula from the Accelerator Physics Handbook.

Beware 1: using Mathematica your results will be noisy because when you create the w-of-z error function there you will intuitively produce something that becomes numerically unstable very easily. Unfortunately AFAIK there is no Mathematica-native version of w-of-z (or the Faddeeva function). Luckily libcerf provides the needed w-of-z so in c/c++ it is not an issue.

Beware 2: I am quite sure there is a 2 missing in the denominator in the fraction in the exponential of the formula for the fields of a round Gaussian beam in the book. Without the 2 (I added it my code below in roundGauss(x,y,a,b)), the results for an ever-so-slightly elliptical beam and a round beam are off by a lot. Also the two is found in other sources and if I remember correctly I re-derived the formula once and arrived at a result with the two in the denominator. TL;DR here is the code:

#include <math.h>
#include <cerf.h>
#include <complex.h>

typedef long double _Complex dcmplx;
dcmplx i = 0 + 1 * I;
#define w(z) w_of_z(z)
#define max(a,b) (a&lt;b)?b:a
#define sab(a, b) sqrt(2. * (pow(a, 2) - pow(b, 2)))
#define rsq(x, y) (pow(x, 2) + pow(y, 2))
#define prefakt(a, b) sqrt(2. * PI / (pow(a, 2) - pow(b, 2)))
#define roundGauss(x, y, a, b) 2. * (i * x + y) / rsq(x, y) * (1. - exp(-rsq(x, y) / (2 * pow(a, 2))))
#define P0(x, y, a, b) prefakt(a, b) * (w((x + i * y) / sab(a, b)) - exp(-(pow(x / a, 2) + pow(y / b, 2)) / 2.) * w((x * b / a + i * y * a / b) / sab(a, b)))
#define P1(x, y, a, b) ((a > b) ? ((y >= 0) ? P0(x, y, a, b) : -1 * conj(P0(x, -y, a, b))) : roundGauss(x, y, a, b))
#define P2(x, y, a, b) (a < b ? i * conj(P1(y, x, b, a)) : P1(x, y, a, b))

The function you want to call is P2(x,y,a,b) and should be equivalent to the function Π given in the Handbook.