def direction(theta, phi):
return np.array([np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta),<EOL>np.cos(theta)])<EOL>
Return the direction vector of a cylinder defined
    by the spherical coordinates theta and phi.