Code covered by the BSD License

# Generation of Random Variates

### James Huntley (view profile)

generates random variates from over 870 univariate distributions

complexerrorfunc(Z)
```function [W,failure] =  complexerrorfunc(Z)
%
% This function computes the complex scaled complementary error function
% also known as the FADDEEVA function
% W(Z) =  exp(-Z^2) * (1-erf(-i*Z)).
% Note that W(i*X) = =  erfcx(X) for real-valued X. erfcx is a =
% build-in MatLab function.
%
% Inputs and outputs are complex matrices of any size.
% The failure flag indicates numerical overflow.
%
% Source: http://www.netlib.org/toms/680.
% Work published in TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOL. 16, NO. =
% 1, PP. 47.
%
% Translated from F77 to MatLab by Thomas Winiecki 10/05/2006.

% some constants to prevent numerical overflows
RMAX     =  1e308;
RMAXREAL =  sqrt(RMAX);
RMAXEXP  =  log(RMAX) - log(2);
RMAXGONI =  1e10*pi;
FACTOR   =  2/sqrt(pi);

% split input into real/imag
XI =  real(Z);
YI =  imag(Z);

% reset outputs
sizein  =  size(Z);
failure =  repmat(0  ,sizein);
W       =  repmat(NaN,sizein);

for idx =  1:prod(sizein)

XABS =  abs(XI(idx));
YABS =  abs(YI(idx));
X    =  XABS/6.3;
Y    =  YABS/4.4;

if ((XABS > RMAXREAL)|(YABS > RMAXREAL))
failure(idx) =  1;
break;
else
QRHO  =  X^2 + Y^2;
XABSQ =  XABS^2;
A     =  QRHO < 0.085264;
if A
QRHO  =  (1-0.85*Y)*sqrt(QRHO);
N     =  round(6 + 72*QRHO);
J     =  2*N+1;
XSUM  =  1.0/J;
YSUM  =  0.0;
for I =  N:-1:1
J    =  J - 2;
XSUM =  XAUX + 1.0/J;
end
U1   =  -FACTOR*(XSUM*YABS + YSUM*XABS) + 1.0;
V1   =   FACTOR*(XSUM*XABS - YSUM*YABS);
U    =  U1*U2 - V1*V2;
V    =  U1*V2 + V1*U2;
else
if (QRHO > 1.0)
H    =  0.0;
KAPN =  0;
QRHO =  sqrt(QRHO);
NU   =  floor(3 + (1442/(26*QRHO+77)));
else
QRHO =  (1-Y)*sqrt(1-QRHO);
H    =  1.88*QRHO;
H2   =  2*H;
KAPN =  round(7  + 34*QRHO);
NU   =  round(16 + 26*QRHO);
end
B =  H > 0;
if B
QLAMBDA =  H2^KAPN;
end
RX =  0.0;
RY =  0.0;
SX =  0.0;
SY =  0.0;
for N =  NU:-1:0
NP1 =  N + 1;
TX  =  YABS + H + NP1*RX;
TY  =  XABS - NP1*RY;
C   =  0.5/(TX^2 + TY^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 =  QLAMBDA/H2;
end
end
if (H ==  0.0)
U =  FACTOR*RX;
V =  FACTOR*RY;
else
U =  FACTOR*SX;
V =  FACTOR*SY;
end
if (YABS ==  0)
U =  exp(-XABS^2);
end
end
if (YI(idx) < 0.0)
if A
U2    =  2*U2;
V2    =  2*V2;
else
failure(idx) =  1;
break
end
end
U =  U2 - U;
V =  V2 - V;
if XI(idx) > 0
V =  -V;
end
else
if XI(idx) < 0
V =  -V;
end
end
end

W(idx) =  U + i*V;

end

return

```