Is this the correct solution for finding the samples of phase delay at a given frequency of an IIR BiQuad filter? I want to calculate the number of samples of phase delay at a given frequency of an IIR BiQuad Filter. I believe I have solved this, but I would like to confirm I have done it correctly. My work is as follows.

Background:

The transfer function of a BiQuad filter is:

$$
H(z) = g \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{a_0 + a_1z^{-1} + a_2z^{-2}}
$$
ω is normalized by the sample rate (samples/sec) with units of radians/sample. The normalized Nyquist frequency is π radians/sample and the normalized sample rate is 2π radians/sample.

$$
\omega = \frac{2 \pi f}{f_s}
$$
The frequency response of a system can be calculated by letting $z = e^{j\omega}$, where ω is a normalized frequency in the range (−π,π):

$$
H(e^{j\omega}) = g\frac{b_0 + b_1e^{-j\omega} + b_2e^{-j2\omega}}{a_0 + a_1e^{-j\omega} + a_2e^{-j2\omega}}
$$

The value of $H(e^{j\omega})$ will be complex. The phase response at the frequency ω is just the phase angle of the resulting complex number. The magnitude response at the same frequency is likewise equal to the magnitude of the number.

Euler's Formula for solving complex exponents is:
$$e^{ix}=cos(x)+isin(x)$$

Solution:
$$
H(e^{j\omega}) = g \frac{b_0 + b_1e^{-j\omega} + b_2e^{-j2\omega}}{a_0 + a_1e^{-j\omega} + a_2e^{-j2\omega}}
$$
$$
H(e^{j\omega}) = g \frac{b_0 + b_1(cos(-ω) + jsin(-ω)) + b_2(cos(-2ω) + jsin(-2ω))}{a_0 + a_1(cos(-ω) + jsin(-ω)) + a_2(cos(-2ω) + jsin(-2ω))}
$$
$$
H(e^{j\omega}) = g \frac{b_0 + b_1cos(ω) + b_1jsin(-ω) + b_2cos(2ω) + b_2jsin(-2ω)}{a_0 + a_1cos(ω) + a_1jsin(-ω)) + a_2cos(2ω) + a_2jsin(-2ω)}
$$
$$
H(e^{j\omega}) = g \frac{[b_0 + b_1cos(ω) +b_2cos(2ω)] + j[b_1sin(-ω) +  b_2sin(-2ω)]}{[a_0 + a_1cos(ω) + a_2cos(2ω)] + j[a_1sin(-ω))  + a_2sin(-2ω)]}
$$
$$
H(e^{j\omega}) = \frac{[g(b_0 + b_1cos(ω) +b_2cos(2ω))] + j[g(b_1sin(-ω) +  b_2sin(-2ω))]}{[a_0 + a_1cos(ω) + a_2cos(2ω)] + j[a_1sin(-ω))  + a_2sin(-2ω)]}
$$

Simplify:

$$m = g(b_0 + b_1cos(ω) +b_2cos(2ω))$$

$$n = g(b_1sin(-ω) +  b_2sin(-2ω))$$

$$p = a_0 + a_1cos(ω) + a_2cos(2ω)$$

$$q = a_1sin(-ω))  + a_2sin(-2ω)$$

$$H(e^{j\omega}) =\frac{m+in}{p+iq}$$

Find the Angle (delay in samples?):

For a given complex equation of:

$z= x+iy$

Where φ is the angle between the x axis and the vector z measured counterclockwise in radians:

$φ = atan2(y, x)$

Additionally, the phase of the ratio of two complex numbers is the difference of the phases of the numerator and denominator (numerator - denominator), as noted by Matt L in his answer below:

$$φ(ω) = atan2(n,m) - atan2(q,p)$$

To get the final phase delay in samples, we must take the negative of the phase we have solved and then divide that by the normalized frequency:

$$\tau_p(\omega)=-\frac{φ(\omega)}{\omega}$$

Therefore the final equation for phase delay for any normalized angular frequency ω would be:

$$\tau_p(\omega)= \frac{-atan2(n,m) + atan2(q,p)}{\omega}$$

(Updated from Matt L's feedback.)
You almost got it, but the final result is wrong because the phase of the ratio of two complex numbers is the difference of the phases of the numerator and denominator, respectively, not the ratio of the phases:

$$z=\frac{u}{v}=\frac{|u|e^{j\phi_u}}{|v|e^{j\phi_v}}=\frac{|u|}{|v|}e^{j(\phi_u-\phi_v)}\Longrightarrow \phi_z=\phi_u-\phi_v$$

Furthermore, for finding the phase delay, you need to divide the phase by the (normalized) frequency in radians (as shown in this answer). The result is the phase delay in samples.