Expanded from comments:

The function you're after is the (standard) normal cumulative distribution function. It's not specific to the application domain you mention, it's just the probability that a standard normal random variable is less than or equal to some fixed value. You're answering the question: if these "z-scores" were drawn from a standard normal distribution, what would be the probability of observing a z-score which is smaller or equal to this one?

This function is defined mathematically by an integral:

$$P(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^x e^{-y^2/2}dy$$

This has no closed form. That is, there is no exact formula for $P(x)$ involving only usual operations like sums and products or functions like exponentials, logarithms, trigonometric functions, etc. To implement this function in software, you need an approximate solution which you can actually compute and which is close enough to the true solution for your needs.

This function is ubiquitous in many fields and the problem of approximating it in software has been extensively studied. The classic reference on the topic of approximating mathematical functions is Abramowitz & Stegun (1964), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.

The specific approximation you show is given there as equation 26.2.16. Note that the approximation is defined only for $x \geq 0$, but you actually want to be able to compute $P(x)$ for any real $x$. You can extend it for negative $x$ by taking advantage of the symmetry of the distribution, which implies (26.2.5-26.2.6):

$$P(-x) = 1-P(x)$$

So, if the input is negative, flip the sign, compute the approximation on that, then return one minus the result (that's what AbsValue and the if/else at the end are doing in your implementation).

In their notation, they first define (the normal probability density function):

$$Z(x) = \frac{1}{\sqrt{2\pi}}e^{-x^2/2}$$

and the transformation:

$$t=\frac{1}{1+px}$$

The approximation is then given by:

$$1-Z(x)(a_1t+a_2t^2+a_3t^3)$$

$p$, $a_1$, $a_2$, $a_3$ are specific numerical constants given in the text. They have been carefully chosen so as to minimize some error criterion between the approximation and the true value. An absolute error bound is also given:

$$|\epsilon(x)| = |P(x) - \left[1-Z(x)(a_1t+a_2t^2+a_3t^3)\right]| < 10^{-5}$$

That is, this approximation always yields an answer which is within $10^{-5}$ of the true value. Note that this is an absolute bound, and it may not be very useful to have such a bound if you are looking at very small probabilities. More modern statistical software will use a more accurate approximation in general, or combine several approximations where each is better adapted to different ranges of the input parameter to get accurate results across all possible inputs.