You are given that $\left\langle a+3b,2a-3b\right\rangle=0$ and $\left\langle a-4b,a+2b \right\rangle=0$. Hence
$$\left\langle a+3b,2a-3b\right\rangle=2\left\langle a,a\right\rangle-9\left\langle b,b\right\rangle+3\left\langle a,b\right\rangle=0$$ and
$$\left\langle a-4b,a+2b \right\rangle=\left\langle a,a \right\rangle-8\left\langle b,b \right\rangle-2\left\langle a,b \right\rangle=0.$$(Here I assumed that you are working with a real inner-product space).
By subtracting the second equation twice from the first, we obtain
$$\left\langle b,b \right\rangle+\left\langle a,b \right\rangle=0.$$
Plugging the previous equation into the second yields
$$\|a\|^2=6\|b\|^2.$$
Hence $\frac{\|b\|}{\|a\|}=\sqrt{\frac{1}{6}}.$
Now you know that $\left\langle a,b \right\rangle=\cos(\theta)\|a\|\|b\|.$
Thus, after plugging the previous results in the third equation, we get
$$\cos(\theta)=-\frac{\|b\|}{\|a\|}=-\sqrt{\frac{1}{6}}.$$
So, $\vec{p}_0$, $\vec{p}_u$, and $\vec{p}_v$ define a plane.
If we use $\hat{e}_u$ as the unit $u$ axis vector, then
$$\hat{e}_u = \frac{\vec{p}_u - \vec{p}_0}{\left\lVert \vec{p}_u - \vec{p}_0 \right\rVert} = \frac{\vec{u}}{\left\lVert\vec{u}\right\rVert} $$
If we do one step of Gram–Schmidt process, we can use $\vec{v} = \vec{p}_v - \vec{p}_0$ to obtain the other unit axis vector for the 2D plane:
$$\begin{aligned}
\vec{e}_v &= \vec{v} - \left(\hat{e}_u \cdot \vec{v}\right) \hat{e}_u \\
\hat{e}_v &= \frac{\vec{e}_v}{\left\lVert\vec{e}_v\right\rVert} \\
\end{aligned}$$
We can extend this trivially to a 3D basis via $$\hat{e}_w = \hat{e}_u \times \hat{e}_v$$
Unit vectors $\hat{e}_u = (u_x, u_y, u_z)$, $\hat{e}_v = (v_x, v_y, v_z)$, and $\hat{e}_w = (w_x, w_y, w_z)$ now form an orthonormal basis for our 2D (and 3D) coordinate system corresponding to the plane. On the plane, $w = 0$.
Point $(u, v, w)$ on the plane coordinate system corresponds to point $\vec{p} = (x, y, z)$ in the original 3D coordinate system via
$$\vec{p} = \vec{p}_0 + u \hat{e}_u + v \hat{e}_v + w \hat{e}_w \quad \iff \quad \left\lbrace \begin{aligned}
x &= x_0 + u u_x + v v_x + w w_x \\
y &= y_0 + u u_y + v v_y + w w_y \\
z &= z_0 + u u_z + v v_z + w w_z \\
\end{aligned} \right .$$
which we can solve for $u$, $v$, and $w$. (Note that $\det \left[ \begin{matrix} \hat{e}_u & \hat{e}_v & \hat{e}_w \end{matrix} \right] = 1$, because the three vectors are orthonormal – orthogonal and of unit length – which simplifies the solution quite a bit.)
$$\left\lbrace ~ \begin{aligned}
u &= (x - x_0)(v_y w_z - v_x w_y) + (y - y_0)(v_z w_x - v_x w_z) + (z - z_0)(v_x w_y - v_y w_x) \\
v &= (x - x_0)(u_z w_y - u_y w_z) + (y - y_0)(u_x w_z - u_z w_x) + (z - z_0)(u_y w_x - u_x w_y) \\
w &= (x - x_0)(u_y v_z - u_z v_y) + (y - y_0)(u_z v_x - u_x v_z) + (z - z_0)(u_x v_y - u_y v_x) \\
\end{aligned} \right .$$
Note that given an arbitrary point $(x, y, z)$, $(u, v)$ is the location on the plane, and $w$ is the (signed) distance to the plane. In other words, if you ignore the $w$ coordinate (by treating it as zero), you project the point to the closest point on the plane.
(Note that you can define three vectors, $\vec{u}_e$ and $\vec{v}_e$ (and optionally $\vec{w}_e$), so that $u = \vec{u}_e \cdot (\vec{p} - \vec{p}_0)$ and $v = \vec{v}_e \cdot (\vec{p} - \vec{p}_0)$, (and optionally) $w = \vec{w}_e \cdot (\vec{p} - \vec{p}_0)$), speeding up the 3D-to-2D conversion process.)
In the 2D coordinate system, we can define $\theta = \operatorname{atan2}(v, u)$ (using the two-argument version of arcus tangent, covering full $360°$ or $2 \pi$ in radians, instead of half the plane as $\arctan(v/u)$ would).
Because of the way we chose $\hat{e}_u$, $\vec{p}_u$ is at $\left(\left\lVert\vec{u}\right\rVert, 0\right)$ in the 2D coordinates; i.e., on the positive $u$ axis. This means that in the 2D coordinate system, its angle $\theta_u = 0$. For $\vec{p}_v$ and $\vec{p}_s$, we need to use the above formula to find their 2D coordinates, and then their respective angles (wrt. the positive $u$ axis), $\theta_v$ and $\theta_s$.
This, finally, leads to a solution of the original question. The angle $\varphi$ between the two vectors $\vec{u}$ and $\vec{v}$, on the same side as point $\vec{p}_s$, projected to the plane passing through points $\vec{p}_0$, $\vec{p}_0+\vec{u}$, and $\vec{p}_0+\vec{v}$, is
$$\varphi = \left\lbrace ~ \begin{aligned}
\theta_v, & \quad 0 \le \theta_s \le \theta_v \\
2 \pi - \theta_v, & \quad \theta_s \lt 0 \le \theta_v \\
2 \pi - \theta_v, & \quad 0 \le \theta_v \lt \theta_s \\
-\theta_v, & \quad \theta_v \le \theta_s \le 0 \\
2 \pi + \theta_v, & \quad \theta_v \le 0 \lt \theta_s \\
2 \pi + \theta_v, & \quad \theta_s \lt \theta_v \le 0 \\
\end{aligned} \right . $$
assuming a typical implementation of $\operatorname{atan2}$ that yields results $\pm \pi$. (To convert from radians to degrees, multiply by $180 / \pi$.)
The first three cases are when $\theta_v$ is in the positive $v$ half of the 2D plane. The first case is $\vec{p}_s$ between the two angles, the second case is $\vec{p}_s$ in the negative $v$ half of the 2D plane, and the third case is $\vec{p}_s$ in the positive $v$ half of the 2D plane but outside $\vec{v}$.
The last three cases are when $\theta_v$ is in the negative $v$ half of the 2D plane, so we need to negate $\theta_v$ in the result. The fourth case is $\vec{p}_s$ between the two angles. The fifth case is $\vec{p}_s$ in the positive $v$ half of the 2D plane. The last case is $\vec{p}_s$ being in the negative $v$ half of the 2D plane, but outside $\vec{v}$.
Best Answer
You could apply the dot product formula $a \cdot b=|a|\ |b| \cos(30)$ where $|a|$ and $|b|$ are the lengths of the two vectors (square root of sum of squares of components). The vector $b$ is not determined by your condition but is only restricted to lie in a cone making all 3-d locations of vectors $b$ with axis $a$ and apex of cone angle 30 degrees.