function/d fVoigtPoppeWijers(x,GamL,GamG)
	variable/d x,GamL,GamG;
//	GamL: FWHM of Lorentzian component
//	GamG: FWHM of Gaussian component
	variable/d gL,gG;
	gL=0.5*GamL;
	gG=0.6005612043932249*GamG;
		// 1/(2*sqrt(ln(2))) = 0.6005612043932249
	if (gG==0.0)
		return 1/(PI*gL*(1+(x/gL)^2));
	endif;
	return 0.5641895835477562 * real(wofz(cmplx(abs(x)/gG,gL/gG)))/gG;
		// 1/sqrt(pi)=0.5641895835477562
end;

function/c wofz(z)
	variable/c z;
//  GIVEN A COMPLEX NUMBER Z = (XI,YI), THIS SUBROUTINE COMPUTES
//  THE VALUE OF THE FADDEEVA-FUNCTION W(Z) = EXP(-Z*Z)*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 (COS, SIN, ...)
//  THE REASON WHY THESE PARAMETERS ARE NEEDED AS THEY ARE DEFINED WILL
//  BE EXPLAINED IN THE CODE BY MEANS OF COMMENTS
//
//
	variable/d FACTOR, RMAXREAL, RMAXEXP, RMAXGONI;
	variable FLAG;
	variable/d XI,YI,XABS,YABS,X,Y;
	variable/d QRHO,XABSQ,XQUAD,YQUAD;
	variable N,J,I
	variable/d XSUM,YSUM;
	variable/d U1,V1,DAUX,U2,V2,U,V,C;
	variable/d H,H2,RX,RY,SX,SY,XAUX,QLAMBDA,TX,TY,W1;
	variable KAPN,NU,NP1;
//
	FACTOR=1.12837916709551257388;
	RMAXREAL = 0.5e154;
	RMAXEXP = 708.503061461606;
	RMAXGONI = 3.53711887601422e15;
	
      FLAG = 0

	XI = real(z);
	YI = imag(z);
      XABS = ABS(XI)
      YABS = ABS(YI)
      X    = XABS/6.3
      Y    = YABS/4.4

//
//     THE FOLLOWING if-STATEMENT PROTECTS
//     QRHO = (X*X + Y*Y) AGAINST OVERFLOW
//
	if ((XABS > RMAXREAL) %| (YABS > RMAXREAL))
		return cmplx(inf,inf);
	endif;

      QRHO = X*X + Y*Y

      XABSQ = XABS*XABS
      XQUAD = XABSQ - YABS*YABS
      YQUAD = 2*XABS*YABS

	if (QRHO < 0.085264)
//
//  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-0.85*Y)*sqrt(QRHO)
		N     = IDNINT(6 + 72*QRHO)
		J     = 2*N+1
		XSUM  = 1.0/J
		YSUM  = 0.0
		I=N;
		do
			J    = J - 2
			XAUX = (XSUM*XQUAD - YSUM*YQUAD)/I
			YSUM = (XSUM*YQUAD + YSUM*XQUAD)/I
			XSUM = XAUX + 1.0/J
			I = I-1
		while (I>0);
        U1   = -FACTOR*(XSUM*YABS + YSUM*XABS) + 1.0
        V1   =  FACTOR*(XSUM*XABS - YSUM*YABS)
        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 FRACTION
//  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 > 1.0) 
			H    = 0.0
			KAPN = 0
			QRHO = SQRT(QRHO)
			NU   = IDNINT(3 + (1442/(26*QRHO+77)))
		else
			QRHO = (1-Y)*SQRT(1-QRHO)
			H    = 1.88*QRHO
			H2   = 2*H
			KAPN = IDNINT(7  + 34*QRHO)
			NU   = IDNINT(16 + 26*QRHO)
		endif

		if (H > 0.0) 
			QLAMBDA = H2^KAPN;
		endif;
//*
		RX = 0.0
		RY = 0.0
		SX = 0.0
		SY = 0.0
//*
		N = NU;
		do
			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 ((H > 0.0) %& (N <= KAPN))
				TX = QLAMBDA + SX
				SX = RX*TX - RY*SY
				SY = RY*TX + RX*SY
				QLAMBDA = QLAMBDA/H2
			endif
			N = N-1;
		while (N>=0);
//*
		if (H == 0.0)
			U = FACTOR*RX
			V = FACTOR*RY
		else
			U = FACTOR*SX
			V = FACTOR*SY
		endif
//*
		if (YABS == 0.0) 
			U = EXP(-XABS^2)
		endif;
//*
	endif
//*
//*
//
//  EVALUATION OF W(Z) IN THE OTHER QUADRANTS
//
//*
	if (YI < 0.0)
//*
		if (QRHO < 0.085264)
			U2    = 2*U2
			V2    = 2*V2
		else
			XQUAD =  -XQUAD
//*
//
//         THE FOLLOWING IF-STATEMENT PROTECTS 2*EXP(-Z*Z)
//         AGAINST OVERFLOW
//
			if ((YQUAD > RMAXGONI) %| (XQUAD > RMAXEXP)) 
				return cmplx(inf,inf);
			endif;
//*
			W1 =  2*EXP(XQUAD)
			U2  =  W1*COS(YQUAD)
			V2  = -W1*SIN(YQUAD)
		endif;
//*
		U = U2 - U
		V = V2 - V
		if (XI > 0.0) 
			V = -V
		else
			if (XI < 0.0) 
				V = -V
			endif;
		endif
//*
	endif;
      RETURN cmplx(U,V);
end

function IDNINT(x)
	variable/d x;
	return round(x);
end;