Jean Pellegri suggested another, simpler, approach, namely to solve for
tan(x):
(%i2) load(to_poly_solver)$
...
(%i3) eq: 2*cos(x)^2*x0-2*sin(x)*cos(x)*y0=x0;
(%o3) 2*cos(x)^2*x0-2*cos(x)*sin(x)*y0 = x0
(%i4) subst(atan(q),x,%);
(%o4) 2*x0/(q^2+1)-2*q*y0/(q^2+1) = x0
(%i5) solve(%,q);
(%o5) [q = -(sqrt(y0^2+x0^2)+y0)/x0,q = (sqrt(y0^2+x0^2)-y0)/x0]
(%i6) subst(tan(x),q,%);
(%o6) [tan(x) = -(sqrt(y0^2+x0^2)+y0)/x0,tan(x) = (sqrt(y0^2+x0^2)-y0)/x0]
To put this in terms of x:
(%i7) map(atan,%);
(%o7) [atan(tan(x)) = -atan((sqrt(y0^2+x0^2)+y0)/x0),
atan(tan(x)) = atan((sqrt(y0^2+x0^2)-y0)/x0)]
Note that Maxima is cautious, and does not simplify atan(tan(x)) to x; For
example, atan(tan(%pi)) is 0, not %pi.
-s