But I'm not quite done with it yet...

As noted by Bastin and Heijligers in Electron Probe Quantitation, p. 160, the equations they give for beta/(2*alpha) only provide an approximation of the beta parameter. In order to maintain internal consistency, the value of this parameter must then be refined. Below is some Fortran 90/95 code (i = emitter and j = absorber) in which I implement a simple iterative procedure taken and slightly modified from GMRFILM that forces the integral of phi(rho*z) to be equal (within specified tolerance) to (R/S)*(1/Q(U0)) determined from PAP. I'm still fiddling with it, but I think it's working alright:

! Approximation of beta:

R_b_2a = ( gamma(j) - 2.D0 * alpha(i,j) * FZ(i,j) / SQRT( Pi ) ) &

/ ( gamma(j) - phi_0(i,j) )

IF ( ( R_b_2a >= 1.D0 ) .OR. ( R_b_2a < 0.D0 ) ) THEN

R_b_2a = 0.5D0

alpha(i,j) = ( ( phi_0(i,j) + gamma(j) ) * SQRT( Pi ) ) / ( 4.D0 * FZ(i,j) )

END IF

IF ( ( R_b_2a < 1.D0 ) .AND. ( R_b_2a >= 0.9D0 ) ) THEN

beta(i,j) = 0.9628832D0 - 0.9642440D0 * R_b_2a

ELSE IF ( ( R_b_2a < 0.9D0 ) .AND. ( R_b_2a >= 0.8D0 ) ) THEN

beta(i,j) = 1.122405D0 - 1.141942D0 * R_b_2a

ELSE IF ( ( R_b_2a < 0.8D0 ) .AND. ( R_b_2a >= 0.7D0 ) ) THEN

beta(i,j) = 13.43810D0 * EXP( -5.180503D0 * R_b_2a )

ELSE IF ( ( R_b_2a < 0.7D0 ) .AND. ( R_b_2a >= 0.57D0 ) ) THEN

beta(i,j) = 5.909606D0 * EXP( -4.015891D0 * R_b_2a )

ELSE IF ( ( R_b_2a < 0.57D0 ) .AND. ( R_b_2a >= 0.306D0 ) ) THEN

beta(i,j) = 4.852357D0 * EXP( -3.680818D0 * R_b_2a )

ELSE IF ( ( R_b_2a < 0.306D0 ) .AND. ( R_b_2a >= 0.102D0 ) ) THEN

beta(i,j) = ( 1.D0 - 0.5379956D0 * R_b_2a ) / ( 1.685638D0 * R_b_2a )

ELSE IF ( ( R_b_2a < 0.102D0 ) .AND. ( R_b_2a >= 0.056D0 ) ) THEN

beta(i,j) = ( 1.D0 - 1.043744D0 * R_b_2a ) / ( 1.604820D0 * R_b_2a )

ELSE IF ( ( R_b_2a < 0.056D0 ) .AND. ( R_b_2a >= 0.03165D0 ) ) THEN

beta(i,j) = ( 1.D0 - 2.749786D0 * R_b_2a ) / ( 1.447465D0 * R_b_2a )

ELSE IF ( ( R_b_2a < 0.03165D0 ) .AND. ( R_b_2a >= 0.0D0 ) ) THEN

beta(i,j) = ( 1.D0 - 4.894396D0 * R_b_2a ) / ( 1.341313D0 * R_b_2a )

END IF

beta(i,j) = beta(i,j) * ( 2.D0 * alpha(i,j) )

! Refinement of beta(i,j), adapted from GMRFILM by R. Waldo:

beta0 = beta(i,j)

y1 = beta(i,j) / ( 2.D0 * alpha(i,j) )

DO

beta1 = beta0

y2 = ARFC( y1 ) / R_b_2a * y1

beta0 = y2 * ( 2.D0 * alpha(i,j) )

IF ( ABS( ( beta1 - beta0 ) / beta1 ) < 1.D-4 ) EXIT

y1 = y2

END DO

beta(i,j) = y1 * ( 2.D0 * alpha(i,j) )

For mixtures, FZ(i,j) is replaced by FZ_mix(i), gamma(j) is replaced by gamma_mix, phi_0(i,j) is replaced by phi_0_mix(i), alpha(i,j) is replaced by alpha_mix(i), and beta(i,j) is replaced by beta_mix(i).

The ARFC function is equivalent to ZAFErrorFunction in CalcZAF and evaluates the polynomial that multiplies the exponential factor in the approximation for the complementary error function, as in Abramowitz and Stegun, Handbook of Mathematical Functions, p. 299, eq. 7.1.26.