Is this formula for percentile from z-score correct? Preface: I am not a mathematician. I have a basic understanding of statistics but am not familiar with mathematical notation beyond high school algebra.
I am a software developer maintaining medical records forms for a clinic. We have a form that calculates pediatric patients' height, weight, BMI, and other growth-related percentiles. In the process of troubleshooting a probably unrelated issue with the percentile calculations, I encountered this function that I do not understand.
This function takes only the z-score as an argument and returns a number percentile. I was able to confirm that the method to derive the z-score is correct (it is based on the LMS method outlined in this nature.com article), but I don't understand what is going on in the percentile calculation.
Code in MEL, a scripting language created in the mid 90's by GE:
local p10
local p9
local p4
local p3
local p2
local p1
local p

p10 = 1 + (0.33267 * AbsValue(zScore))
p9 = 1 / p10
p4 = (0.4361836 * p9) - (0.1201676 * (p9 ^ 2)) + (0.937298 * (p9 ^ 3))
p1 =  2 * 3.14159265
p3 = (AbsValue(zScore)^2)
p2 = -1 * p3 / 2

local pp
pp = 1 / (p1 ^ 0.5) * (2.71828182845905 ^ p2) * p4
p = 1 - pp

if zScore > 0 then
    p = p * 100
else
    p = 100 - p * 100
endif

p = div((p + 0.5), 1)

return p

Is the calculation correct? (I know it's at least close, I can cross-check against a table.) If it is correct, is it specific to the LMS growth chart tables (from the CDC/WHO) in some way, or can it be generalized? Can anyone offer any insight into what is going on here? Even if it is all correct and should not be changed, I would like to name the variables more descriptively and add some comments about what calculations are being done and why.
I have been searching on this question all day and can find only either a) explanations on how to find a percentile by looking a z-score up on a table and b) explanations involving formulas full of unfamiliar notation that I don't even know where to start to understand.

*

*I don't think the issue I was troubleshooting is related to this code, but it's hard to be sure without understanding the parts.

*I found this question but looking through the links in the answer and comments there did not lead me to anything I found useful.

 A: 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.
