function phi = phi_exact(P, sg_t, al, r)

% Find A by solving the equation that defines it.  It's easier to bound
% 1/A rather than A so I solve for 1/A instead.
invA = fzero(@(z) al / z * log((1 + z) / (1 - z)) - 2, [0.001,.99999]);
A = 1/invA;

% Compute P_d using the formula.
P_d = (2 * (A^2 - 1)) / (al * A^2 * (1 - A^2 * (1 - al)));

% Compute the fluence.  The integral is written as an integration from 1 to
% infinity, but I change variables to u = 1/t so as to get a finite
% interval [0, 1] for the numerical integration.
phi = P * sg_t / (4*pi * r) .* ...
      (P_d * exp(-sg_t * r / A) + ...
       quad(@(u) g(1./u, al) ./ u.^2 .* exp(-sg_t * r ./ u), 0, 1));

end

function g = g(t, al)

% Compute the g function using the formula.
g = 1 ./ ((1 - al./t .* atanh(1./t)).^2 + (pi/2 * al./t).^2);

end
