From Cartesian to Spherical and back Cartesian again

I was writing a Python class to represent an atomic object. The various attributes are id, type, x, y, z and such. Inside the class I defined three methods to obtain the spherical coordinates.
While doing some tests, I discovered something crazy.
If you use the usual definitions for coordinates transformations (like these), take a single point of coordinates (1, 1, 1). Now transform it in spherical coordinates. And then back to Cartesian. Everything ok.
Now take (-1, 1, 1) or (1, -1, 1) or (-1, -1, 1). And have fun discovering that, depending on the choice for the azimuthal angle one of the coordinates (either x or y) will return with reversed sign.

Try it.

It took me some time to (re)discover that I need to change the phi definition according to the value of x and y.

Here's a simple python code on how it should be done:

import numpy as np
import sys

try:
 x= float(sys.argv[1]); y=float(sys.argv[2]) ; z=float(sys.argv[3])
except:
 print ("Error\nUsage:", sys.argv[0], "x y z"); sys.exit(1)

print ("you inserted\nx\ty\tz\n%f\t%f\t%f\n" %(x,y,z))

rho = np.sqrt((x*x)+(y*y)+(z*z))
theta= np.arccos(z/rho)
phi = np.arctan(y/x)

if y<0 and x>=0:
 phi= np.arcsin(y/np.sqrt(x*x+y*y))
elif x<0 and y>=0:
 phi=np.arccos(x/(np.sqrt(x*x+y*y)))
elif x<0 and y<0:
 phi+=np.pi

print ("rho\ttheta\tphi\n%f\t%f\t%f\n" %(rho, theta, phi))

x1= rho*np.sin(theta)*np.cos(phi)
y1= rho*np.sin(theta)*np.sin(phi)
z1= rho*np.cos(theta)

print ("x1\ty1\tz1\n%f\t%f\t%f\n" %(x1, y1, z1))

The cool part is to add pi to the arctan definition when both x and y are negative, and to use the definitions with arcsin and arccos when only y or only x is negative. If you do like that, things never loose sign.
I don't know if it is me or what, but I find irritating that nobody ever told me that if you go like Cartesian->Spherical->Cartesian you have to be so careful.
Well, now it works.