It's an approximation as also admitted by your statement "satisfactory performance". Here is why.
Let $a[n]=\cos(w_0 n)$ and $b[n]=\cos(w_0 n + \theta)$ be two sequences of size $N$, represented by $1 \times N$ row matrices $\bar{a}$ and $\bar{b}$ in a program (such as Octave or Matlab)
Treating those two sequences $a[n]$ and $b[n]$ which are represented by two row matrices $\bar{a}$ and $\bar{b}$ as components of vectors $\vec{a}$ and $\vec{b}$, we can take advantage of the dot product stated for vectors and (approximately) compute the phase angle $\theta$ in between those two sinusoidal sequences $a[n]$ and $b[n]$
Geometric evaluation of the dot product between two vectors of same size is:
$$ \vec{a} \cdot \vec{b} = |\vec{a}||\vec{b}| \cos(\phi)$$
And the angle between those two vector therefore is:
$$\cos(\phi) = \frac{\vec{a}\cdot \vec{b}}{|\vec{a}||\vec{b}|}$$
Where $\phi$ is the angle between those vectors, which will be equal to $\theta$ as shown:
The dot product between the vectors $\vec{a}$ and $\vec{b}$ can also be algebrically computed from their MAC sum of components which are contained in the elements of the matrices $\bar{a}$ and $\bar{b}$:
$$ \vec{a} \cdot \vec{b} = \sum_{i=1}^{N} \bar{a}(i)\bar{b}(i) $$
The product $\bar{a}(i)\bar{b}(i)$ can be shown to be equal to the following (using trigonometry):
$$\bar{a}(i)\bar{b}(i) = \cos(w_0 i) \cos(w_0 i + \theta)= 0.5\cos(2 w_0 i + \theta) + 0.5 \cos(\theta)$$
And the dot-product summation becomes:
$$ \vec{a} \cdot \vec{a} = \sum_{i=1}^{N} a(i)b(i) = 0.5 N \cos(\theta) + \sum_{i=1}^{N} 0.5 \cos(2 w_0 i + \theta)$$
The absolute values $|\vec{a}|$ and $|\vec{b}|$ of the vectors $\vec{a}$ and $\vec{b}$ can also be computed from their norms computed over the matrices $\bar{a}$ and $\bar{b}$ as follows:
$$ |\vec{a}| = ||\bar{a}||= ( \sum_{i=1}^{N} a(i)^2 )^{0.5}$$
$$ |\vec{b}| = ||\bar{b}||= ( \sum_{i=1}^{N} b(i)^2 )^{0.5}$$
Again expanding the squared terms result in:
$$ |\vec{a}| = ||\bar{a}||= ( 0.5 \sum_{i=1}^{N} 1 + \cos(2 w_0 i) )^{0.5}$$
$$ |\vec{b}| = ||\bar{b}||= ( 0.5 \sum_{i=1}^{N} 1 + \cos(2 w_0 i + \theta) )^{0.5}$$
Now, it can easily be shown that the following summations are small compared to large $N$:
$$ N+\sum_{i=1}^{N} \cos(2 w_0 i) \approx N , ~~~ N+\sum_{i=1}^{N} \cos(2 w_0 i + \theta) \approx N $$
for large $N$ and suitable $w_0$:
Simplify the above equations to get the angle:
$$\cos(\phi) = \frac{\vec{a}\cdot \vec{b}}{|\vec{a}||\vec{b}|}$$
$$\cos(\phi) \approx \frac{ 0.5 N \cos(\theta)}{|\vec{a}||\vec{b}|}$$
$$\cos(\phi) \approx \frac{ 0.5 N \cos(\theta)}{0.5 N}$$
$$\cos(\phi) \approx \cos(\theta)$$
$$\phi \approx \theta$$
The approximation depends on the length $N$.
