How to rotate a 3D spot with numpy?

Question

I am trying to distribute spots on a sphere. There's a function to define a spot with rho, theta, phi its spherical coordinates:

def make_spot_3d_spherical(bright, spread, rho, theta, phi):


    x0 = int(rho*np.sin(theta)*np.cos(phi))
    y0 = int(rho*np.sin(theta)*np.sin(phi))
    z0 = int(rho*np.cos(phi))

    # Create x and y indices
    x = np.linspace(-50, 50, 200)
    y = np.linspace(-50, 50, 200)
    z = np.linspace(-50, 50, 200)

    X, Y, Z = np.meshgrid(x, y, z)

    Intensity = np.uint16(bright*np.exp(-((X-x0)/spread)**2
                                        -((Y-y0)/spread)**2
                                        -((Z-z0)/spread)**2))

    return Intensity

Two sets of spots are defined (S_t or S_p) by varying theta or phi for a fixed value of rho:

#set of Spots defined by varying theta or phi
S_t = np.asarray([make_spot_3d_spherical(1000,2, 30,t,0) for t in [0,np.pi/6,np.pi/3,2*np.pi/3]])
S_p = np.asarray([make_spot_3d_spherical(1000,2, 30,0,phi) for phi in [0,np.pi/6,np.pi/3,2*np.pi/3]])

Then the set of spots are summed to make a 3D array containing the different spots with np.sum(S_t, axis =0) . Then the 3D array is projected along one of the three axis:

set_of_St0 = np.sum(np.sum(S_t, axis =0), axis = 0)
set_of_St1 = np.sum(np.sum(S_t, axis =0), axis = 1)
set_of_St2 = np.sum(np.sum(S_t, axis =0), axis = 2)

set_of_Sp0 = np.sum(np.sum(S_p, axis =0), axis = 0)
set_of_Sp1 = np.sum(np.sum(S_p, axis =0), axis = 1)
set_of_Sp2 = np.sum(np.sum(S_p, axis =0), axis = 2)

Finally, the different projections are displayed:

plt.subplot(131, xticks=[], yticks=[])
plt.imshow(set_of_Sp0, interpolation = 'nearest')

plt.subplot(132, xticks=[], yticks=[])
plt.imshow(set_of_Sp1, interpolation = 'nearest')

plt.subplot(133, xticks=[], yticks=[])
plt.imshow(set_of_Sp2, interpolation = 'nearest')
plt.show()

In the resulting images,I was waiting to see the spots distributed along some circles, which is not the case, when phi varies:

or when theta varies:

Thanks for advices.

Answer 1

In the function make_spot_3d_spherical you got the sin and cos mixed up in your definition of x0:

x0 = int(rho*np.sin(theta)*np.cos(phi))

should be

x0 = int(rho*np.cos(theta)*np.sin(phi))

Now it works. Here's if you vary phi:

S_p = np.asarray([make_spot_3d_spherical(1000,2, 30,0,phi) for phi in np.linspace(0, 2*np.pi, 20)])
set_of_Sp0 = np.sum(np.sum(S_p, axis =0), axis = 0)
plt.imshow(set_of_Sp0, interpolation = 'nearest')