Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

erfz(zz)
```function f = erfz(zz)
%ERFZ  Error function for complex inputs
%   f = erfz(z) is the error function for the elements of z.
%   Z may be complex and of any size.
%   Accuracy is better than 12 significant digits.
%
%   Usage:  f = erfz(z)
%
%   Ref: Abramowitz & Stegun section 7.1
%   equations 7.1.9, 7.1.23, and 7.1.29
%
%   Tested under version 5.3.1
%

%   Main author Paul Godfrey <pgodfrey@intersil.com>
%   Small changes by Peter J. Acklam <jacklam@math.uio.no>
%   09-26-01

error(nargchk(1, 1, nargin));

% quick exit for empty input
if isempty(zz)
f = zz;
return;
end

twopi = 2*pi;
sqrtpi=sqrt(pi);

f = zeros(size(zz));
ff=f;

az=abs(zz);
p1=find(az<=8);
p2=find(az> 8);

if ~isempty(p1)
z=zz(p1);

nn = 32;

x = real(z);
y = imag(z);
k1 = 2 / pi * exp(-x.*x);
k2 = exp(-i*2*x.*y);

s1 = erf(x);

s2 = zeros(size(x));
k = x ~= 0;          % when x is non-zero
s2(k) = k1(k) ./ (4*x(k)) .* (1 - k2(k));
k = ~k;              % when x is zero
s2(k) = i / pi * y(k);

f = s1 + s2;

k = y ~= 0;          % when y is non-zero
xk = x(k);
yk = y(k);

s5 = 0;
for n = 1 : nn
s3 = exp(-n*n/4) ./ (n*n + 4*xk.*xk);
s4 = 2*xk - k2(k).*(2*xk.*cosh(n*yk) - i*n*sinh(n*yk));
s5 = s5 + s3.*s4;
end
s6 = k1(k) .* s5;
f(k) = f(k) + s6;
ff(p1)=f;
end

if ~isempty(p2)
z=zz(p2);
pn=find(real(z)<0);

if ~isempty(pn)
z(pn)=-z(pn);
end

nmax=193;
s=1;
y=2*z.*z;
for n=nmax:-2:1
s=1-n.*(s./y);
end

f=1.0-s.*exp(-z.*z)./(sqrtpi*z);

if ~isempty(pn)
f(pn)=-f(pn);
end

pa=find(real(z)==0);
%  fix along i axis problem
if ~isempty(pa)
f(pa)=f(pa)-1;
end

ff(p2)=f;
end

f=ff;

return

%a demo of this function is
x = -4:0.125:4;
y = x;
[X, Y] = meshgrid(x,y);
z = complex(X, Y);
f = erfz(z);
af = abs(f);
%let's truncate for visibility
p = find(af > 5);
af(p) = 5;
mesh(x, y, af);
view(-70, 40);
rotate3d on;

return
```