is there any analytical solution here?

I'm sorry to say that I don't know whether analytical solution exists. It exists under some assumptions. It could also be quite easy to solve it numerically.

I'm going to start with the numerical solution as it is what you probably need. One could solve the system numerically as it is given, but the dimension of the problem can be easily reduced. That can be done by using multiplying the equations by $\sin(\theta)$ and $-\cos(\theta)$ respectively and adding them together which gives an equation of a single variable ($\theta$): $$\left[ \frac{a}{c}(\theta-z_0)+x_0\right]\sin(\theta)=\left[ \frac{b}{c}(\theta-z_0)+y_0\right]\cos(\theta) $$ This equation can already be solved numerically pretty easily by Newton's method for example.

It can also be rearanged to this form: $$\tan(\theta)=\frac{a(\theta-z_0)+cx_0}{b(\theta-z_0)+cy_0}= \frac{a\theta-az_0+cx_0}{b\theta-bz_0+cy_0}$$ This form can be used for numerical evaluation aswell though I doubt there is much advantage in terms of efficiency. But from this form it is clear we want to find the intersections of tangentoid and hyperbola. This can be used for graphical solution. It may also be used for spotting important roots, which can then be evaluated numerically.

It may also be of importance to you, that for a larger $|\theta|$ the solutions are closer to a solution of: $$\tan(\theta)\approx\frac{a}{b} \iff \theta\approx\arctan\left(\frac{a}{b}\right)+k\pi$$ It is easy to show: $$\lim_{\theta\to\infty}\frac{a\theta-az_0+cx_0}{b\theta-bz_0+cy_0}=\frac{a}{b}$$ This is nice because for any finite precision needed there is only finite number of solutions which need be solved numerically and the rest of large enough $|\theta|$ is given by this simple equation.

Finally this problem has analytical solution at least for the case the line goes through the origin. Then the aproximate expression above holds exactly: $$\tan(\theta)=\frac{a}{b} \iff \theta=\arctan\left(\frac{a}{b}\right)+k\pi$$

I wish you luck with your research and hope I have not done a mistake in my reasoning here. In the worst case I believe you'll get your answer in accordance with Cunningham's Law.

Strating from @Lukáš J.'s answer, consider the function

$$f(\theta)=\left(\frac ac(\theta -z)+x\right)\sin (\theta )-\left(\frac bc(\theta -z)+y\right)\cos (\theta )$$ for which you want the first zero.

For easier manipulations, let $$\alpha=\frac ac \qquad \beta=x-\frac acz\qquad \gamma=\frac bc \qquad \delta=y-\frac bcz$$ to write $$f(\theta)=(\alpha\, \theta+\beta)\sin (\theta )-(\gamma\, \theta+\delta)\cos (\theta )$$

Assuming that $\theta$ could be small, expand $f(\theta)$ as a series around $\theta=0$. This will give $$f(\theta)=-\delta + (\beta -\gamma )\theta +\frac{2 \alpha +\delta }{2} \theta ^2 +\frac{3 \gamma -\beta}{6} \theta ^3 -\frac{4 \alpha +\delta}{24} \theta ^4 +\frac{\beta -5 \gamma}{120} \theta ^5 +O\left(\theta ^6\right)$$ Now, using series reversion, the approximation

$$\theta =t-\frac{2 \alpha +\delta }{2(\beta - \gamma) }t^2+\frac{3 (2 \alpha +\delta )^2+(\beta -3 \gamma ) (\beta -\gamma )}{6 (\beta -\gamma )^2}t^3+\cdots +O\left(t^6\right)$$ where $t=\frac{f(\theta)+\delta }{\beta -\gamma }$; since we want $f(\theta)=0$, then $t=\frac{\delta }{\beta -\gamma }$