Seriously, it has no solution? Gosh, am I surprised. So what does this do?
eqn = (a*sin(theta0)-b*sin(thetaz))./(a*cos(theta0)+b*cos(thetaz)) == z;
thsol = solve(eqn,thetaz,'returnconditions',true)
So there are infinitely many solutions. Each a multiple of 2*pi away from each other. We can arbitrarily choose the primary solutions, thus with k==0.
thsol = subs(thsol.thetaz,k,0)
Still not very interesting, but now if we do this:
thsol = matlabFunction(thsol);
now we find three trivially easy to obtain solutions, two of them are real, one is complex.
But I guess no analytical solution exists, because you said so. :)