I want to calculate Mercator projection of points inside of sphere. First, the points are transformed to local system. Then the latitude and longitude are calculated based on following equations
$\Large r=\sqrt{x^2+y^2+z^2}\\\Large\theta=\cos^{-1}\frac zr=\cos^{-1}\frac{z}{x^2+y^2+z^2}\\\Large\psi=\tan^{-1}\frac yx$
and the Mercator projection for each point is calculated based on:
$\Large(x,y)=(\lambda, \ln[\tan(\frac{\phi}{2}+\frac{\pi}{4})])$
but it has error in some points because of $\ln$ of negative numbers and when I check the latitudes in these cases, they are more than $90^\circ$. What's wrong?
In your second equation, you compute the coordinates $(x,y$) of the projected point, that is, on the planar map. Here $y$ can be negative, but there is no log of a negative value, and latitudes $\varphi\in]-\pi/2,\pi/2[$ (in radians, not degrees) are mapped to $y\in]-\infty,+\infty[$.
The problem lies in your first system of equations. In this system, the $\color{red}{\mathrm{co}}$latitude is $\theta$, and the longitude $\varphi$. In cartography the longitude is $\lambda$, and the latitude is $\varphi$, so I will use $\varphi$ for the latitude here.
First, be aware that the notation $\arctan \frac yx$ should really be $\varphi=\mathrm{atan2}(y,x)$ (see the link on Wikipedia, where you picked the equations). Here Wikipedia perpetuates a long tradition of incorrectly writing the formula: the correct value of the longitude is not $\arctan \frac yx$ in general. (additional caution: some programming languages write it
atan2(y,x), othersatan2(x,y), so you will have to check)But the problem really is the latitude: in mathematics and physics, it's common to measure the colatitude, that is the angle between the $z$ axis and the vector $\vec r$, which is then in $[0,\pi]$: that's what is given on Wikipedia, and it's clearly stated in the article that "$\theta$ is inclination from the $z$ direction". In cartography, you want the latitude, measured from the equator. Here the latitude is $\varphi=\frac\pi2-\theta$.