Express $(\theta, \phi)$ of a torus using $(x,y,z)$ I know how to write the $(x,y,z)$ coordinates of a torus using $(\theta, \phi)$:
$$
\begin{align}
x(\theta, \phi) &= (R+r\cos(\theta))\cos(\phi) \\
y(\theta, \phi) &= (R+r\cos(\theta))\sin(\phi) \\
z(\theta, \phi) &= r\sin(\theta)
\end{align}
$$
For $\theta, \phi\in [0, 2\pi]$. How can I find the inverse relation? I would like to find expressions
$$
\theta(x,y,z) = \ldots \\
\phi(x,y,z)=\ldots
$$

Given that I have a point on the torus in cartesian coordinates $(x, y,z)$. How can I find its expression in terms of angles $(\theta, \phi)$?

From the last expression I can get
$$
\theta = \arcsin\left(\frac{z}{r}\right)
$$
and then I guess from the first equation one could find
$$
\phi = \arccos\left(\frac{x}{R + r\cos(\theta)}\right)
$$
which works when the denominator is not zero. To avoid having the denominator being zero we need
$$
R+r\cos(\theta) \neq 0 \implies \theta \neq \arccos\left(-\frac{R}{r}\right)
$$
 A: Note that you have $$ \sqrt{x^2+y^2} = R+r\cos\theta$$
We have then
$$ \cos\phi = \frac{x}{\sqrt{x^2+y^2}} $$
$$ \sin\phi = \frac{y}{\sqrt{x^2+y^2}} $$
$$ \cos\theta = \frac{\sqrt{x^2+y^2}-R}{r}$$
$$ \sin\theta = \frac{z}{r}$$
Which can be used to recover $\theta$ and $\phi$. You should be careful with using $\arcsin$ though, at it will give you the correct value of the angle only if said angle is in the interval $[-\frac{\pi}{2},\frac{\pi}{2}]$.
To deal with that you need to consider different cases. That is because the equation $\sin\alpha = t$ has many solutions; one is $\alpha=\arcsin t$, but $\alpha=2\pi n+\arcsin t$, $n\in\mathbb{Z}$ is also a solution, as is $\alpha=2\pi n+\pi-\arcsin t$, $n\in\mathbb{Z}$. Out of all of these solutions the correct one is the one that gives the correct sign of $\cos \alpha$ and gives angle $\alpha$ in the correct predefined interval. For example, if you want $\phi$ to be in the interval $[0,2\pi)$, the full formula for $\phi$ will be
$$ \phi = \left\{\begin{array}{ll}\arcsin\frac{y}{\sqrt{x^2+y^2}} & \text{if } x\ge 0,y\ge 0 \\ \pi-\arcsin\frac{y}{\sqrt{x^2+y^2}} &\text{if } x \le 0\\ 2\pi+\arcsin\frac{y}{\sqrt{x^2+y^2}} &\text{if } x\ge 0, y<0\end{array}\right. $$
On the other hand, if you wanted the values of $\phi$ to be in the range $(-\pi,\pi]$, the formula would be
$$ \phi = \left\{\begin{array}{ll}\arcsin\frac{y}{\sqrt{x^2+y^2}} & \text{if } x\ge 0 \\ \pi-\arcsin\frac{y}{\sqrt{x^2+y^2}} &\text{if } x \le 0, y\ge 0\\ -\pi-\arcsin\frac{y}{\sqrt{x^2+y^2}} &\text{if } x\le 0, y<0\end{array}\right. $$
It can be checked that if $\phi$ is definied like this, then $\sin\phi$ and $\cos\phi$ satisfy the previous equations regardless of the signs of $x$ and $y$. $\theta$ can be properly defined in an analogous way, with different formulas depending on the signs of $z$ and $\sqrt{x^2+y^2}-R$.
A: If you remember what $r$ and $R$ are (inner and outer radius), you see that $R > r$ and so your additional assumption is automatically satisfied:
$$ 
R > r 
\quad \Longrightarrow \quad 
R > r \cos \theta 
\quad \Longrightarrow \quad 
R - r \cos \theta \neq 0 
\quad \text{for any } \theta.
$$
To be precise, these formulas only work in some region, in the same way that $\alpha = \arcsin(\sin \alpha)$ is only sometimes true. To fix this, you can either treat each case separately or compute both $\sin$ and $\cos$ for the angles you need (which is better - see the explanation below and Adam Latosiński's answer).

Added later. A general remark on using $\arcsin$ and $\arccos$, especially in coding context. Knowing $\cos \alpha$ or $\sin \alpha$ is not enough to recover $\alpha$. You could artificially divide the domain into several parts and use a different formula (involving $\arcsin$ or $\arccos$) in each region, but even then, the method is numerically bad.
In contrast, if you know both $\sin \alpha$ and $\cos \alpha$, the angle $\alpha$ is uniquely determined (up to a multiple of $2\pi$). Programming languages provide two-variable functions that allow you to recover $\alpha$ from these two numbers, see here for Python.
A: Your torus $T \subset \mathbb R^3$ is parameterized by
$$f : [0,2\pi] \times [0,2\pi] \to \mathbb R^3, f(\theta,\phi) = ((R+r\cos(\theta))\cos(\phi), (R+r\cos(\theta))\sin(\phi), r\sin(\theta)) .$$
Note that this function is not injective. We have $f(\theta,\phi) = f(\theta',\phi')$ iff both $\theta - \theta', \phi - \phi' \in \{0, 2\pi\}$. This reflects the "usual" construction of $T$ by identifying opposite edges of the square $[0,2\pi] \times [0,2\pi]$: First identify $[0,\phi]$ with $[2\pi,\phi]$ which yields a tube. In this tube identify $[\theta,0]$ with $[\theta,2\pi]$ which yields a torus.
You can consider the restriction $g = f\mid_{[0,2\pi) \times [0,2\pi)}$ to get a bijective parameterization of $T$. You want to determine $g^{-1}$.
In his answer Adam Latosiński has given four equations which uniquely determine $\theta,\phi \in [0,2\pi)$ in terms of $x, y, z$. Unfortunately $\sin, \cos : [0,2\pi] \to [-1,1]$ are not injective, thus we cannot naively take the functions $\arcsin, \arccos$ to compute  $\theta,\phi$. In fact we have bijections $\arcsin : [-1,1] \to [-\frac{\pi}{2},\frac{\pi}{2}]$ and $\arccos : [-1,1] \to [0,\pi]$ (principal values). Using them correctly gives ugly formulae involving many cases. For example we get
$$\phi  = \begin{cases} \phantom{-}\arccos(\frac{x}{\sqrt{x^2+y^2}}) & y \ge 0 \\ 2\pi -\arccos(\frac{x}{\sqrt{x^2+y^2}}) &  y \le 0 \end{cases}$$
It is much better to work with complex numbers. We know that $\exp : [0,2\pi) \to S^1 =  \{ w \in \mathbb C \mid \lvert w  \rvert = 1\}$ is a bijection. Since $\exp(t) = \cos(t) + i \sin(t)$ we get from Adam Latosiński equations (with $w = x + iy$)
$$\exp(\phi) = \frac{w}{\lvert w  \rvert} ,$$
$$\exp(\theta) = \frac{\lvert w  \rvert - R +iz}{r} .$$
Update:
An alternative approach is this.
We can identify $\mathbb C$ and $\mathbb R^2$. Doing so, $S^1 = \{ (u,v) \in \mathbb R^2 \mid u^2 + v^2 = 1\}$.
The function $\exp : [0,2\pi] \to S^1, \exp(t) =(\cos(t), \sin(t))$, is a surjective closed map between compact Hausdorff spaces, hence a quotient map. This implies that $f$ induces a homeomorphism $\bar f : S^1 \times S^1 \to T$ given by
$$\bar f((a,b),(c,d)) = ((R+rc)a, (R+rc)b, rd).$$
Its inverse $\bar f^{-1} : T \to S^1 \times S^1$ is given by
$$a(x,y,z) = \frac{x}{\sqrt{x^2 + y^2}} ,$$
$$b(x,y,z) = \frac{y}{\sqrt{x^2 + y^2}} ,$$
$$c(x,y,z) = \begin{cases} \phantom{-}\sqrt{1 - (\frac{z}{r})^2} & x^2 +y ^2 \ge R^2 \\ -\sqrt{1 - (\frac{z}{r})^2} &  x^2 +y ^2 \le R^2 \end{cases}$$
$$d(x,y,z) = \frac{z}{r} .$$
Note that  $x^2 +y ^2 = R^2$ means $z = \pm r$ for $(x,y,z) \in T$.
This gives you the "$S^1$-coordinates" of $\bar f^{-1}(x,y,z)$. But assocating to a point $(u,v) \in S^1$ an angle $\psi$ such that $(cos(\psi), \sin(\psi)) = (u,v)$ is easy. Either take
$$\psi  = \begin{cases} \phantom{-}\arccos(u) & v \ge 0 \\ 2\pi -\arccos(u) &  v \le 0 \end{cases}$$
or
$$\psi  = \begin{cases} \phantom{-}\arcsin(v) & u \ge 0 \\ \pi - \arcsin(v) &  u \le 0 \end{cases}$$
Note that the first formula gives you $\psi \in [0,2\pi]$ and the second $\psi \in [-\frac{\pi}{2},\frac{3\pi}{2}]$.
A: I'm going to answer the question you meant to ask (and indeed, did ask in a comment to Michał Miśkiewicz's answer), namely:

*

*And basically all I want to do now is, given my point (0,0,0) I want to find (0,0) so that I can plug them in the two expressions above to find the basis vectors of the tangent plane in terms of (,,)

I'm going to first rephrase that:
Given a point $P = (x, y, z)$ on a torus with main radius $R$ and tube radius $r$ (both known), find "the" basis for the tangent plane at $P$, where I mean the basis where one arrow points "around the tube" and the other points "around the big circle".
I confess that I'm only conjecturing that this is the basis you want, but it is the one given by the vectors $\frac{\partial{x,y,z}}{\partial \phi}$ and $\frac{\partial{x,y,z}}{\partial \theta}$, so that seems likely.
The second one is easiest: a vector pointing "around the big circle" at location $(x, y, z)$ is given by
$$
{\mathbf v}_\theta = \pmatrix{-y\\ x\\ 0}.
$$
That's easy to see if you look from above (i.e., along the $z$-axis) and consider the problem as being just in the $xy$-plane.
The first is slightly messier, but not much. If we let $\mathbf u$ be the unit vector pointing from the origin to $(x, y, 0)$, then the "tube slice" at $P$ lies in the $uz$ plane (with positive $u$ coordinate!). In the $uz$-coordinate system, your point $P$ lies at $(u, z) = (R + r\cos \phi, r \sin \phi)$ and by taking the derivative of this with respect to $\phi$, we get the vector (again in $uz$ coordinates)
$$
\pmatrix{-r\sin \phi \\ r \cos \phi}.
$$
The first coordinate here is just the negative $z$ coordinate of our point $P$; the second is the $u$-coordinate of $P$ (namely $\sqrt{x^2 + y^2}$), minus $R$, divided by $r$. So because
$$
\mathbf u = \frac{1}{\sqrt{x^2 + y^2}}\pmatrix{x\\y\\0},
$$
we get that
\begin{align}
\mathbf v_\phi 
&= -z\mathbf u + \left(\frac{\sqrt{x^2 + y^2} - R}{r}\right)\hat{\mathbf z}\\
&= \frac{-z}{\sqrt{x^2 + y^2}}\pmatrix{x\\y\\0} + \left(\frac{\sqrt{x^2 + y^2} - R}{r}\right)\pmatrix{0\\0\\1}\\
\end{align}
Because you didn't ask for a unit vector, I guess I can multiply through by the constant $\sqrt{x^2 + y^2}$ and say instead that
\begin{align}
\mathbf v_\phi 
&= {-z}\pmatrix{x\\y\\0} + \left(\frac{(x^2 + y^2) - R\sqrt{x^2 + y^2}}{r}\right)\pmatrix{0\\0\\1}\\
&= \pmatrix{-zx\\-zy\\ \frac{(x^2 + y^2) - R\sqrt{x^2 + y^2}}{r}}
\end{align}
instead.
Notice that the numbers $\theta$ and $\phi$ don't enter into this anywhere, so you don't need to worry about the domain of arctan, arccos, etc.
